Concept: Bulk, Skin, Contact, Metal and Pad

In Quokka3 a solar cell is conceptualized into the following regions:

The bulk: it is the main absorber, and (for Si solar cells) can well be assumed to be everywhere quasi-neutral if defined to exclude the near-surface regions. Only the bulk is discretized in the depth coordinate.
The skins: skins are region-wise homogenous areas which cover everything between the quasi-neutral bulk and either the actual surface or the contact to the metal. Typical skins are contacted or non-contacted diffused surfaces, or passivated surfaces (including the inversion or accumulation layer). Quokka3 supports skins to be defined by lumped inputs (lumped skin) e.g. sheet resistance and \(J_0\), or detailed inputs (detailed skin) e.g. doping profile and surface SRH.
The contacts: contacts in Quokka3 define the interface where current can flow between the skin and the metal. Where no contact is defined between the skin and the metal layer, they are assumed to be perfectly isolated.
The metal: represents finger and busbar geometry and accounts for lateral current flow within them (if it is set to be modelled by ‘finite-differences’)
The pads: pads represent the probes or solder pads to which the voltage is applied and current is extracted. They are important only when solving current transport in the metal layer, because otherwise the entire metal geometry effectively represents a probe applying a ‘constant potential’.

*geometric layers of an exemplary PERC ‘busbar enhanced’ unit cell*

The key of Quokka3 to make it numerically efficient and specifically useful for wafer-based silicon solar cells, is to numerically treat the skins always as a parameterized boundary condition when solving 2D or 3D carrier transport in the bulk, which is called the skin concept (Fell et al. 2017). The main lumped parameters to describe a skin as a boundary condition are its recombination property \(J_{0,skin}\), and its sheet resistance \(R_{sheet}\), which are well-established to sufficiently describe of a Si solar cell in many cases. Quokka3 notably can also model skins in detail and then use a more general parameterization, which is applied in the multiscale modeling, ensuring high accuracy also for arbitrary skin properties.

Meshing

Quokka3 uses a structured, orthogonal and non-equidistant mesh to discretize the solution domain and build the finite differences expressions. This means that the solution domain is always strictly cuboidal, and that the geometry definition of the device needs to consist of primitive rectangular geometric features aligned to the coordinate axes. The geometric features are arbitrary in number, position and size, resulting in a generic geometry definition within this rectangular restriction. As an exemption from this a circle as a primitive shape is also supported. It is however recommended in most cases to instead use a square shape with equal area, which requires a much smaller mesh for accurate discretization. The orthogonal mesh is well suited for silicon solar cells, as most cases of interest can be approximated by rectangular device geometries.

This type of mesh comes with the benefits of a rapid and fully automated meshing algorithm (seconds for millions of elements) requiring minimal attention from the user, and of fast and robust electrical solver numerics based on finite differences.

*sketch of solution domain showing structured orthogonal mesh, coordinate system and side name convention*

The 3D solution domain of Quokka3 for above mentioned reasons is strictly cuboidal. To distinguish the boundary planes they are named ‘front’, ‘rear’, ‘east’, ‘west’, ‘north’ and ‘south’ side. For a 2D simulation, the y-coordinate and consequently the ‘north’ and ‘south’ sides are disregarded. For a 1D simulation also the x-coordinate is disregarded. The origin of the coordinate system is defined at the ‘rear’-‘west’-‘south’ corner, which notably means that \(z=0\) corresponds to the rear side, and not to the front surface as in many other tools.

*sketch of solution domain for skins with different depths on the front side*

The bulk mesh covers only the bulk excluding the skin depth, which precisely is the skin’s near-surface layer thickness only. A problem arises if on one side skins with different depths are defined, as the bulk mesh can only have a unique size in z-coordinate. In such a case, Quokka3 calculates the area-average of the skin-depths, and subtracts this average from the user-defined bulk thickness. As far as possible Quokka3 ensures consistency, e.g. still using the correct individual skin depth for optical and detailed electrical modeling. However, for bulk carrier transport, only the single average skin depth can be assumed. This is a good assumption if the differences in skin depths are small compared to the device thickness, but may lead to some inaccuracy otherwise.

Optical modeling

Optical modeling in Quokka3 can be separated into two cases: modeling optics for determining the generation rate in the device, and modelling luminescence and reabsorption to determine the emitted luminescence spectrum and intensity, as well as photon-recycling effect.

Solving generation: the Text-Z model

Besides directly setting the generation via user-defined generation rate or profile, which is called the defined-generation model type, Quokka3 supports spectrally resolved determination of the generation profile via the Text-Z model. It’s a further simplification of the Basore-model (Rand and Basore 1991; Basore 1993), requiring only the (spectrally dependent) input parameters external front surface transmission \(T_{ext}\) and pathlength enhancement factor \(Z\) which quantifies the light-trapping performance of the device. The determination of the generation profile via \(T_{ext}\) and \(Z\) is described in detail in (Fell, McIntosh, and Fong 2016).

The main features of the Text-Z model are:

Being (semi-) analytical, it is rapid for all purposes
The determined generation profiles are accurate for non-extreme wafer-based silicon solar cells
The input parameters can be derived either from experiments or detailed optical simulations (e.g. ray tracing + transfer matrix), providing the most suitable interface to detailed optical modelling tools
Being spectrally resolved, it enables rapid optical modeling for quantum efficiency (QE) simulations
It is compatible with the skin concept, i.e. accounts for generation and collection efficiency in the skin, and also for active top cell generation within a tandem cell simulation.
The input parameters are, to good approximation, independent of temperature and the wafer thickness (the latter if additionally a parameterization of \(Z\) is used (McIntosh and Baker-Finch 2015)), meaning that those quantities can be varied within Quokka3 without re-doing detailed optical modeling or characterization.

Relation to spectral response (SR)

*exemplary spectral response (SR) of a silicon solar cell; the black line divides external and internal optics*

To understand the spectral optical modelling approach of Quokka3, it is useful to distinguish between external and internal optics, which are separated by the black line in the figure above. External optics means everything which happens until the incident photons reach the active absorber material of the solar cell, namely reflection from the surface, i.e. the external reflection \(R_{ext}\), and parasitic absorption in the skin on the illuminated side, e.g. short-wavelength absorption within an anti-reflection coating \(A_{ext}\). All effects happening after this are called internal optics. The internal optics is dominated by active absorption in the bulk material up to a wavelength where the absorption becomes moderate relative to the device thickness. Then light is able to reach the opposite non-illuminated side of the device, for typical c-Si this starts around 950 nm. Between this wavelength and the band-gap (~1200 nm for c-Si), light is reflected multiple times between the surfaces. The ability to actively absorb light in this near-infrared (NIR) wavelength regime is called light-trapping. Light-trapping is reduced by parasitic absorption both in the bulk material (e.g. FCA) and also by NIR absorption in the surface films, i.e. within the skins. Light-trapping can never be 100% effective as at least at the illumated side NIR light can escape the device, which is called escape reflection \(R{esc}\). The same goes for transmission at the rear side of a bifacial cell. Theoretical limits for the light-trapping capabilities have been proposed, see e.g. (Green 2002).

Within the Text-Z model, \(T_{ext}\) lumps all external effects, which is defined as the fraction of incident photons neither reflected nor absorbed before reaching the active absorber: \[ T_{ext}=1-R_{ext}-A_{ext} \] If a subcell is present at the illuminated side within a tandem simulation, then \(T_{ext}\) is further reduced by the active absorption in the top cell \(A_{gen}\): \[ T_{ext}=1-R_{ext}-A_{ext}-A_{gen} \] The active absorption \(A_G\) of transmitted photons in the light-trapping regime is fully quantified by \(Z\), the absorption coefficient of the absorber \(\alpha\) and its width by virtually extending the width to match \(A_G\) with using the Lambert-Beer law: \[ A_G(\lambda)=1-\exp\left(-\alpha(\lambda)Z(\lambda)W\right) \] optical parameters of an example silicon solar cell

\(T{ext}\) and \(Z\) fully determine the optical influence on the measurable external quantum efficiency (EQE) of the device. The electrical influence on EQE is quantified by the collection efficiency \(\eta_c\), which is the fraction of collected over generated charge carriers actually contributing to the device current. The input parameter \(Z\) is directly related to the absorbed fraction of transmitted photons, \(A_G\), and the internal quantum efficiency (IQE). \[ EQE(\lambda)=T_{ext}(\lambda)IQE(\lambda) \]

\[ IQE(\lambda)=A_G(\lambda)\eta_c(\lambda) \]

Notably, for a given spectrum and device thickness \(W\), the device’s generation current \(J_G\) is fully defined by the Text-Z model’s input parameters. Consequently, this most important quantity is ensured to be consistent, e.g. when the inputs are derived from ray tracing software. \[ J_G=q\int_{0}^{\infty}{I(\lambda)T_{ext}(\lambda)\left(1-\exp\left(-\alpha(\lambda)Z(\lambda)W\right)\right)d\lambda} \]

Calculating the generation profile

*sketch for the Text-Z model, a) assumed (textured) geometry, b) quasi-1D representation in Quokka3*

The Text-Z model is based on the faceted surface model (Basore 1993). It accounts for different effective path angles in the near-surface region (i.e. in the skin) and in the bulk. This way, both the short-wavelength and long-wavelength IQE is correctly modelled. Differently to the Basore model, light reflected at the rear-surface is not followed further by internal optical properties of the surfaces. Rather, the remaining generation after the first path \(l_1\) is distributed equally over the device thickness: \[ G(\zeta,\lambda)=G_1(\zeta,\lambda)+W^{-1}I(\lambda)T_{ext}(\lambda)\left(\exp\left(-\alpha(\lambda)l_1\right)-\exp\left(-\alpha(\lambda)Z(\lambda)W\right)\right) \]

Different ways to define light-trapping / \(Z\)

Quokka3 supports various options to define or model light-trapping:

Directly user-defined, i.e. the user provides \(Z(\alpha)\) values as a function of absorption coefficient. Such \(Z\) values can be derived e.g. via a ray-tracing software, or via other optical modeling external to Quokka3.
Choosing a limit for light-trapping
- \(4n^2\): well-known limit where \(Z(\lambda)=4(n(\lambda))^2\); notably it assumes zero absorption, and is thus valid for very weakly absorbing wavelengths only.
- (Green 2002): Green did improve over the \(4n^2\) limit by considering the effect of absorption, thus providing a more realistic upper limit of light trapping in a device: \(Z(\lambda) = 4 + \frac{\ln\left((n(\lambda))^2 + \left(1 - (n(\lambda))^2\right)\exp(-4\alpha (\lambda)W)\right)}{\alpha (\lambda)W}\)
A parameterization of \(Z\) into three parameters \(Z_\infty\), \(Z_0\) and \(Z_p\) (McIntosh and Baker-Finch 2015)
Using the extended Basore model, an extension of PC1D’s light-trapping model (Basore 1993). It models light trapping using lumped internal optical properties of the skins, which is described in detail here. This option enables to define a spatial variation in the light trapping properties of the device by assigning different internal optical properties to the skins.

