A new theoretical approach for the performance simulation of multijunction solar cells

A new theoretical approach is proposed for the performance simulation of multijunction (MJ) solar cells, starting from the weakness and strength of the Hovel model and of the transfer matrix method for describing the propagation of electromagnetic waves inside the solar cell structure. It is based on the scattering matrix method (SMM) and on a simplified generation function that allow describing with good accuracy the propagation of electromagnetic waves in the solar cell device, preserving, at the same time, the possibility of getting simple analytical solutions of the continuity equations. The numerical stability of the new theoretical approach is first demonstrated on triple junction InGaP/InGaAs/Ge solar cells, in which the Ge substrate is considered as the last layer (layer N) and then as the N‐1 layer. Further, the new theoretical approach is applied to simulate the performance of thin quadruple junction (QJ) InGaP/InGaAs/SiGeSn/Ge solar cells, in two‐ and three‐terminal configurations. Efficiency values of up to 45.1% and 44.9%, respectively, have been simulated at 1000× concentration, by considering the MJ limited by the InGaAs subcell. Finally, it is estimated that the QJ InGaP/InGaAs/SiGeSn/Ge solar cell has the potential to reach efficiencies over 50% by assuming proper antireflective coatings.


| INTRODUCTION
The simulation of the multijunction (MJ) solar cell performances requires several calculation steps. In Figure 1, for example, we show a possible scheme of the necessary steps. Among them, the choice of the generation function, which has to be inserted in the continuity equation to get the generation rates of holes/electrons, is of particular importance. According to the solar cell structures whose performances have to be simulated, different mathematical models have been implemented in order to calculate the propagation of radiation inside the solar cell. Finite element analysis has been applied to take into account 2D effects, 1 while rigorous coupled-wave analysis has been considered for solar cells with nanophotonic designs. 2 For MJ solar cells, the most applied mathematical models are the Hovel 3 and the transfer matrix method (TMM). [4][5][6] In particular, the analytical Hovel model is widely and successfully considered for calculating the performance of optically thick solar cell structures, while TMM has become the most applied mathematical method to evaluate the performances of MJ cell structures 7-10 and of some organic thin films, 11 where the thickness of the layers is comparable with the wavelength of radiation and where, therefore, interference effects have to be considered. However, if the advantage of TMM over the Hovel model is the ability to calculate more accurately the electromagnetic radiation inside the solar cell structure, its disadvantage is the higher complexity. In general, when the electromagnetic radiation in the solar cell layers is calculated by considering the interference of the forward and backward waves and inserted in the semiconductor continuity equation, it gives rise to solutions that are prohibitively complicated even for relatively few layers. [12][13][14] In particular, TMM presents an important weak point, as it suffers from numerical instability. For example, it has been reported that TMM fails in the determination of the reflection coefficient in the short wavelength range, when applied to the most recently designed MJ, in which the substrate is thinned and a mirror is deposited on the back side. 15 It would be therefore very useful for the performance simulation