Z parameterization

A parameterization to describe the dependence of \(Z\) on the wavelength, or more fundamentally on the absorption-coefficient, was derived in (McIntosh and Baker-Finch 2015): \[ Z(\alpha)=Z_\infty+\frac{1}{\alpha Z_pW}\ln\left(\frac{Z_0}{Z_\infty}-\left(\frac{Z_0}{Z_\infty}-1\right)\exp\left(-\alpha Z_\infty Z_pW\right)\right) \] The benefits arising from the \(Z\) parameterization are the following:

Light trapping is accurately described by the 3 parameters \(Z_0\), \(Z_\infty\) and \(Z_p\) only, instead of requiring a full spectrally dependent curve.
The parameters are specific to the optical properties of the surfaces only, i.e. independent of the device thickness \(W\). As also \(T_{ext}\) is independent of \(W\), this means that \(W\) can be varied in Quokka3 without adjusting optical inputs, yielding correct changes in generation.
The parameters are, to good approximation, independent of temperature, allowing variation of temperature without adjusting optical inputs. See next section for details.

Temperature (in-)dependence

In (Fell, McIntosh, and Fong 2016) it was shown that both \(T_{ext}\) and \(Z\) can be approximately assumed temperature independent for typical silicon solar cell operating conditions. When using Quokka3’s Text-Z model it is therefore conveniently possible to vary the temperature without having to adjust any optical inputs, and get a correct prediction of the change in generation current.

This is due to the fact that temperature dependence of the optical properties is largely dominated by changes in the absorption coefficient \(\alpha(\lambda)\), which is a different function of wavelength for different temperatures. As mentioned in the Z parameterization, \(Z\) fundamentally is a function of \(\alpha\), which is the only wavelength- and temperature-dependent quantity. Therefore, the input \(Z\) needs to be defined as a function of \(\alpha\), or by its parameterization, and can then be considered independent of temperature. The reasoning is similar for the long-wavelength region of \(T_{ext}\). For the short-wavelength region of \(T_{ext}\) however, the wavelength-dependence can be considered temperature independent instead. Often for typical silicon solar cells the entire \(T_{ext}(\lambda)\) curve can well be considered temperature-independent.

Modeling light-trapping and luminescence: the extended Basore model

Model description

sketch for the extended Basore model, a) assumed (textured) geometry, b) quasi-1D representation in Quokka3; the model inputs are given in blue (and green / red for Appp), the outputs are given in grey

Note that for simplicity the following description assumes front side illumination, meaning that “front” means the illuminated side, and “back” the non-illuminated side. However, the model can equivalently also be used for back side illumination, and even bifacial illumination (Fell et al. 2019).

The extended Basore model is based on (Basore 1993). It is an analytical model describing light-trapping, i.e. the internal optics, of a solar cell based on the main absorber properties (thickness, absorption coefficient) and internal optical properties of the surfaces, mainly internal reflection. It is able to calculate the spectrally resolved probability of absorptance \(A\), escape reflection \(R_{esc}\) and transmission \(T\) of photons which did reach the main absorber, i.e. have not yet been lost by external optical effects (external reflection and first-pass absorption in front side films).

The required internal optical input parameters, which in Quokka3 are defined as the internal optical properties of the skins, are:

the first and \(n^{th}\)-pass internal reflection of the front and rear skin \(R_{f1}\), \(R_{fn}\), \(R_{b1}\) and \(R_{bn}\); note that often it is useful to simplify the first and \(n^{th}\) pass reflection to have the same value, e.g. for calibrating the model via fitting to reflectance and transmittance measurements (Fell et al. 2019).
The absorbed fraction per perpendicular pass \(A_{ppp}\) representing parasitic absorption in the skins; this is an extension to the original model of Basore, see next section.
whether light is reflected specular or diffuse, the latter assuming an ideal Lambertian scatterer; while it makes a significant difference whether the sample / cell is assumed full planar (both sides specular) or not, it makes only little difference whether one or both sides are set to diffuse.
the facet angle of the front side texture, which determines the angle of the first pass through the sample / cell.

Additionally one needs to input the external optical properties of the front side, namely the external reflectance \(R_{ext}\) and external absorptance \(A_{ext}\), which relate to the external transmission \(T_{ext}\):_ \[ T_{ext}=1-R_{ext}-A_{ext} \] For each wavelength, the model analytically traces a representative ray along the multiple passes and bounces within the sample / cell. It discriminates the first, second and \(n^{th}\) pass of the through the thickness \(W\) with the angle \(\theta_1\), \(\theta_2\) and \(\theta_n\), respectively. The angles determine the transmission of each pass with the total absorption coefficient \(\alpha\): \[ T_{i=1,2,n}=\exp\left(\frac{-\alpha W}{\cos(\theta_i)}\right) \] \(\theta_1\) is determined using Snell’s law and a facet angle representing the surface texture, \(\theta_2\) and \(\theta_n\) take on the same value after reflection from a specular back side. In case any side is defined to be diffusive, \(T_n\), and subsequently \(\theta_n\), is determined by assuming a Lambertian angular distribution resulting in (Green 2002): \[ T_n=\exp(-\alpha W)(1-\alpha W)-(\alpha W)^2Ei(-\alpha W) \] with \(Ei(x)\) being the exponential integral function. \(\Theta_2\) and \(T_2\) take the same values as \(\Theta_n\) and \(T_n\)if the back side is diffusive.

Using a geometric sum for the infinite number of \(n^{th}\) passes one can derive the following formulas for the spectral absorbance \(A\), transmittance \(T\) and escape reflectance \(R_{esc}\): \[ R_{esc}=T_{ext}\left[1-R_{b1}T_1T_2(1-R_{f1})-\frac{R_{b1}R_{f1}T_1T_2R_{bn}(1-R_{fn})T_n^2}{1-R_{bn}R_{fn}T_n^2}\right] \]

\[ T=T_{ext}\left[T_1(1-R_{b1})+\frac{R_{b1}R_{f1}T_1T_2T_n}{1-R_{bn}R_{fn}T_n^2}(1-R_{bn})\right] \]

\[ A=T_{ext}\left[1-T_1+R_{b1}T_1(1-T_2)+\frac{R_{b1}R_{f1}T_1T_2R_{bn}(1-T_n)(1+R_{bn}T_n)}{1-R_{bn}R_{fn}T_n^2}\right] \]

The total absorbance in the main absorber can be further broken down into active (=generation) \(A_{gen}\) and parasitic absorption, namely free carrier absorption FCA \(A_{FCA}\), where the total absorption coefficient equals the sum of the band-to-band and FCA absorption coefficient, \(alpha_{bb}\) and \(alpha_{FCA}\), respectively: \[ A_{gen}=A\frac{\alpha_{bb}}{\alpha_{bb}+\alpha_{FCA}}\\ A_{FCA}=\frac{\alpha_{FCA}}{\alpha_{bb}+\alpha_{FCA}} \]

Accounting for skin parasitic absorption via Appp

As an extension to the original Basore model, parasitic absorption in the skin can be accounted for by the parameter \(A_{ppp}\), which is the absorbed fraction per perpendicular pass. It is described in detail in (Andreas Fell, Greulich, et al. 2022).

Essentially, at each bounce the intensity of the representative ray is not only reduced by the internal reflection \(1-R\), but by considering the pass through the absorptive skin to the reflective surface and back again with the angle \(\theta\): \(A_{ppp}\cos(\theta)(1-R)A_{ppp}\cos(\theta)\).

The main application of \(A_{ppp}\) is to represent parasitic absorption in the light-trapping regime, in contrast to the first-pass parasitic absorption losses within a skin which is represented by \(A_{ext}\).

One should be careful on the correct differentiation when the skin is on an illuminated side: the same absorption mechanism can lead to a non-negligible first-pass, i.e. external absorption loss, but also to a light-trapping loss resulting from multiple internal bounces at that skin. This is for examples the case for a front TCO layer, which shows significant absorption in the NIR, meaning in the light-trapping regime. It can still be useful to discriminate this absorption into it’s external (\(A_{ext}\)) and internal (\(A_{ppp}\)) part in order to quantify the light-trapping losses as explained in (Andreas Fell, Greulich, et al. 2022) and on the corresponding HJT optics example on this website.

Luminescence, photon reabsorption and photon recycling

The extended Basore model with the identical inputs used for modeling light-trapping for external illumination is also used in Quokka3 to model photon escape and reabsorption, which impacts luminescence and photon recycling, respectively. The details are described in (Fell et al. 2021).

After the electrical solution Quokka3 knows the Fermi level split at each mesh element, which makes it straightforward to calculate the radiative recombination (=spontaneous emission) rate \(R_{rad}\) at each element: \[ R_{rad}=B_{rad,low}n_{i,0}^2\exp\left(\frac{\Delta\eta}{kT}\right) \] where \(B_{rad,low}n_{i,0}^2\) was shown to be invariable for carrier densities typically occurring in the silicon bulk in (Fell et al. 2021), and was parameterized in (Nguyen, Baker-Finch, and Macdonald 2014). It is thus accurate without the need to account for injection-dependence via a \(B_{rel}\) or band-gap-narrowing / \(n_{i,eff}\) model.

For spectrally resolved emission, as required for spectral luminesce modeling, Quokka3 rather uses the equivalent generalized Planck law to calculate the spontaneous emission rate \(r_{spont}\) for each wavelength from the band-to-band absorption coefficient: \[ r_{spont}(\lambda)=\frac{(\hbar\omega)^2n^2\alpha_{bb}}{\pi^2\hbar^3c_0^2}\exp\left(\frac{\Delta\eta}{kT}\right) \] The emitted photons can then either escape at the detection side in form of the luminescence signal, or be reabsorbed, the respective probabilities being \(f_{esc}\) and \(f_{reabs}\). The band-to-band part of the reabsorption leads to additional electrical generation reducing the effective radiative recombination rate, which is called photon-recycling.

*sketch for using the extended Basore model for left: escape probability, right: reabsorption probability*

For the case of escape probability, i.e. for modelling the luminescence signal, it has been shown that Basore’s light-trapping model can be used via the reciprocity of absorption and emission, see for example (Schinke et al. 2013). It assumes a small solid angle of detection. If applied to a 3D solution, it notably results in a hyperspectral image of the detectable luminescence signal, which in Quokka3 is converted into the luminescence image by applying the spectral transmission of the detection optics and also the spectral sensitivity of the camera, both being definable by the user. The formula for the escape probability can be derived to: \[ f_{esc}(\lambda)=\frac{\Omega}{4\pi n^2}T_{ext}\left[\frac{1}{\cos\theta_1}\exp\left(\frac{-\alpha z}{\cos\theta_1}\right)+\frac{T_1R_{b1}}{\cos\theta_2}\exp\left(\frac{-\alpha (W-z)}{\cos\theta_2}\right)+\frac{T_1R_{b1}T_2R_{f1}}{(1-T_n^2R_{bn}R_{fn})\cos\theta_n}\left(\exp\left(\frac{-\alpha z}{\cos\theta_n}\right)+T_nR_{bn}\exp\left(\frac{-\alpha (W-z)}{\cos\theta_n}\right)\right)\right] \]

For the case of reabsorption, and also if at least one side is diffuse / scattering, the reabsorption process is well represented by the \(n^{th}\) pass of the light-trapping model only. This is because the radiative emission is already isotropic, and the scattering at the surface leads to approximately the same randomized distribution of path angles for all passes, which situation is represented by the \(n^{th}\) pass in the light-trapping model. The resulting reabsorption probability is: \[ f_{reabs}(\lambda)=1-\frac{1-T_n}{\frac{\alpha W}{\cos\theta_n}}\left(1-\frac{(1-T_n)\frac{R_{bn}+R_{fn}}{2}}{1-T_n\frac{R_{bn}+R_{fn}}{2}}\right) \] Quokka3 supports an ideal option for internal optics, which can be applied to model both luminescence and photon recycling. The ideal options means that the front and rear is assumed to be diffusive, and the the internal reflections are set to their Lambertian limit: \(R=1-1/n^2\) (Green 2002).

Note that the model for the planar case described in (Fell et al. 2021) is not yet available in Quokka3.

To derive the absolute emission and reabsorption rates, the probabilities are multiplied with the spontaneous emission rate and integrated over the wavelength, for the luminescence case multiplying further with the optics transmission and sensor sensitivity before the integration.

If injection dependence of the reabsorption probability can be neglected, photon recycling can simply be accounted for multiplying the radiative recombination rate with one minus the probability of reabsorption for all wavelengths, which is determined by: \[ f_{reabs}=\frac{\int f_{reabs}(\lambda)\frac{\alpha_{bb}}{\alpha_{bb}+\alpha_{FCA}}r_{spont}(\lambda)d\lambda}{\int r_{spont}(\lambda)d\lambda} \] Evaluating the above equation only once in low-injection conditions can become inaccurate if the reabsorption itself is dependent on carrier densities. This can happen in the case for very high lifetime samples showing a substantial fraction of radiative recombination of the total recombination, and if high injection conditions are investigated so that FCA of excess carriers becomes relevant, see (???). This would require iterative coupling of the electrical and optical modeling to re-evaluate above equation, which is not yet implemented in Quokka3.

Note that even though not shown in the equations above, parasitic skin absorption via \(A_{ppp}\) is also accounted for in the escape and reabsorption modeling within Quokka3.

Inhomogeneous properties

The extended Basore model is the only option to account for optical inhomogeneities of the cell in the light-trapping regime. Similar to the external optics properties, i.e. \(T_{ext}\), the internal optical properties can be set for each skin individually. However, differently to the external optics, before evaluating the extended Basore model, all internal input parameters are area-averaged to derive a single quantity for each input on each side. For the case of the back-side reflection \(R_{bn}\) this is has been shown to be accurate in (Sabater et al. 2021), opposed to area-averaging the results of a separate Basore model evaluation for each region. The same reasoning is applied to the other internal optical properties, therefore Quokka3 area-averages all inputs before a single Basore model evaluation.

Notably this area-averaging is valid if the scale of the inhomogeneities is within the lateral length scale of the light-trapping, which is well the case for e.g. local contacts. In case there are large-scale systematic inhomogeneities within a solution domain, it must be kept in mind that this inhomogeneity is removed via the area-averaging for the light-trapping modeling, and thus also luminescence modeling.

Electrical modeling

General semiconductor model equations

Current transport in a semiconductor is well established to be modelled by coupling the continuity equations for charge carrier transport with the Poisson equation, also known as the drift-diffusion model. In Quokka3, up to 4 charge carriers are considered in the 1D detailed solver, that is anions and cations additionally to electrons and holes. The following equations are given in the most general form to cover all effects considered in Quokka3. Note however that those equations are never solved in such generality, meaning a fully coupled transient solution of continuity and Poisson equations in 3D, including ion transport. Refer to the individual solver section to see under what specific assumptions those equations are solved.

Volume carrier transport

Electron and hole carrier transport is described by the respective continuity equation: \[ \frac{dn}{dt}=q(G-R_{el})-\nabla\vec{J}_{el} \]

\[ \frac{dp}{dt}=q(G-R_{hol})+\nabla\vec{J}_{hol} \]

Quokka uses the electrochemical or quasi-Fermi potentials \(\phi_{el}\) and \(\phi_{hol}\) as solution variables, from which the current density is calculated as: \[ \vec{J}_{el/hol}=\kappa_{el/hol}\nabla\phi_{el/hol} \] The conductivities \(\kappa_{el}\) / \(\kappa_{hol}\) are calculated from the respective mobility \(\mu_{el}\) / \(\mu_{hol}\) and carrier density: \[ \kappa_{el} = q\,n\,\mu_{el}, \kappa_{hol}=q\,p\,\mu_{hol} \] Ion transport can be modeled equivalently, with anions and cations as charge carriers, however without generation and recombination: \[ \frac{dc_{an}}{dt}=-\nabla\vec{J}_{an} \]

\[ \frac{dc_{cat}}{dt}=\nabla\vec{J}_{cat} \]

\[ \vec{J}_{an/cat}=\kappa_{an/cat}\nabla\phi_{an/cat} \]

The Poisson equation relates the electric field \(\vec{E}\) to the net charge density: \[ \nabla\varepsilon_r\varepsilon_0\vec{E}=q(N_D-N_A+p-n+c_{cat}-c_{an}-n_t) \] The electric field is the negative gradient of the electrostatic potential \(\varphi_e\): \[ \vec{E}=-\nabla\varphi_e \]

Carrier statistics

Electrons and holes

Carrier statistics relate the semiconductor band-structure (i.e. energy levels) to the carrier densities. Generally valid is the use of Fermi-Dirac statistics: \[ n=N_c\,F_{1/2}\left(\frac{\phi_{el}-E_c/q}{V_t}\right) \]

\[ p=N_v\,F_{1/2}\left(-\frac{\phi_{hol}-E_v/q}{V_t}\right) \]

Here \(F_{1/2}()\) denotes the half-order Fermi integral. The conduction and valence band edges can be calculated from the band gap and the electrostatic potential: \[ E_c=-q\,\varphi_e+E_{g,0}/2-\delta E_{BGN,c} \]

\[ E_v=-q\,\varphi_e-E_{g,0}/2+\delta E_{BGN,v} \]

An important quantity is the intrinsic carrier density \(n_i\), which is related to the density of states and the band gap: \[ n_{i,0}^2=n_0p_0=N_cN_v\exp\left(-\frac{E_{g,0}}{q\,V_t}\right) \] If band gap narrowing (BGN) becomes significant, the effective intrinsic carrier density is defined by: \[ n_{i,eff}^2=n_{i,0}^2\exp\left(\frac{\delta E_{BGN,c}+\delta E_{BGN,v}}{q\,V_t}\right) \]

If carrier densities (i.e. doping densities) are below \(10^{19} cm^{-3}\), Fermi-Dirac statistics can be simplified to the mathematically easier Boltzmann statistics, where the Fermi integral is replaced by an exponential function (Altermatt 2011). A well known expression derived from Boltzmann statistics is that the product of the carrier densities (pn-product) is directly related to the Fermi level split: \[ pn=n_{i,eff}^2\exp{\left(\frac{\phi_{el}-\phi_{hol}}{V_t}\right)} \]

Ions

For ions, Quokka3 uses Boltzmann statistics: \[ c_{an}=N_{0,an}\,\exp\left(\frac{\phi_{an}+\varphi_e}{V_t}\right) \]

\[ c_{cat}=N_{0,cat}\,\exp\left(-\frac{\phi_{cat}+\varphi_e}{V_t}\right) \]

Generation and recombination

The generation rate of electron-hole pairs \(G\) is the result of optical modeling.

Recombination of electrons and holes happens via 3 different mechanisms:

Radiative recombination, \(R_{rad}\)
Auger recombination, \(R_{aug}\)
Defect recombination via Shockley-Read-Hall (SRH) formalism, \(R_{SRH}\), as described in (Shockley and Read Jr 1952)

The first two are unavoidable contributions due to the intrinsic material properties, and are thus called intrinsic recombination, \(R_{intrinsic}=R_{rad}+R_{aug}\). They are assumed fast enough to be considered (quasi-)steady-state, and the recombination rates are equal for electrons and holes.

General SRH

This section describes the general SRH formalism described in (Shockley and Read Jr 1952). General SRH can currently be enabled for a single bulk defect within the 1D-detailed solver only. The default way to model SRH for all other cases is by its simplified form.

SRH recombination happens via defects having their energy level \(E_t\) placed within the band gap, called traps in the general formalism. In Quokka3, \(E_t\) is defined as the value relative to the middle of the band gap. The total defect / trap density \(N_t\) is differentiated from the density of occupied traps \(n_t\), and its fraction is defined as \(f_t=\frac{n_t}{N_t}\). The recombination rates for electrons and holes are their net rates of capture and emission, and can become negative if emission is larger: \[ R_{SRH,el}=\frac{\left(1-f_t\right)(n-n_0)-(n_0+n_1)\left(f_t-f_{t0}\right)}{\tau_{n0}} \]

\[ R_{SRH,hol}=\frac{f_t(p-p_0)+(p_0+p_1)\left(f_t-f_{t0}\right)}{\tau_{p0}} \]

The equilibrium occupied fraction is defined as: \[ f_{t0}=\frac{1}{1+\frac{n_1}{n_0}} \] \(n_1\) and \(p_1\) are the carrier concentrations for the case that the respective Fermi level coincides with the defect level: \[ n_1=n_{i,eff}\exp\left(\frac{E_t}{q\,V_t}\right) \]

\[ p_1=n_{i,eff}\exp\left(-\frac{E_t}{q\,V_t}\right) \]

The fundamental electron and hole lifetimes \(\tau_{n0}\) and \(\tau_{p0}\) can be calculated from the capture cross sections \(\sigma_n\) and \(\sigma_p\) as: \[ \tau_{n0/p0}=\frac{1}{v_{th}\sigma_{n/p}N_t} \] (Only) for the general formalism, the continuity equation for the trap density needs to be solved: \[ \frac{dn_t}{dt}=R_{SRH,el}-R_{SRH,hol} \]

Simplified SRH

Commonly in silicon solar cell modeling the simplified form of the general SRH formalism is used, which assumes (quasi-)steady-state and that the occupied trap density is small compared to the excess carrier density, i.e. \(\frac{n-n_0}{p-p_0}\simeq1\), see discussion in (Macdonald and Cuevas 2003). This means essentially that the effect of trapping is negligible. Then the electron and hole recombination rates are equal and become the well-know expression: \[ R_{SRH,el}=R_{SRH,hol}=R_{SRH}=\frac{p\,n-n_{i,eff}^2}{\tau_{n0}(p_1+p)+\tau_{p0}(n_1+n)} \]

Solving the bulk with lumped skins: qn-bulk solver