| GENERAL ASSUMPTIONS
The mathematical models presented hereafter have all in common the following assumptions: 1. One dimension (edge effects are neglected).
2. Boltzmann approximation for the calculation of the carriers distribution. Nondegenerate semiconductors.
3. Low injection conditions. 4. All parameters, mobility, diffusion coefficient, and lifetime are assumed constant within each semiconductor material. Since the concentrations of dopant atoms in the base and emitter are assumed uniform (quasi-neutral region), the electric field value is zero outside the depletion region.
5. The mobility of minority carriers is assumed to depend on the doping of the side where majority carries are located. It should be pointed out that repulsive potentials scatter less than attractive potentials of the same strength 17 ; therefore, the minority carrier F I G U R E 1 Main steps for the calculation of the output performance of an MJ solar cell (as a particular case, a triple junction solar cell is considered) [Colour figure can be viewed at wileyonlinelibrary.com] mobility is higher than the majority carrier mobility. If the majority carrier mobility values obtained by Hall measurements are assumed for minority carriers, the related minority diffusion length values will be lower than the actual ones.
6. The solar cell window absorbs the solar light; however, the carriers here produced by absorption are neglected for the computation of the photovoltaic current. This assumption can lead to underestimating the photovoltaic current, as reported in Létay et al. 18 The same assumption is applied for the back surface field (BSF) layers. The reason behind these assumptions is that such layers are very thin compared with base and emitter layers of the solar cell.
7. The tunnel current peak value of each tunnel diode (TD) is supposed to be higher than the solar cell short-circuit current value.
Therefore, TD layers are considered only with respect to the light absorption. The series resistance of the TD is included in the total series resistance of the MJ device.
8. For each absorbed photon, an electron-hole pair is produced, that is, the photon absorption rate [n.ph s −1 m −3 nm −1 ] is equal to the carriers generation rate [n.e − /h + s −1 m −3 nm −1 ]. 9. Each subsolar cell of the MJ structure is divided into three parts: the two quasi-neutral regions of the emitter and the base and the space charge region. The continuity equations are solved in the quasi-neutral region to get the emitter and the base current densities. In the depletion region, owing the strong electric field, it is assumed that all the generated carries are swept in opposite directions and collected. 10. The generation terms owing to self-excitation 19 and photon coupling effects 20 are in first approximation neglected.

| HOVEL MODEL
The great advantage of the Hovel mathematical model is the simplicity of its generation function. In this model, multiple light reflections inside the photovoltaic active layers are neglected; therefore, the generation term, G, is calculated only considering a positive light flux, F (λ,x), absorbed inside the cell's structure according to the Beer-Lambert law. The reflection coefficient, R(λ), which takes into account the solar light reflected at the solar cells surface, is calculated only by considering the interference effects produced inside the coating and the solar cell window layer. By considering assumption 8 of Section 2, it holds the following:

| THE TMM
The TMM is a powerful mathematical method that takes into account the complex nature of the wave vector and allows considering the multiple reflections and interference phenomena of the forward and backward travelling waves. The calculated electromagnetic radiation by TMM can then be inserted in the semiconductor continuity equations to determine the photocarrier generation, replacing the flux term used by Hovel. A schematic representation of the MJ solar cell structure, generally composed of N layers, is reported in Figure 2, along with the possible coordinate systems.
The following general assumptions are usually considered: 1. the solar radiation propagates like a plane wave coming from the #0 layer (air); 2. the solar intensity is not polarized, and it can be thought as equally distributed in the TM and TE polarization modes; 3. the last layer (#N layer) is considered optically thick, and the reflected radiation at its back surface can be neglected; 4. each layer is homogeneous and isotropic; 5. source-free layers and interfaces (current density, J, and charge density, ϱ, are set to zero).
The electromagnetic radiation inside the MJ structure is calcu-  Table 1.
The solution of the linear equations system given by Equation Electromagnetic field at the beginning of layers j + 1, j Electromagnetic field at the end of layers j, j + 1 Notes. For the layer j, k j,x = x component of the complex wave number; n j complex refractive index; ε j = complex permittivity; μ j = complex permeability. In the S-matrix method, the outgoing waves are connected to the incoming waves. In particular, the layer scattering matrix S j is built up in such a way to connect the outgoing and incoming waves of the generic j layer with the outgoing and incoming waves of the last layer, N: The S j -matrix elements are calculated by forcing the compatibility of the system of Equation (13.3) with the system of Equation (5.1), therefore considering the following: Developing Equation (13.3), we have the following: while Equation (5.1) can be written as follows: By considering the above relationships, the compatibility of the system of Equation (5.3) with the system of Equation (5.4) allows calculating the scattering matrix elements: The outgoing and incoming waves in each layer, for TM and TE polarization, can be easily calculated, considering that from Equation (5.4)D, for j = 0, it holds the following: The electric field ℵ 0 + propagating in the layer "0" (air) in the positive direction (towards the solar cell surface) can be calculated from the flux of the Poynting vector of the solar radiation. Therefore, from Equation (5.8), we can get ℵ + N and, therefore, from Equation (5.4) and Equation (5.7) all ℵ + j and ℵ − j .