The qn-bulk solver solves 3D carrier transport in the bulk, as well as 2D current transport in up to two layers on each side: the skin layer and the metal layer. The following assumptions are made, which are specifically scoped to be valid for typical silicon solar cells:

bulk semiconductor transport with lumped skins only: 3D semiconductor carrier transport is only modelled in the bulk. Anything happening in the skins is “lumped away” into boundary conditions, including the lateral current transport. Note that the user can still provide detailed inputs for the skins, in which case multiscale modeling is used.
quasi-neutrality: as the space-charge regions are “outsourced” into the skins, the bulk region where semiconductor carrier transport is modelled can well be assumed to be quasi-neutral. This should NOT be confused with low-injection. It also does NOT set charge density within the Poisson equation to zero meaning zero electric field. Rather solving the Poisson equation can be omitted, and from the solution of the qn-bulk solver the correct electric field can be derived. Quasi-neutrality in the bulk holds very well for any typical silicon solar cell condition.
Boltzmann statistics: as within the bulk both doping and excess carrier densities are well below the critical values, Boltzmann carrier statistics can be used without loss of accuracy.
steady-state only
simplified SRH only, no ion transport

Volume equations

The main benefit from the quasi-neutrality assumption is that the carrier densities can be calculated without knowledge of the electrostatic potential. This way, the continuity equations can be solved without solving the Poisson equation, which greatly reduces numerical challenges and thus increases speed and robustness. That means NOT that quasi-neutrality is assumed for the Poisson equation, which would render the electric field to zero. Rather the electric field causing the drift currents is implicitly solved by solving the continuity equations, which is accurate if the quasi-neutrality assumption holds. \[ \nabla\vec{J}_{el}=q(G-R) \]

\[ \nabla\vec{J}_{hol}=-q(G-R) \]

Applying quasi-neutrality and Boltzmann statistics, the carrier densities can be calculated as: \[ n-n_0=p-p_0=\sqrt{\frac{(n_o+p_0)^2}{4}-\left(p_0n_0-n_{i,eff}^2\exp{\frac{\phi_{el}-\phi_{hol}}{V_t}}\right)}-\frac{n_0+p_0}{2} \] The electric field can be post-calculated once the solution for the Fermi levels has been found, either by the electron or hole density: \[ \varphi_e=-\left(\phi_{el}-V_t\ln{\left(\frac{n}{n_{i,eff}}\right)}\right)=-\left(\phi_{hol}+V_t\ln{\left(\frac{p}{n_{i,eff}}\right)}\right) \]

Boundary conditions

Symmetry / perfect edge boundary

By default, that is if no skin layer is defined on a specific side of the solution domain, a symmetry boundary condition is applied. Notably, the symmetry condition is equivalent to a perfect edge of the device, as it represents no current through and no recombination at this edge. \[ \vec{n}\vec{J}_{el/hol}=0 \]

Lumped skin and metal layer boundary

If skins are defined on a side of the solution domain, a single-element-deep mesh layer is added to the bulk’s mesh. Within that layer, lateral current transport is modelled by solving the continuity equation for the skin’s potential \(\varphi_{skin}\), where the conductivity is defined by the skin’s sheet resistance \(R_{sheet}\): \[ \nabla_t\left(\frac{1}{R_{sheet}}\nabla\varphi_{skin}\right)=\left(\vec{n}\vec{J}_{el}+\vec{n}\vec{J}_{hol}\right)-J_{cont} \] The right hand side is a source term, with \(\vec{n}\vec{J}_{el}+\vec{n}\vec{J}_{hol}\) being the net current density from the bulk into the skin, and \(J_{cont}\) the current density from the skin into the metal (i.e. through the contact).

Optionally on the front and rear side, an additional mesh layer can be added to solve current transport in the metal, which is called the ‘finite differences’ metal model type. The equation to solve is equivalent as for the skin layer, using the metal potential \(\varphi_{met}\), the metal sheet resistance, and \(J_{cont}-J_{pad}\) as the source term, where \(J_{pad}\) is the current from the metal into the pad.

The current density between the skin and metal layer, as well as between the metal and pad layer, is calculated via the respective contact resistivity: \[ J_{cont}=\frac{\varphi_{met}-\varphi_{skin}}{\rho_{cont}} \]

\[ J_{pad}=\frac{\varphi_{pad}-\varphi_{met}}{\rho_{pad}} \]

If current transport in the metal layer is modelled by ‘finite differences’, the pad potential \(\varphi_{pad}\) equals the applied voltage, which is defined to be zero for a ‘p-type’ polarity, and \(V_{int}\) for a ‘n-type’ polarity pad. If current transport in the metal is not modelled, the metal is assumed to be at ‘constant potential’, meaning that the metal potential is set to the applied voltage with the respective polarity, and the pads are consequently disregarded. The ‘constant potential’ metal model type should be the default choice when modeling common unit cells, as metal current transport within such a small domain size is usually not significant.

The current densities between the bulk and the skin form the actual boundary conditions for solving the bulk continuity equations. They, together with the lateral current transport, represent all physics happing within the skin in an effective way. As discussed in (Fell et al. 2017), for a steady-state solution, a single skin element can be considered a black-box electrical circuit element with a unique operating point defined by the potential differences between the skin and the neighboring elements. If the relation between the current densities \(\vec{n}\vec{J}_{el}\), \(\vec{n}\vec{J}_{hol}\) and \(J_{cont}\) and the potential differences is accurately known for all relevant operating points, the influence of the skin on the device characteristics is exactly described.

To parameterize the mentioned relation, an equivalent circuit model of a skin element is used. The parameterization can represent all possible steady-state current-potential relationships, and is thus generally valid for any properties of typical solar cell skins. The only relevant assumptions are taken for the case that lateral current transport in the skin is significant compared to lateral transport in the bulk: i) lateral current transport can be described via assuming a single skin potential, ii) it can be superimposed to vertical current transport, and iii) it is significant for one type of carrier only (named the majority carrier), i.e. the skin is either ‘p-type’ or ‘n-type’. Those assumptions hold well for most typical silicon solar cell skins, essentially because they are thin. The relevant potential differences reduce in this case to the bulk-side quasi-Fermi level split \(\phi_{split}=\phi_{el}-\phi_{hol}\), and the potential drop of majority carriers between the bulk and the skin layer \(\Delta\varphi_{maj}=\varphi_{skin}-\phi_{el/hol}\).

*single-element-deep skin layer (left), general parameterization model into boundary condition (right)*

The parameters of the lumped skin model are the sheet resistance \(R_{sheet}\), the recombination property \(J_{0,skin}\), a generation current density \(J_{G,skin}\), the collection efficiency \(\eta_c\), and a vertical resistivity \(\rho_{skin}\). It also includes the contact to the metal via the contact resistivity \(\rho_{cont}\). The parameters have a good physical meaning in most cases. However they generally are mathematical parameters which, with full potential dependence, can represent any operating point of a typical skin. For example the vertical resistivity is a meaningful parameter to represent resistance over a thin oxide within a polysilicon passivating contact, however it becomes strongly non-ohmic and complex to understand when it represents the resistance over a pn-junction’s space-charge-region within the full multiscale coupling.

When the vertical resistivity is negligible and the other parameters are assumed to be constant, the general skin parameterization is equivalent to the conductive boundary model, published in (Brendel 2012; Fell 2013) and used by Quokka2.

Quokka3 supports to define the lumped parameters directly as user inputs, with constant vertical resistance and sheet resistance, but allowing injection-dependence of the recombination, i.e. \(J_{0,skin}(\phi_{split})\). Alternatively, the skin solver can parameterize the results from a detailed solution of the skin and automatically apply it as the boundary condition to the qn-bulk solver, which is called multiscale modeling.

In addition to the default ohmic values for the lumped vertical skin resistivity and the contact resistivity between skin and metal, Quokka3 also supports a lumped MIS model to define a non-ohmic, i.e. diode-like contact resistivity. The diode parameter \(J_0\) is thereby calculated via a thermionic emission plus direct tunneling model, neglecting the potential drop over the insulator, for the details see TLS / MIS contact: \[ J_{0,MIS}=A_0T^2 \exp\left(-\frac{\Phi_{B,el/hol}}{V_t}\right)P_{tun} \] Together with a parallel ohmic resistivity the current density is calculated to: \[ J=J_{0,MIS}\exp(\frac{\Delta\varphi}{V_t})+\rho_c\Delta\varphi \] The lumped MIS model might be a useful approach for passivating contacts for the typical case that the non-ohmic resistance effect should be accounted for, however the detailed interface properties are not known precisely and thus modeling recombination by a lumped approach is preferable.

The lumped MIS model is also compatible with the resistive mode described below, enabling to simulate resistance structures with MIS-type contacts.

For tandem simulations, the equivalent-circtiut model defining the lumped skin boundary condition is extended to include a subcell JV-curve, see tandem modeling.

Resistive mode

The qn-bulk solver also supports an resistive mode which is used when solving a resistive device, for which the following additional assumptions are taken:

Solves a single potential only: any conductivity-type setting (‘n-type’ or ‘p-type’) is disregarded but assumed to be the same type for bulk and all skin features.
No excitation: all conductivities and resistivities are constant values, meaning also the resulting device will be ohmic, i.e. having a constant resistance regardless of the applied voltage. No semiconductor physics other than a single continuity equation per element are solved.
MIS contact exception: the only non-ohmic feature within the resistive mode are lumped MIS resistivities between the bulk and the skin, or between the skin and the metal, see previous section.

Equivalent-circuit bulk model type

As an alternative to numerically solve the current transport equations within a 3D bulk mesh, Quokka3 also supports a so-called equivalent-circuit model type to account for generation, recombination and lateral carrier transport in the bulk. Here, the bulk is not discretized at all, meaning the numerical mesh only consists of the 2D sheet meshes in the skins and optionally the metal layers. The bulk is represented via a parallel diode and current source connection of the front and skin layer, and also by contributing to the lateral sheet resistances.

*left: 3D drift-diffusion qn-bulk solver with quasi-neutrality assumption; right: equivalent circuit representation of the bulk*

The figure above shows one lateral mesh column for left: the quasi-neutral transport equations in the discretized bulk region; and right: replacing the bulk mesh by a diode and current source representing recombination (bulk + surface) and generated current, respectively.

Due to not having to solve for two Fermi-potentials in the bulk mesh, the equivalent-circuit mode strongly simplifies the numerical solution effort. It is thus of particular usefulness for large geometries, if the details of the bulk transport mechanism are of lesser importance. The advantages of the equivalent-circuit mode are:

Much higher simulation speed compared to solving the quasi-neutral current transport equations.
- Enables large geometries like e.g. full cell domains to be solved in practical timeframes.
Better numerical robustness.
- Copes well with small and negative voltages, a drawback of the qn-bulk solver.
It is easy to enable in Quokka3 by a single setting option, all definitions of the device properties stay the same.

These advantages come at a price, namely the lack of detail in carrier transport physics in the bulk. In particular, depth-uniform carrier densities are effectively assumed, meaning that all related effects are neglected:

Diffusion-limitation of moderate to high surface recombination is not accounted for.
Current-crowding effects are not considered / inaccurate, e.g. recombination at local contacts are overestimated.
Bulk and rear side recombination does not have an impact on quantum efficiency, the collection efficiency is 100%.
There are no “voltage independent carriers”: at voltages below MPP the carrier density saturates of a solar cell at a minimal non-zero level, which is e.g. observable as a non-zero luminescence signal at short-circuit. This is well reproduced with the qn-bulk solver. For the equivalent-circuit mode however, the bulk’s Fermi level split strictly equals the voltage between the skins, meaning there are zero excess charge carriers and luminescence signal at short circuit.

Tandem modeling

In Quokka3 a tandem cell is modelled by including a lumped subcell layer within the skin on the illuminated side, where the subcell is essentially defined via it’s steady-state JV-curve and EQE. The JV-curve is best parameterized into a two-diode model, which extends the equivalent-circtuit skin boundary condition. A detailed description can be found in (Andreas Fell, Schultz-Wittmann, et al. 2022).

*extension of the equivalent-circuit boundary condition to the 3D bulk solver by a two-diode-model representing the top cell of a tandem cell*

The current source, i.e. the top cell \(J_{SC}\) is either a direct input, or determined from the spectrally resolved optical top cell properties. These are essentially the active absorption \(A_{gen}\) (which equals the EQE), and the parasitic absorption \(A_{gen}\) and reflection \(A_{gen}\) at the first light. Those quantities determine the transmission to the bottom cell \(T_{ext}\) and is thus fully compatible with the TextZ model.

*spectrally resolved optical properties of a top cell*

This equivalent-circuit top cell tandem approach as part of the skin concept is well suited and of particular usefulness for the case of a wafer-based c-Si bottom cell. Having the top cell included in the boundary condition to the 3D bulk solver i.e. without any additional numerical mesh, the computational demand is comparable to single junction cell simulations, meaning that the same large 3D geometries can also be used with tandem simulations.

Some limitations resulting from this concept are:

support for steady-state only, i.e. it is not possible to account for e.g. hysterisis effects of perovksite top cells
restricted to two-terminal tandems.

Notably, as an extension to the implementation presented in originally in (Andreas Fell, Schultz-Wittmann, et al. 2022), Quokka3 can account for lateral conductance in the layers between the bottom cell and top cell, and also their contact resistivity which is called junction resistivity.

*updated Fig. 6 of the original tandem JPV paper, showing that Quokka3 now correctly accounts for lateral conduction effects between the bottom and top cell*

Having the top cell JV-curve and EQE as inputs means that the top cell electro-optical characteristics need to be known a-priori when simulating a tandem cell in Quokka3. This shifts the scope of tandem modeling in Quokka3 away from material properties and layer thickness investigations to basic tandem physics and 3D tandem effects (e.g. distributed current matching effects). The lumped inputs also makes a tandem simulation much easier to setup compared to a full drift-diffusion tandem model requiring many material and interface properties as inputs. Some example modeling tasks for which this equivalent-circuit top cell modeling approach can be useful are:

understanding of basic tandem physics with relatively simple input parameters compared to full drift-diffusion simulations
accurate metal grid desing optimization
bottom cell optimization
investigating edge and perimeter effects
investigating impact of lateral nonuniformities

Notably, the tandem model is compatible with both the drift-diffusion and equivalent-circuit bulk model type. When combined with the equivalent-circuit type a 3D distributed equivalent-ciruit model is created, which 1-2 orders of magnitude higher simulation speed and also greatly increased numerical robustness compared to the drift-diffusion type. Also due to some convergence challanges for the tandem approach with the drift-diffusion type, the equivalent-circuit bulk model tpye is thus the recommended option for tandem simulations.

Solving detailed physics: 1D detailed solver

The 1D detailed solver can account for the full general models, however solving them in 1 dimension only. Meaning that the currents and the electric field have only a z-component, and \(\nabla\) is reduced to the derivation in z-direction. The 1D detailed solver can solve two types of solution domains:

A 1D semiconductor device, i.e. 1D solar cell. This closely resembles the functionality of the popular solar cell simulator PC1D (Clugston and Basore 2010). It requires contacts on the front and rear.
A semiconductor skin, i.e. a detailed skin only. Here on the bulk-side of the skin a quasi-neutral boundary condition is applied. It is called the skin solver and plays a decisive role in multiscale modeling.

Boundary conditions

The boundary conditions for the detailed solver need to define the electron and hole current densities towards the boundary \(J_{el,b}\) and \(J_{hol,b}\), as well as the electric field \(E_b\). Not that this is NOT a lumped boundary or skin, but represents the actual surface or interface to metal.

All boundary types can be assigned surface recombination. The resulting recombination current density \(J_{rec}\) is treated independently (i.e. additionally) to the current transport model over the boundary. This is different to some other semiconductor solvers like for example Sentaurus Device, where recombination is the actual current transport model. The physical meaning of an apparently identical setting can therefore be quite different, and consequently give different results.

If not set directly, the polarity and the potential applied to the contact is defined by the covering metal feature. When solving a cell in 1D, for ‘n-type’ polarity the applied potential is \(\varphi_{met}=V_{int}/2\), for ‘p-type’ it is \(\varphi_{met}=-V_{int}/2\) (different definition to the qn-bulk solver). Unless otherwise enforced, Quokka3 assumes that current transport over the boundary is dominated by one type of carrier, i.e. assumes a pure electron contact or hole contact determined by the polarity. This means current transport is calculated only for the majority carrier type, and set to zero for the minorities. Note that this doesn’t necessarily mean perfect selectivity, as recombination still leads to a non-zero minority current density towards the boundary: \[ J_{min}=J_{hol,b}=J_{rec}\ (electron\ contact) \]

\[ J_{min}=J_{el,b}=-J_{rec}\ (hole\ contact) \]

Surface SRH

Simplified steady-state surface SRH can be defined at all boundary types. The expression is equivalent to the simplified bulk SRH expression, using the fundamental surface recombination velocities for electrons and holes \(S_{n0}\), \(S_{p0}\): \[ J_{rec}=\frac{p\,n-n_{i,eff}^2}{S_{n0}^{-1}(p_1+p)+S_{p0}^{-1}(n_1+n)} \] The user can directly input \(S_{n0}\), \(S_{p0}\), or choose a published parameterization [MATERIAL].

Non-contacted boundary

For a non-contacted boundary no net current can be extracted, i.e. \(J_{el,b}+J_{hol,b}=0\), meaning that only the recombination current density s present: \[ J_{el,b}=-J_{hol,b}=-J_{rec} \] The electric field at the boundary is determined by the surface charge \(Q_s\): \[ E_{b}=\frac{qQ_s}{\varepsilon_r\varepsilon_0} \]

Ideal contact

An ideal contact assumes that the majority carrier type can freely flow over the boundary. It also assumes flat-band conditions, meaning zero electric field. Note however that surface recombination is independently defined, and can thus lead to some bending of the minority Fermi level towards the surface. \[ \phi_{maj}=\varphi_{met} \]

\[ E_b=0 \]

TLS / MIS contact

A more detailed description of Quokka3’s approach to model passivating contacts via metal-insulator-semiconductor (MIS) theory is given in (Fell et al. 2018).

TLS / MIS means a “transport-layer - semiconductor” or “metal - insulator - semiconductor” contact model, respectively. The MIS model can also be applied to materials other than metals, in which case those materials are assumed to be “metal-like”. That means the material is assumed to have a unique potential throughout its thickness, meaning that band-bending effects within the material layer are neglected. The MIS model may thus be applicable e.g. to semiconductor material layers with very high doping densities, like highly doped poly-Si. A “TLS” boundary is more general a contact between a transport-layer material and the semiconductor, and does account for band-bending effects in the transport layer, at the cost of not supporting direct tunneling through a thin insulator. The TLS boundary is enabled when selecting the “analytic band-bending” type of a transport layer.

Quokka3’s MIS boundary condition is scoped to be a compromise between the complex fully detailed approach and the potentially oversimplifying lumped skin approach. It is able to describe some effects specifically observed for passivating contacts without adding too much complexity, i.e. input parameters. Those effects are for example non-ohmic contact resistivity of tunnel-oxide based contacts, non-ideal recombination and \(iV_{oc}\) - \(V_{oc}\) losses.

*MIS model as a compromise between a fully detailed and lumped skin approach to model passivating contacts*

Three cases of the MIS boundary are distinguished in Quokka3: MS, MIS and MIS_simple

MS case (Schottky contact)

For the MS case, the current transport of electrons and holes over the boundary is governed by thermionic emission (TE): \[ J_{TE,el/hol}=A_0T^2 \exp\left(-\frac{\Phi_{B,el/hol}+\Delta\phi_I}{V_t}\right) \left(\exp\left(\frac{\varphi_{met}-\phi_{el/hol}}{V_t}\right)-1\right) \] where \(\Phi_{B,el}+\Phi_{B,hol}=E_g\) are the Schottky barriers for electrons and holes, which sum equals the band-gap. By using the Schottky-Mott rule, the Schottky barrier can be calculated from the difference of the metal workfunction and the absorber’s electron affinity, however note that in reality the actual Schottky barrier is always substantially lower due to various barrier-lowering effects (see(Tung T. Raymond 2014)). The potential drop over the insulator \(\phi_I\) is zero for the MS case.

The current boundary condition is achieved by adding surface recombination: \[ J_{MS,el/hol}=J_{TE,el/hol}\pm J_{rec} \] The boundary condition for the electrostatic potential is given by defining the conduction band edge: \[ E_c=\varphi_{met}+\phi_{B,el} \] Note that Quokka3’s independent treatment of current transport over the boundary and surface SRH is different to other simulation tools, e.g., Sentaurus Device. There no surface SRH is allowed for an MS contact, and the \(S_{n0}\) and \(S_{p0}\) inputs are the “thermionic emission velocities.” The TE current when using an unscaled Richardson constant implies very high (∼\(10^7\) cm/s) thermionic emission velocities, which effectively results in high surface recombination, meaning that surface SRH should not additionally be applied. However, Quokka3 allows to disable the thermionic emission current for the minority carriers, rendering the effective recombination caused by thermionic emission to zero, which makes it then sensible to define surface recombination. While this is a less physical MS model, it is well suited to effectively describe a passivating contact by already resolving the interplay of band bending and surface recombination. A similar approach may fail when “misusing” the thermionic emission velocities for the same purpose, as those are not SRH input parameters. Consider the extreme example: setting \(S_{n0}\) and \(S_{p0}\) to very low values results in a highly non-ohmic thermionic emission contact resistivity for majority carriers, and consequently, S-shaped J–V curves, but may be an excellent contact in Quokka3 when disabling minority thermionic emission.