| A SIMPLIFIED GENERATION FUNCTION
The energy flux through each interface is calculated by means of the flux of the Poynting vector, for TM and TE modes. In the layer j, the Poynting vector at position r j is given by the following: The symbol < > means that we are considering the average of the Poynting vector, since the integration time of the measurement of the electromagnetic radiation is much longer than the period of the wave.
By considering that the Poynting vector does not change on the plane yz (owing to assumption 4 of Section 4), the energy crossing a surface in unit time and unit surface (the flux of the Poynting vector): where n is the normal is simply given by the following: It can be shown that whereS j x j À Á D E + represents the energy current density owing to the electric field propagating in the x positive direction,S j x j À Á D E − represents the energy current density owing to the electric field propagating in the x negative direction, whileS j x j À Á D E AE represents the energy current density oscillatory terms due to interference between the counter-propagating waves, whose expressions for TM and TE mode are reported in Table 2.
In order to get the generation function to be used to solve the continuity equation, we have to show how the Poynting vector is changing within each layer. For the moment, let us assume that we can neglect the energy current density oscillatory term due to the interference between the counter propagating waves (Equations 13.14 and 13.15); this assumption will be later assessed.
Considering the complex amplitudes of the forward and backward waves at the end of the layer, that is, in x j , and changing the position x towards the beginning of a layer (negative direction), it straightforward to show that The generation rate for the positive and negative fluxes is then given by the following: ð6:9Þ therefore, with x 0 < x < x N , it can be written as follows: T A B L E 2 Energy current density owing to the electric field propagating in the positive, negative direction and due to the counter propagating waves interference TM TẼ Notes. The symbol * indicates the complex conjugate; ω = angular frequency.
where x 0 has been assumed equal to zero without losing generality.
By considering assumption 4 of Section 2, the continuity equation can be written as follows: where D p is the hole diffusion coefficient and τ p is the hole lifetime. In Equation (7.2), the equilibrium hole concentration in the emitter, p n0 , has also been considered as position independent.
By setting

Equation (7.2) becomes
The generic solution of Equation (7.5) is as follows: The constants A n , B n can be determined by applying the usual two boundary conditions: dy dX = η p y at X = 0 ð7:7Þ with η p = SpLp Dp , where S p is the holes recombination velocity, and y = 0 at X = X N = X e − w n L p : ð7:8Þ From Equation (7.7), it follows: ð7:9Þ From Equation (7.8) and considering Equation (7.9), it holds the following:

ð7:10Þ
In the quasi-neutral emitter region, the current density is purely diffusive and can be calculated in X N , by the following: where A p and B p are given respectively by Equations (7.10) and (7.9). 7.2 | Determination of J n By considering Equation (6.10), the generation function for the base is as follows: By considering assumption 4 of Section 2, the continuity equation can be written as follows: where D n is the electrons diffusion coefficient and τ n is the electrons lifetime. In Equation (7.14), the equilibrium electron concentration in the base, n p0 , has also been considered position independent.
By setting y = n p −n p0 ð Þ L n 3 ; L n = ffiffiffiffiffiffiffiffiffiffi D n τ n p ;X = x− x p L n ð7:15Þ and The constants A n , B n can be determined by applying the two boundary conditions: where η n = SnLn Dn and S n is the electron recombination velocity. Considering Equation (7.19), it holds the following: Considering Equation (7.20), it holds the following: ð7:22Þ Setting The coefficient B n is given from Equation (7.22) by the following: In the quasi-neutral region of the base, the current density is purely diffusive, and it can be calculated in X=0: where B n is given by Equation (7.26).