The implementation of the MS boundary condition in Quokka3 further allows to consider the influence of a semiconductor buffer layer in a very simplistic way. Assuming that the buffer layer does not cause any transport limitations, the thermionic emission barriers are increased by the respective band offsets, \(\phi_{B,el,offset}\) and \(\phi_{B,hol,offset}\), to account for the shift of the TE barrier interface; see Fig. 2. Those offsets are added to the respective thermionic emission barrier, but not to the electrostatic potential boundary, meaning that the band bending in the absorber is unaffected.

MIS case (direct tunneling contact)

*MIS boundary condition; band diagram (left) and equivalent circuit representation (right)*

Additional to the MS case, the MIS case considers a thin insulator layer between the metal and the absorber via a simple direct tunneling model. A rectangular tunneling barrier is assumed characterized by the barrier thickness \(d_{tun}\) and barrier height \(\Phi_{tun}\). The resulting tunneling probability \(P_{tun}\) scales the TE current: \[ J_{MIS,el/hol}=P_{tun,el/hol}J_{TE,el/hol}\pm J_{rec}\left(+J_{ohmic}\right) \]

\[ P_{tun,el/hol}=\exp\left(-\frac{2d_{tun}\sqrt{2\,q\,m_{tun,el/hol}^*\left(\Phi_{tun,el/hol}-\Delta\phi_I/2\right)}}{\hbar}\right) \] Compared to the MS case, the band bending in the absorber is reduced by the potential drop over the insulator layer \(\Delta\phi_I\) , expressed via the electric field at the absorber surface \(E_{MIS}\) being a function of the permittivity ratio: \[ \vec{E}_{MIS}=\frac{\Delta\phi_I}{d_{tun}}=\frac{\varepsilon_{r,tun}}{\varepsilon_{r,skin}}\times\frac{\varphi_{met}+\Phi_{B,el}-E_{c}}{d_{tun}} \]

The insulator potential drop \(\Delta\phi_I\) also impacts the thermionic emission current by increasing the thermionic emission barrier.

Furthermore, one can define an ohmic contact resistivity \(\rho_{c,ohmic}\) in parallel to the majority carrier thermionic emission / tunneling current path, adding \(J_{ohmic} = \frac{\varphi_{met}-\phi_{el/hol}}{\rho_{c,ohmic}}\) to the majority carrier transport boundary condition. It can account for lateral inhomogeneities and / or pinholes resulting in an ohmic contact fraction. Note that this is not (yet) sufficient to properly model arbitrarily mixed tunneling–pinhole contacts as the influence of the ohmic contact fraction on the recombination is not considered.

TLS case

In case the material in contact is not very highly doped, and band-bending effects within that layer should be considered, Quokka3 offers the “TLS” boundary condition. This “analytic band-bending” type of a transport layer is in particular useful to model electron- and hole transport layers (ETL and HTL) in perovskite solar cells, but may well be useful for many passivating contact materials used for c-Si solar cells. The downside to the MIS model is that the analytical band-bending model does not yet support to account for tunneling over an insulator.

When activating the TLS boundary condition via choosing the “analytic band-bending” model type of a transport layer, some material properties are required as inputs. These are the band-gap \(E_g\), electron affinity \(\chi\), relative permittivity \(\epsilon_r\), activation energy \(E_{act}\) and doping density \(N_{dop}\).

*TLS boundary condition, exemplary for a electron transport layer; the band-bending in the TL is analytically modelled and reduces the band-bending in the absorber compared to the MS case*

The current boundary condition is the same as for the MS case above. The electric field boundary condition is given by (???): \[ \vec{E}_{TLS}=sign(\Phi_{bb,TL})\frac{\sqrt{2}V_t}{L_D}\sqrt{\frac{u_b+\frac{\Phi_{bb,TL}}{V_t}}{\cosh(u_b)}-\frac{\Phi_{bb,TL}}{V_t}\tanh(u_b)-1} \] and \[ \Phi_{bb,TL}=E_{c,A}-\Phi_B-\varphi_{appl}-E_{act} \] where \(L_D=\sqrt{\frac{\epsilon_r\epsilon_0 V_t}{qN_{dop}}}\) is the Debye length, \(u_b=\frac{E_{act}-E_g/2}{V_t}\) and \(E_{c,A}\) is the conduction band edge on the absorber side of the boundary. The model for a hole transport layer is equivalent.

The analytical model is accurate if the thickness of the layer is much larger than the Debye length, and if the majority quasi Fermi level is flat throughout the thickness of the transport layer, i.e. if vertical transport losses through the layer’s bulk are negligible. These are often good assumptions, but may fail e.g. for heterojunction structures where generation and recombination as well as transport effects within the layers are significant.

Quasi-neutral boundary condition

The quasi-neutral boundary condition is applied on the quasi-neutral bulk-side of a skin when using the skin solver to solve a detailed skin only. A given Fermi-level split and the electrical potential calculated numerically by assuming quasi-neutrality is applied: \[ \phi_{el}-\phi_{hol}=\phi_{split} \]

\[ \varphi_{e,b}=\varphi_{e,qn}(\phi_{split}) \]

Parameterizing detailed skins: skin solver

The skin solver is the 1D detailed solver applied to a detailed skin only, using a quasi-neutral boundary condition on the bulk-side. The main purpose is to determine the lumped parameters representing the detailed skin, e.g. determine the \(J_{0,skin}=J_{0e}\) and \(R_{sheet}\) of a phosphorus emitter diffusion.

To solve a skin domain, an operating point needs be defined by the bulk-side Fermi-level split and the net current density through the skin. From the solution, the lumped parameters of the general parameterization which exactly represent the skin’s operating point, are extracted as follows:

\(1/R_{sheet}\) is calculated for each carrier type by integrating the conductivity. The lower \(R_{sheet}\) determines the majority carrier type, i.e. whether the skin is ‘n-type’ or ‘p-type’. Quokka3 issues a warning if both types have low \(R_{sheet}\), which may lead to inaccuracies within multiscale modeling and should not happen in typical silicon solar cell skins.
\(J_{G,skin}\) is found by integrating the generation rate (if the skin is illuminated)
The skin potential \(\varphi_{skin}\) is assigned the majority potential of the mesh element with the highest majority carrier conductivity (usually the mesh element next to the surface).
The vertical resistance \(\rho_{skin}\) is calculated by the difference of the skin potential and the bulk-side majority carrier quasi-Fermi potential over the net current density.
The total recombination current density \(J_{rec,light}\) is determined by integrating volume recombination and adding interface recombination.
If the skin is illuminated a second simulation in the dark simulation is performed.
The total recombination current density \(J_{rec,dark}\) is determined by integrating volume recombination and adding interface recombination (which (only) in the dark equals the minority carrier current density through the bulk-side boundary).
Any additional recombination in the illuminated case compared to the dark case is attributed to carriers generated within the skin, so that the collection efficiency is defined as: \[ \eta_c=1-\frac{J_{rec,light}-J_{rec,dark} }{J_{G,skin}} \]
Consistent with above, \(J_{0,skin}\) for both the illuminated and dark case is defined as: \[ J_{0,skin}= \frac{J_{rec,dark}}{\exp\left(\phi_{split}⁄V_t\right)-1} \] Note that different to some other available software tools providing a similar functionality, Quokka3 does not assume quasi-neutrality within the skin, and thus also includes the space-charge-region (SCR) if present. While this might not make much difference for e.g. typical highly doped near-surfaces like emitter diffusions, the results are more accurate in general.

Multiscale solver

*multiscale concept: detailed skin + 3D lumped skin solver is equivalent to a full detailed 3D simulation*

The multiscale solver couples the “small scales” of the skins using the skin solver with the qn-bulk solver solving the “large scale” 3D geometry. It used whenever a skin is defined via its detailed instead of its lumped properties. Then the skin solver is used to parameterize the respective skin and use the lumped parameters as a boundary condition for the qn-bulk solver subsequently simulating the 3D geometry.

Coupling degree

Quokka3 supports 3 different degrees of coupling a detailed skin to the qn-bulk which enables the user to balance speed vs. accuracy:

single point: simulates and parameterizes the skin once at a single operating point, and assumes the parameters to be constant. This is fastest and most robust, and is well valid e.g. for highly doped near-surfaces like e.g. diffused emitters or Al-BSF, which usually have a constant \(J_{0,skin}\) and a negligible vertical resistivity.
injection dependent recombination: simulates and parameterizes the skin for the required range of bulk-side Fermi-level splits and a at constant (and low) net current density. This is useful e.g. for lowly-doped surfaces passivated by a layer with moderate charge, which can show significant injection dependence of the recombination, but where vertical resistance is not relevant. This option should also be considered for typical doped near-surfaces, like e.g. light emitters, to increase accuracy.
full coupling: simulates and parameterizes the skin for all relevant operating points. This is the generally valid and thus the most accurate option, but at the expense of simulation time. The general coupling is required e.g. to correctly account for resistive losses in the space-charge-region (SCR).

Texture multiplier

In Quokka3 the surfaces of the model domain are always strictly planar, a simplification used also for the vast majority of silicon solar cell device simulations with other tools. A complication arises when the actual surface is textured. This generally increases the effective recombination of such a textured skin due to the surface area enlargement. When assuming perfectly conformal skin properties along the textured surface (which is not necessarily a good assumption, but not to be discussed here), one could expect the recombination to scale with the geometric surface area enlargement (e.g. ~1.7 for random pyramid texture). However, due to current transport limitations within the pyramidal volume, a lower increase of e.g. \(J_{0,skin}\) is typically measured, see e.g. (Baker-Finch et al. 2012).

For lumped skin modeling this states no problem, as the required lumped properties can be measured on textured surfaces. For detailed skin modeling (and subsequent multiscale modeling) using a planar surface however, the recombination needs to be increased to represent the textured case. This is usually done by scaling surface recombination (and sometimes Auger) by such a factor, that the experimentally measured increase of \(J_{0,skin}\) (textured vs. planar) is replicated. While this is suitable to result in the correct effective recombination (i.e. matching \(V_{oc}\)), it leads to an unphysical distribution of the recombination. This is in particular noticeable in the skin’s collection efficiency (i.e. the IQE): due to the conformal assumption the a correct modelling of collection efficiency must use the actual surface recombination, and increasing it leads to an unphysical decrease of skin’s collection efficiency. Therefore a single fully detailed cell simulation can not simultaneously predict \(J_{0,skin}\) (\(V_{oc}\)) and collection efficiency (\(J_{sc}\)).

The multiscale model of Quokka3 however can resolve this problem: as both \(J_{0,skin}\) and collection efficiency are available as separate intermediate results, the texture multiplier is applied to \(J_{0,skin}\) skin only. Note this definition of the texture multiplier is different to the common definition (which increases physical recombination parameters). Next to the advantage described above, Quokka3’s definition can be considered the practically more meaningful one, as it is the directly measurable increase of effective surface recombination, i.e. \(\frac{J_{0,skin,textured}}{J_{0,skin,planar}}\).

Loss analysis

Background

One of the major purposes of a solar cell simulation is to establish a loss analysis. A loss analysis provides a quantitative breakdown of avoidable, and also unavoidable, loss mechanisms, to identify bottlenecks and improvement potentials of the investigated cell design and properties.

There are several established methods how to perform a loss analysis, all with it’s own benefits and limits. The easiest one typically performed in numerical simulations is the current loss analysis (CLA), which lists optical and recombination losses in units of \(mA/cm^2\). The CLA is straightforward to interpret, but has the disadvantage of excluding transport losses. An approach to include both recombination and transport losses in units of power \(mW/cm^2\) is the free energy loss analysis (FELA) as presented in (Brendel et al. 2008). An extension of the FELA to accurately include optical losses has been presented in (Greulich et al. 2013). While being able to compare all losses in same power units, and to identify the main loss mechanisms, the detailed interpretation of the FELA is often difficult due to non-intuitive meaning of free energy, and the associated transport losses (see FELA section below). Finally, the synergistic efficiency gain analysis (SEGA) has been suggested in (Brendel et al. 2016). In the SEGA the various loss mechanisms are subsequently disabled in several simulation runs to quantify the associated efficiency gain, and can accurately quantify the path towards an ideal Auger-limited device, and thereby overcomes the shortcomings of the FELA. However, synergistic effects come in to play, which complicate the interpretation due to the results being dependent on the chosen order of disabling the individual losses.

In Quokka3, the CLA and a PLA (power loss analysis) are implemented, where the PLA is the FELA with an approximate inclusion of optical losses.

Unique operating point

The losses depend on the operating point of the solar cell, i.e. the illumination intensity and terminal current / voltage. Therefore, naturally, a loss analysis can only be performed if a defined unique operating point exists. This means that the chosen solution type must either be a unique operating point (e.g. the ‘single JV-point’ solution type), or must otherwise define a unique operating point as part of the simulation. This for now is the case for the automatic ‘light JV-curve’ solution type (i.e. if the Vterm steps are not user defined): here the unique operating point is the maximum power point (MPP) which means that the loss analysis given in the results is valid at MPP only.

CLA: current loss analysis

The CLA in Quokka3 lists all optical and recombination losses in units of current density \(mA/cm^2\) which are available depending on how the solar cell properties have been defined. For example, if a generation profile file is used to directly define the generation, no optical modeling is performed and therefore no optical losses are listed except for shading by the metallization. The most detailed optical loss breakdown in Quokka3 is achieved when using the extended Basore model, which then comprises external and escape reflection, transmission and various parasitic absorption contributions.

Balances

The sum of the optical losses determine the generated current density \(J_{gen}\) from the usable incident illumination \(J_{illum}\):

\[ J_{gen}=J_{illum}-\sum J_{loss,optical} \]

The usable illumination is determined by integrating the incident spectrum up to the absorption cut-off (~band gap) of the main absorber, i.e. bulk, material. For the default Si material the cut-off is set at 1200 nm. For the default AM1.5g spectrum this results in \(46.43 mA/cm^2\) of usable illumination.

Care needs to be taken when calculating the balance equation above, as some loss contributions are redundant and must be included only once:

The total absorption \(A\) is the sum of active band-to-band absorption in the bulk \(A_{bb,bulk}\) and all parasitic absorption losses
the external transmission loss \(T_{ext}\) is the sum of external reflection \(R_{ext}\) and external absorption \(A_{ext}\) losses
the light-trapping losses are the sum of escape reflection \(R_{esc}\), transmission \(T\) and all parasitic absorption losses

The electrical losses and the external circuit shunt loss (if the external circuit is enabled) \(J_{shunt}\) determine the difference between the generated current density \(J_{gen}\) and terminal current density \(J_{term}\):

\[ J_{term}=-(J_{gen}-J_{shunt}-\sum J_{loss,recombination}) \]

avoidable and non-avoidable light-trapping losses

Quokka3 differentiates light-trapping into the avoidable and non-avoidable part. It does so by effectively performing a second internal optical simulation in which all light-trapping properties are set to their theoretical Lambertian limit, i.e. disregarding all internal parasitic absorption and setting ideal internal Lambertian reflection at the surfaces. This can be considered the practically achievable limit. The unavoidable part are then the remaining light-trapping losses of this ideal device, which physically can not be zero but are in fact substantial. The purpose is to more easily interpret the large light-trapping losses, by directly seeing which part of it is avoidable / practically improvable.

Notably, the differentiation is valid for the given cell thickness only. An improvement of light-trapping by increasing the cell thickness is also possible, but not considered in the avoidable part.

PLA: power loss analysis

FELA

The power loss analysis (PLA) is based on the free energy loss analysis (FELA), see (Brendel et al. 2008), which calculates recombination and transport losses in units of power density \(mW/cm^2\).

For generation and recombination this is done by locally multiplying the generation or recombination rates with the quasi Fermi level split and subsequently integrate over the solution domain. This leads to the non-intutive issue that the contributions are dependent on the Fermi level split, i.e. dependent on the device voltage at the operating point. This explains why in the FELA the terminal power and losses sum up to the generated power only, which is the generated free energy and significantly lower than e.g. the Auger-limit. The generated power should thus not be considered an achievable upper limit of the device, as the generated power would itself increase when improving the device and thereby also the terminal voltage. However, the relative comparison of individual loss contributions is still valid and meaningful to identify bottlenecks of the investigated cell design and properties.

Another complexity of the FELA regards the transport losses in the bulk, see also explanation in (Brendel et al. 2008). The FELA transport losses is total free energy transport loss, which is different from the normal understanding of transport losses, namely the Joule’s heat / resistive losses. The total transport loss is the sum of the resistive and diffusive loss, which is equivalent to the differentiation in current contributions when speaking oft drift-diffusion transport. The resistive part can be interpreted as the one observable as a ‘series resistance’ in a diode-model representation of the solar cell’s JV-curve, and is typically dominated by the majority carrier transport. The diffusive part is less obvious to interpret. It is commonly related to diffusion of carriers from the generation site to a recombintion site, and typically significant for minority carriers, or both carriers in high-injection conditions. E.g. for rear-junction devices and in particual IBC devices diffusive losses can be significant by minority carriers having to diffuse from the front side to the rear emitter.

A notable wrong conclusion can be that wafer doping should be increased when high transport losses are observed. If for a particular device the transport losses are dominated by diffusive transport losses, the opposite may be true in order to increase minority carrier density and thus reduce the diffusive transport loss. It is therefore advisable to interpret high bulk transport losses carefully for correct conclusions.

Optical power losses

The original FELA does not include optical losses. In order to include all electro-optical losses in a single comparable list, Quokka3 includes them by an approximate approach. Namely, the optical loss contributions from the CLA are all multiplied by the internal voltage, i.e. the external series resistance free voltage, to achieve the same units of \(mW/cm^2\). This effectively assumes that the quasi Fermi level split is constant throughout the solution domain and equals the internal voltage, which is a reasonable approximation for well-working solar cells. A more accurate approach would be to also consider the local quasi Fermi level split, as in (Greulich et al. 2013), which is however not implemented in Quokka3.

Due to the reasons explained in the FELA section above, and also because of the mentioned approximation, the sum of generated power and optical power loss do NOT result in the incident power (\(100 mW/cm^2\) for AM1.5g).

The optical power loss balances are the same as for the CLA.

PLA equals efficiency loss

Coincidentally, when using the AM1.5g spectrum with \(100 mW/cm^2\) illumination intensity, the PLA numbers can be interpreted also in units [%_{abs}] efficiency loss, making it most convenient to interpret.

Material properties

TODO

List of symbols

Note that the units are chosen to be consistent, but are not necessarily equal to the units of the corresponding input parameters of the settingsfile.

symbol	unit	description
\(A_0\)	\(Acm^{-2}K^{-2}\)	Richardson constant
\(A_G\)		actively absorbed fraction of photons after being transmitted through the illuminated surface
\(A_{ppp}\)		absorbed fraction per perpendicular pass
\(c_{an}\), \(c_{cat}\)	\(cm^{-3}\)	anion, cation density
\(\vec{E}\)	\(V/cm\)	electric field
\(E_b\)	\(V/cm\)	electric field at boundary
\(E_g\)	\(eV\)	band gap
\(E_c\), \(E_v\)	\(eV\)	conduction band edge and valence band edge
\(E_{act}\)		activation energy
\(f_t\)		occupied trap density fraction
\(G\)	\(cm^{-3}s^{-1}\)	generation rate of electron-hole pairs
\(\vec{J}_{el}\), \(\vec{J}_{hol}\), \(\vec{J}_{an}\), \(\vec{J}_{cat}\)	\(A/cm^2\)	charge carrier current density
\(J_{cont}\), \(J_{pad}\)	\(A/cm^2\)	current density through contact and between metal and pad layer
\(k_B\)	\(VAs/K\)	Boltzmann constant
\(l_1\)	\(cm\)	optical pathlength of first path through the device
\(m_{tun}^*\)	\(kg\)	electron or hole tunneling mass
\(N_{0,an}\), \(N_{0,cat}\)	\(cm^{-3}\)	average / initial homogenous anion and cation concentration
\(N_A\), \(N_D\)	\(cm^{-3}\)	acceptor and donor density
\(N_c\), \(N_v\)	\(cm^{-3}\)	density of states in conduction band, valence band
\(N_t\)	\(cm^{-3}\)	total trap or defect density
\(n\)	\(cm^{-3}\)	electron density
\(\vec{n}\)		normal vector on boundary
\(n_i\)	\(cm^{-3}\)	intrinsic carrier density
\(n_t\)	\(cm^{-3}\)	occupied trap density
\(p\)	\(cm^{-3}\)	hole density
\(Q_s\)	\(cm^{-2}\)	surface charge
\(q\)	\(As\)	elementary charge
\(R\)	\(cm^{-3}s^{-1}\)	recombination rate of electron-hole pairs
\(R_{front,avg}\), \(R_{rear,avg}\)		average (meaning n-th pass) internal reflectance
\(R_{front,spec}\)		specular internal reflectance at front side
\(R_{front,1}\)		first pass internal reflectance at front side
\(R_{sheet}\)	\(\Omega\)	sheet resistance
\(S_{n0}\), \(S_{p0}\)	\(cm/s\)	fundamental electron and hole surface recombination velocity
\(T\)	\(K\)	temperature
\(T_1\), \(T_2\), \(T_n\)		(effective) transmission of first, second and n-th pass (internal optics)
\(T_{ext}\)		transmitted fraction of external light through the surface
\(V_{int}\)	\(V\)	internal voltage = ‘n-type’ polarity applied voltage
\(V_t\)	\(V\)	thermal voltage, equals \(k_BT/q\)
\(W\)	\(cm\)	device thickness
\(Z\)		(spectrally dependent) pathlength enhancement factor
\(Z_0\), \(Z_\infty\), \(Z_p\)		\(Z\) parameters
\(\alpha\)	\(cm^{-1}\)	absorption coefficient
\(v_{th}\)	\(cm/s\)	thermal velocity
\(\delta E_{BGN,c}\), \(\delta E_{BGN,v}\)	eV	band gap narrowing of conduction band and valence band
\(\varepsilon_0\)	\(AsV^{-1}cm^{-1}\)	vacuum permittivity
\(\varepsilon_r\)		relative permittivity
\(\zeta\)	\(cm\)	closest distance to surface coordinate
\(\eta_c\)		collection efficiency
\(\theta_{surf}\), \(\theta_{bulk}\)		angle of optical path relative to (textured) surface, and (planar) bulk, respectively
\(\theta_1\) (\(=\theta_{bulk}\)), \(\theta_2\), \(\theta_n\)		first, second and n-th pass path-angle (internal optics)
\(\kappa_{el}\), \(\kappa_{hol}\), \(\kappa_{an}\), \(\kappa_{cat}\)	\(AV^{-1}cm^{-1}\)	charge carrier conductivity
\(\Lambda_{rear}\)		Lambertian fraction of rear side (internal optics)
\(\lambda\)	\(nm\)	wavelength
\(\mu_{el}\), \(\mu_{hol}\), \(\mu_{an}\), \(\mu_{cat}\)	\(cm^2V^{-1}s^{-1}\)	charge carrier mobility
\(\rho_{cont}\), \(\rho_{pad}\)	\(\Omega cm^2\)	skin - metal and metal - pad contact resistivity
\(\sigma_n\), \(\sigma_p\)	\(cm^{-2}\)	electron and hole capture cross section
\(\tau_{n0}\), \(\tau_{p0}\)	s	fundamental SRH electron and hole lifetime
\(\Phi_B\)	\(eV\)	Schottky barrier height
\(\phi_{el}\), \(\phi_{hol}\), \(\phi_{an}\), \(\phi_{cat}\)	\(V\)	charge carrier electrochemical potential (=quasi-Fermi potential)
\(\phi_{maj}\)	\(V\)	quasi-Fermi potential of majority carriers
\(\phi_{split}\)	\(V\)	quasi-Fermi potential split
\(\varphi_e\), \(\varphi_{skin}\), \(\varphi_{met}\), \(\varphi_{pad}\)	\(V\)	electrostatic, skin, metal and pad potential
\(\Delta\varphi_{maj}\)	\(V\)	difference between skin potential and bulk-side majority quasi-Fermi potential