| Determination of J dr
With the assumption that all the generated carries are swept in opposite directions and collected without losses (assumption 9 of Section 2), the current density generated in the depletion region is simply given by the following: where w = w n +Δ+w p (see Figure 3).

| THE TRIPLE JUNCTION INGAP/ INGAAS/GE CASE STUDY
A InGaP/InGaAs/Ge triple junction solar cell, whose structure is detailed in Figure 4, is considered as a case study to demonstrate that the photon flux oscillatory terms calculated from Equations  Table 3). This case study confirms that, with good precision, it is possible to solve the continuity equation The very low contribution of counter propagating waves in the overall photon flux is not a characteristic of specific solar cell structure: It has been also evidenced by Deparis, for a Si-based thin-film tandem solar cell. 23 Therefore, the proposed simplification of the generation function can be applied to any cell structure.

| ROBUSTNESS OF SMM VSTMM
In order to demonstrate the numerical stability of the SMM vs TMM, such mathematical methods have been applied to simulate the performance of a 3J InGaP/InGaAs/Ge solar cell in which the Ge substrate is first considered as the last layer (layer N) and then considered as the N-1 layer. The second case allows increasing the total thickness of the cell structure and then checking for any numerical instability. In the simulation, the current-voltage (I-V) curves have been obtained neglecting the series and shunt resistance (ideal I-V curves). Figure 6 shows that, when the Ge substrate is considered as the N-1 layer, the TMM fails in the determination of the reflectance for all the wavelength values <900 nm; on the other hand, the SMM is stable in all the wavelength range. The simulations of the TJ external quantum efficiency (EQE) and IV curve, carried out by applying SMM, give the same results regardless the Ge substrate is considered as the last layer (layer N) or as the N-1 layer (see Figure 7).   Table 4) The utilized grid masks for the 2T and 3T devices are reported in Figure 10. In Table 5, the series resistance calculations are reported, according to following formula 32 : where R sh is the InGaP or Ge emitter sheet resistance (window's contribution not considered); R c is the contact resistance; R m is the sheet resistance of the metal ρ m /w, ρ m and w are respectively resistivity F I G U R E 9 Three-terminal InGaP/InGaAs/SiGeSn/Ge four junction solar cells structure and related layer thicknesses. In the top cell, the BSF is constituted by a high P-doped InGaP layer, with a thickness of 10 nm. In the SiGeSn cell, the BSF is constituted by a high P-doped SiGeSn layer, with a thickness of 20 nm. In the figure, the top cell total thickness and the SiGeSn total thickness are indicated. The substrate has been removed [Colour figure can be viewed at wileyonlinelibrary.com] and thickness of the metal; a is the grid pitch; t is the grid width; b is the average grid length; S = t/a is the grid shadowing factor.
For the 3T device, the series resistance contribution has been calculated for the three junctions connected in series (top, middle, and intermediate) and for the Ge bottom junction. Note that in order to mitigate the series resistance problem of the 3T device, the cell area has been reduced to 1 mm 2 with respect to the value of 5.76 mm 2 related to the 2T device. Figure 11 shows the EQE of the 2T and of the 3T InGaP/InGaAs/SiGeSn/Ge QJ cell, while in Figure 12

| CONCLUSION
For developing next generation MJ solar cells, it is very useful to provide robust mathematical methods, which adopt simplified but sufficiently precise analytical solutions. Owing to its simplicity, the analytical Hovel model has been widely and successfully applied; however, it is not suitable in the performance simulation of QJ solar cells where the interference effects can be important. The most used TMM can present numerical instability when applied to MJ solar cells. Furthermore, the electromagnetic radiations calculated by this method, when inserted in the semiconductor continuity equations to determine the photocarrier generation gives rise to complicated solutions even for relatively few layers. In this work, a simplified SMM has been successfully implemented to overcome the TMM numerical instability, which arises when the MJ solar cell reaches a certain thickness. In order to demonstrate the robustness of the method, the SMM has been applied to TJ solar cell structures in which even a hundred micron thick layer has been included. We have also demonstrated that