References

Altermatt, Pietro P. 2011. “Models for Numerical Device Simulations of Crystalline Silicon Solar Cells—a Review.” Journal of Computational Electronics 10 (3): 314–30. https://doi.org/10.1007/s10825-011-0367-6.

Baker-Finch, Simeon C., Keith R. McIntosh, Mason L. Terry, and Yimao Wan. 2012. “Isotextured Silicon Solar Cell Analysis and Modeling 2: Recombination and Device Modeling.” IEEE Journal of Photovoltaics 2 (4): 465–72.

Basore, Paul A. 1993. “Extended Spectral Analysis of Internal Quantum Efficiency.” In 23rd Photovoltaic Specialists Conference, 147–52. IEEE; IEEE.

Brendel, R., S. Dreissigacker, N.-P. Harder, and P. P. Altermatt. 2008. “Theory of Analyzing Free Energy Losses in Solar Cells.” Applied Physics Letters 93 (17): 173503. https://doi.org/10.1063/1.3006053.

Brendel, Rolf. 2012. “Modeling Solar Cells with the Dopant-Diffused Layers Treated as Conductive Boundaries.” Progress in Photovoltaics: Research and Applications 20 (1): 31–43. https://doi.org/10.1002/pip.954.

Brendel, Rolf, Thorsten Dullweber, Robby Peibst, Christopher Kranz, Agnes Merkle, and Daniel Walter. 2016. “Breakdown of the Efficiency Gap to 29% Based on Experimental Input Data and Modeling.” Progress in Photovoltaics: Research and Applications 24 (12): 1475–86. https://doi.org/10.1002/pip.2696.

Clugston, D. A., and P. A. Basore. 2010. “PC1D Version 5: 32-Bit Solar Cell Modeling on Personal Computers.” In 35th Photovoltaic Specialists Conference, 207–10. IEEE. https://doi.org/10.1109/pvsc.1997.654065.

Fell, A. 2013. “A Free and Fast Three-Dimensional/Two-Dimensional Solar Cell Simulator Featuring Conductive Boundary and Quasi-Neutrality Approximations.” IEEE Transactions on Electron Devices 60 (2): 733–38. https://doi.org/10.1109/ted.2012.2231415.

Fell, A., F. Feldmann, C. Messmer, M. Bivour, M. C. Schubert, and S. W. Glunz. 2018. “Adaption of Basic Metal–Insulator–Semiconductor (Mis) Theory for Passivating Contacts Within Numerical Solar Cell Modeling.” IEEE Journal of Photovoltaics. https://doi.org/10.1109/JPHOTOV.2018.2871953.

Fell, Andreas, Johannes Greulich, Frank Feldmann, Christoph Messmer, Jonas Schoen, Martin Bivour, and Martin C. Schubert, and Stefan W. Glunz. 2022. “Modelling Parasitic Absorption in Silicon Solar Cells with a Near-Surface Absorption Parameter.” Solar Energy Materials and Solar Cells 236: 111534. https://doi.org/10.1016/j.solmat.2021.111534.

Fell, Andreas, Keith R. McIntosh, and Kean C. Fong. 2016. “Simplified Device Simulation of Silicon Solar Cells Using a Lumped Parameter Optical Model.” IEEE Journal of Photovoltaics, 1–6. https://doi.org/10.1109/JPHOTOV.2016.2528407.

Fell, Andreas, Tim Niewelt, Bernd Steinhauser, Friedemann D. Heinz, Martin C. Schubert, and Stefan W. Glunz. 2021. “Radiative Recombination in Silicon Photovoltaics: Modeling the Influence of Charge Carrier Densities and Photon Recycling.” Solar Energy Materials and Solar Cells 230: 111198. https://doi.org/10.1016/j.solmat.2021.111198.

Fell, Andreas, Jonas Schoen, Martin C. Schubert, and Stefan W. Glunz. 2017. “The Concept of Skins for Silicon Solar Cell Modeling.” Solar Energy Materials and Solar Cells, 173 (Supplement C): 128–33. https://doi.org/10.1016/j.solmat.2017.05.012.

Fell, Andreas, Oliver Schultz-Wittmann, Christoph Messmer, Martin C. Schubert, and Stefan W. Glunz. 2022. “Combining Drift-Diffusion and Equivalent-Circuit Models for Efficient 3D Tandem Solar Cell Simulations.” IEEE Journal of Photovoltaics 12 (6): 111534. https://doi.org/10.1109/JPHOTOV.2022.3209608.

Fell, Andreas, Wiebke Wirtz, Hannes Höffler, and Johannes Greulich. 2019. “Determining the Generation Rate of Silicon Solar Cells from Reflection and Transmission Measurements by Fitting an Analytical Optical Model.” In Silicon Pv 2018, 3037–41. https://doi.org/10.1109/PVSC40753.2019.8980730.

Green, Martin A. 2002. “Lambertian Light Trapping in Textured Solar Cells and Light–Emitting Diodes: Analytical Solutions.” Progress in Photovoltaics: Research and Applications 10 (4): 235–41.

Greulich, Johannes, Hannes Höffler, Uli Würfel, and Stefan Rein. 2013. “Numerical Power Balance and Free Energy Loss Analysis for Solar Cells Including Optical, Thermodynamic, and Electrical Aspects.” Journal of Applied Physics 114 (20): –. https://doi.org/10.1063/1.4832777.

Macdonald, Daniel, and Andrés Cuevas. 2003. “Validity of Simplified Shockley-Read-Hall Statistics for Modeling Carrier Lifetimes in Crystalline Silicon.” Physical Review B 67 (7): 075203.

McIntosh, K. R., and S. C. Baker-Finch. 2015. “A Parameterization of Light Trapping in Wafer-Based Solar Cells.” IEEE Journal of Photovoltaics PP (99): 1–8. https://doi.org/10.1109/JPHOTOV.2015.2465175.

Nguyen, Hieu T., Simeon C. Baker-Finch, and Daniel Macdonald. 2014. “Temperature Dependence of the Radiative Recombination Coefficient in Crystalline Silicon from Spectral Photoluminescence.” Applied Physics Letters 104 (11): 112105. https://doi.org/10.1063/1.4869295.

Rand, James A., and Paul A. Basore. 1991. “Light-Trapping Silicon Solar Cells-Experimental Results and Analysis.” In 22nd Photovoltaic Specialists Conference, 192–97. IEEE; IEEE.

Sabater, Aina Alapont, Andreas Fell, Andreas A. Brand, Matthias Müller, and Johannes M. Greulich. 2021. “Parameterization of the Back-Surface Reflection for Perc Solar Cells Including Variation of Back-Contact Coverage.” IEEE Journal of Photovoltaics 11 (5): 1136–40. https://doi.org/10.1109/JPHOTOV.2021.3082402.

Schinke, Carsten, David Hinken, Jan Schmidt, Karsten Bothe, and Rolf Brendel. 2013. “Modeling the Spectral Luminescence Emission of Silicon Solar Cells and Wafers.” IEEE Journal of Photovoltaics 3 (3): 1038–52. https://doi.org/10.1109/JPHOTOV.2013.2263985.

Shockley, We, and W. T. Read Jr. 1952. “Statistics of the Recombinations of Holes and Electrons.” Physical Review 87 (5): 835.

Tung T. Raymond. 2014. “The Physics and Chemistry of the Schottky Barrier Height.” Applied Physics Reviews 1 (1): 011304. https://doi.org/10.1063/1.4858400.

Quokka3 Modelling Guide

Concept: Bulk, Skin, Contact, Metal and Pad

Meshing

Optical modeling

Solving generation: the Text-Z model

Relation to spectral response (SR)

Calculating the generation profile

Different ways to define light-trapping / \(Z\)

Z parameterization

Temperature (in-)dependence

Modeling light-trapping and luminescence: the extended Basore model

Model description

Accounting for skin parasitic absorption via Appp

Luminescence, photon reabsorption and photon recycling

Inhomogeneous properties

Electrical modeling

General semiconductor model equations

Volume carrier transport

Carrier statistics

Electrons and holes

Ions

Generation and recombination

General SRH

Simplified SRH

Solving the bulk with lumped skins: qn-bulk solver

Volume equations

Boundary conditions

Symmetry / perfect edge boundary

Lumped skin and metal layer boundary

Resistive mode

Equivalent-circuit bulk model type

Tandem modeling

Solving detailed physics: 1D detailed solver

Boundary conditions

Surface SRH

Non-contacted boundary

Ideal contact

TLS / MIS contact

MS case (Schottky contact)

MIS case (direct tunneling contact)

TLS case

Quasi-neutral boundary condition

Parameterizing detailed skins: skin solver

Multiscale solver

Coupling degree

Texture multiplier

Loss analysis

Background

Unique operating point

CLA: current loss analysis

Balances

avoidable and non-avoidable light-trapping losses

PLA: power loss analysis

FELA

Optical power losses

PLA equals efficiency loss

Material properties

List of symbols

References