TCAD Modeling and Simulation of Dark Current-Voltage Characteristics in High-Periodicity InGaN/GaN Multiple-Quantum-Wells (MQWs) Solar Cells

Gallium nitride (GaN) based high periodicity indium/GaN multiple quantum wells solar cells have recently been proposed for different applications: from operation in harsh environments, like space applications, to their use in concentrator solar harvesting system and wireless power transfer system. In this article, we extensively investigate the mechanisms, which influence the in-dark current-voltage characteristics of these device under no excitation by means of experimental measurements and numerical simulations performed with the technology computer-aided design (TCAD) Sentaurus suite from Synopsys. Two different nonlocal models are introduced: the first one considers intrabarrier tunneling, which was found to be the dominant conduction mechanism for voltages above the main diode turn-on; the second model considers trap-assisted tunneling, implemented through traps uniformly distributed in space and energy, to reproduce the current flow in the sub turn-on region. A good matching is obtained between the experimental and simulated electrical characteristics, by considering these nonlocal models in addition to the thermionic escape.


I. INTRODUCTION
G ALLIUM nitride (GaN) and its alloys with indium GaN (InGaN) and aluminum GaN (AlGaN) are wide band-gap semiconductors that find wide applications in high power/high frequency transistors and short-wavelength optoelectronic devices.In the last few years, GaN-based structures have been proposed as additional components of multijunction solar cells, in concentrator solar harvesting systems [1], [2], for wireless power transfer [3] and space applications [4], thanks to their reliability and efficiency in harsh environments [5].Early devices were based on the simple p-n and p-i-n approach, but recently InGaN/GaN multiple quantum wells (MQWs) solar cells started to be investigated.In particular, device performances have been studied as a function of the thickness of the barriers and of the wells [6], [7], [8], [9], or the periodicity of the quantum wells [10]; additionally, also the reliability of these device was evaluated under different conditions, such as light concentration and temperature [11], [12], [13], [14], [15].It is well known that both optical and electrical characteristic of MQW solar cells are very sensitive to the type and density of point defects and extended defects in the material, such as pure edge, pure screw, and mixed dislocations, and to the deep levels associated with these dislocations [16], [17], [18].With regard to the optical properties, the defects associated to threading dislocations may behave as Shockley-Read-Hall (SRH) recombination centers, decreasing the efficiency at low excitation regimes, where SRH recombination plays the strongest role [19], [20], [21].With regard to the electrical properties of these devices, these deep levels can favor trap-assisted tunneling (TAT) of carriers, which leads to an increase in the leakage current below the main diode turn-ON voltage [22], [23].However, a comprehensive modeling of the current-voltage characteristic of InGaN/GaN high periodicity MQWs solar cells in a wide current density range, from pA/cm 2 to mA/cm 2 , is still missing in literature.This model would be of extreme importance to understand the role of different conduction mechanisms, such as tunneling through barriers and trap-assisted tunneling through defects in the active regions [12], [13], [24], [25], [26].In this article, we reproduced the current-voltage characteristics of InGaN/GaN MQWs solar cells through numerical simulations performed with the technology computer-aided design (TCAD) Sentaurus suite from Synopsys Inc [27], [28].This goal was achieved by tuning the different conduction mechanisms that can influence device operation [12]: for instance, intrabarrier tunneling gives the most relevant contribution to current flow above the main diode turn-ON, especially for high current levels, whereas trapassisted tunneling mediated by impurities located in the active region gives the most relevant contribution at low current levels.
In this study, trap-assisted tunneling through traps located in the active region of the devices is also used to model the contribution of threading dislocations as leakage current pathways, which are especially significant for current flow under the main diode turn-ON [16].Considering these tunneling mechanisms, we were able to reproduce the in-dark current-voltage characteristics over a wide current range.

II. DEVICES UNDER TEST AND EXPERIMENTAL DETAILS
InGaN/GaN MQWs solar cells were grown on c-plane (0001) sapphire by conventional metal-organic chemical vapor deposition (additional growth details can be found in [29]).The devices under test, whose structure is presented in Fig. 1(a), consist of 2 μm n-GaN (Si doped, 3 × 10 18 cm −3 ) and a 125 nm n + -GaN (Si doped, 2 × 10 19 cm −3 ) layer.Above the n-type layer, an active region composed by a periodic structure of 30 pairs of undoped In 0.15 Ga 0.85 N quantum wells (3 nm) and GaN barriers (7 nm) was grown.On top of the active region, the devices feature a p-region consisting of a 5 nm p-Al 0.15 Ga 0.85 N electron blocking layer (EBL) (Mg doped, 2 × 10 19 cm −3 ), a 100 nm p-GaN layer (Mg doped, 2 × 10 19 cm −3 ), and a 10 nm p + -GaN contact layer (Mg doped, >2 × 10 19 cm −3 ).A semitransparent 130 nm indium-tin oxide (ITO) layer acting as transparent current spreading layer is deposited by dc-sputtering on top of the mesa with postannealing in N 2 /O 2 at 500 °C.Devices were then processed by standard lithography into 1 mm × 1 mm solar cells and Ti/Al/Ni/Au ring contacts and Ti/Pt/Au grid contacts were deposited via electron beam evaporation around the perimeter and on the top of the mesa, respectively, to form cathode and anode contacts.The performance, reliability, and efficiency of devices under different illumination conditions are discussed in [24], [25], [26], and [29].In Fig. 1(b), the energy band diagram at equilibrium is presented: details regarding the simulation framework are presented in the following section.In the band diagram, highlighted with a red shaded background, it is possible to note the 5 nm p-Al 0.15 Ga 0.85 N EBL, which introduces a barrier of 0.32 eV in the conduction band for electrons: this enables a more efficient extraction of holes from the p-GaN by blocking the overflowing electrons; thus, reducing the rate of holes loss due electron-holes recombination at the p side of the device [29].A HP 4145B semiconductor parameter analyzer is used to measure experimental I-V curves.

III. MODELING
Numerical simulations were carried out by TCAD Sentaurus suite from Synopsys Inc [27], [28].For the intentional doping, the different layers have been doped by placing magnesium (Mg) traps at the p side, and silicon (Si) traps at the n side, with the concentrations presented in the previous section and ionization energy levels reported in Table I: this allowed to obtain an overall activation of Mg dopants of 1%, which is compatible with an experimental finding of the dopant activation performed by the manufacturer of the cells [29].SRH, radiative and Auger generation-recombination mechanisms, and thermionic emission processes were considered with the main parameters reported in Table I.
IV. MODELING OF CONTACT AND PARASITIC RESISTANCES Fig. 2 reports the experimental current-voltage characteristic of a representative device (blue dotted triangles line).From the data it is possible to distinguish four different conduction regimes, also marked in Fig. 2.
1) A first region (R.1), up to 0.5 V, that is dominated by parasitic conduction mechanisms that can be modeled through a shunt component.2) A second region (R.2), between 0.5 V and 2 V, with an exponential behavior (linear slope in logarithmic scale) below diode turn-ON voltage.3) A third region (R.3), between 2 V and 2.5 V, corresponding to the onset of conduction of the main diode.4) A fourth region (R.4), above 2.5 V, in which the diode current is limited by the series resistance due to nonideal contacts, buffer layer resistivity, partial activation of doping in the quasi-neutral regions and other nonidealities [34].A first simulation was performed considering the abovedescribed structure with a shunt resistance of R P = 1 × 10 11 Ω and a series resistance of R S = 36 Ω, extrapolated from experimental measurements, and by enabling SRH, radiative and Auger generation-recombination, and thermionic escape processes.With regard to contacts, they were considered as ideal ohmic contact, since highly doped p and n layers were used to obtain a lower resistivity on the p-contact and n-contact of the device obtaining an ohmic contact [29], allowing to consider separately the values of the parasitic resistance.Finally, the ITO layer was not considered in the device, since it is introduced as current spreading layer [29] and, therefore, is considered as part of the ideal ohmic contact at the p-side of the device, whose resistivity is assumed to be position independent.From this simulation, the black curve of Fig. 2 was obtained, matching the turn-ON voltage and the slope of the experimental curve in R.3 and R.4, but the simulated current is at least 2 orders of magnitude lower than the experimental one, because of the nonconsideration of fundamental conduction processes for these devices, considered in the following paragraphs.

V. MODELING OF CARRIER TRANSPORT THROUGH THE MQW REGION
The low value of the forward current simulated without considering any additional transport/conduction models was ascribed to the presence of the 30 GaN barriers in the active region, which introduces an offset of 300 meV in the conduction band and of 180 meV in the valence band with respect to the quantum wells [12].Such offsets act as potential barriers that partially impede electrons and holes injection towards device regions of opposite polarity.As reported in other works, at room temperature thermionic emission alone is not sufficient to model the actual current flow occurring in real devices, therefore, intrabarrier tunneling must be considered [12], [35].The magnitude and relevance of this kind of process at a given bias voltage depends on the width of the barriers and of the wells, on the potential barrier height and also on the number of quantum barriers [9].In particular, the tunneling rate through a single barrier can be calculated and extended to the case of multiple (n) barrier obtaining [36] where L b is the barrier width, L w is the well width, m i is the effective mass in the well, m b,i is the effective mass in the barrier, and H i (F ) is the field-dependent effective barrier height for tunneling under an electric field F .Therefore, in addition to thermionic emission contribution, a nonlocal tunneling model has been included in our simulation framework: for each barrier in the active region, tunneling was implemented using the Wentzel-Kramers-Brillouin (WKB) tunneling probability model, according to [27], [28].By considering this additional transport model, the corresponding simulated device current (see Fig. 2, red solid line) increased by two orders of magnitude with respect to simulation without the intrabarrier tunneling model implemented (see Fig. 2, black solid line), matching the experimental values in R.3 and R.4.The only "free" parameter in this model is represented by the effective tunneling mass: this is a factor of merit that indicates the probability that tunneling takes place in a specific region of the device, i.e., in a specific quantum-defined environment, and using the WKB tunneling probability model in TCAD Sentaurus, it is defined with respect to the free electron mass m o in the material through a constant α: [27], [28].Fig. 3(a) reports the outcome of analysis of the sensitivity of the model with respect to the tunneling mass (the same value for hole and electrons effective tunneling mass was assumed, i.e., m e = m h ): with m e = m h = 0.3 ×m o the current increases by an order of magnitude with respect to the simulation without intrabarrier tunneling, and gains another two orders of magnitude when we set the tunneling mass to m e = m h = 0.003 × m o .These results confirm that tunneling through barriers represents the main transport mechanism at high forward voltages (see region 4 in Fig. 2) in p-i-n devices whose intrinsic region is composed by multiple QWs (30 in this case), also confirming the outcome of other works [35], [36].
In particular, to fit the Experimental I-V characteristics (see red curve in Fig. 2), we considered a tunneling mass of m  stated, the tunneling mass accounts for the efficiency of the tunneling process, whose value strongly depends on the local potential profile and, therefore, cannot be easily confirmed with experimental measurements; considering this, tunneling masses need to be adjusted [38], [39] according to the structure under investigation and to the implemented model.Nevertheless, in the case of simple direct tunneling through barriers, the effective tunneling mass is obtained with a good accuracy with respect to the effective mass of electrons and holes in GaN, highlighting the fit quality of the intraband tunneling conduction mechanism with tunneling mass values of the involved carriers close to the effective mass in GaN.Finally, in a first approximation, tunneling through barrier model does not depend on traps defined in the active region, but mainly depends on barrier thickness and number of multiple quantum well [6], [7], [8], [9].For the simulations reported in Fig. 3(a), we considered an equal effective tunneling mass for both electrons and holes.In fact, the simulated dark current-voltage characteristic is insensitive to the value of the hole tunneling mass for intraband tunneling, even if this parameter is varied by orders of magnitude.This can be seen in Fig. 3(b), where the black solid line represents simulations where the intrabarrier tunneling mass of electron is kept constant while varying that of the hole, and no changes are observed in the voltage-current characteristics.On the contrary, the dark I-V characteristic is very sensitive to the electron tunneling mass for intrabarrier tunneling, as shown in Fig. 3(b) where the dotted lines represent simulations where the intrabarrier tunneling mass of hole is kept constant while varying that of the electron: by decreasing m t,e by two orders of magnitude, we obtain a substantial increase in the current, especially at higher voltages, indicating that the transport of electrons through the MQW region represents a determining factor of the voltagecurrent characteristics.Analyzing the electron/hole densities at different operating voltages presented in Fig. 4(a), (b), and (c), the dominant role of electrons in modeling the forward device current is ascribed to a poor hole injection efficiency in the device [38].Fig. 4(a), (b), and (c), shows the hole (blue line) and electrons (red line) density concentrations at different operating voltages of 0.5, 2, and 2.5 V.In these graphs the energy band diagram is calculated at 0 V and is used as indicator of the position in the device, not being representative of the applied voltage.We can observe that electrons reach the other side of the junction much more quickly than holes, even if we considered an equal effective tunneling mass for electrons and holes: at 2.5 V, electrons reach the p side of the device (with density up to 1 × 10 14 ) whereas holes reach only half of the active region (with lower density also), as already observed in other previous characterizations [39] confirming a more efficient injection of electrons than of holes in the device [38].An effect of this poor hole injection can be experimentally seen by electroluminescence, where the QWs nearest to the p-side dominate the light emission as demonstrated in [39], [40], [41], and [42].The low hole injection efficiency arises from the polarization-related effects that exist in conventional InGaN/GaN devices grown on polar c-plane substrates [43], since the discontinuities in both spontaneous and piezoelectric polarization at the c-plane axis InGaN/GaN heterointerfaces generate large internal electric fields in the quantum wells (QWs), which may be detrimental for device operation [44], [45].For solar cells, the polarization field reduces the motion of carriers, lowering the collection efficiency, and therefore, the power conversion efficiency of solar cells [46], by decreasing the intrabarrier tunneling probability and in particular the hole injection efficiency [38], [47].Semipolar and nonpolar devices reduce the polarization-related electric fields originated inside the QWs, favoring hole injection efficiency and increasing tunnel probability, leading to a higher collection efficiency in these solar cells [48], maintaining high temperature performance [49], and favoring also indium incorporation [50].As can be observed in Fig. 4(d), the simulation confirms that by decreasing the polarization-related electric fields we obtain a higher hole injection efficiency, as just discussed.

VI. MODELING OF CURRENT CONDUCTION IN THE SUBTURN ON REGIME
The authors in [22] and [51] showed that the main mechanism determining the sub turn-ON leakage in forward bias is TAT.Through this process, electrons and holes tunnel, respectively, from the n-side and p-side of the device, toward defect states located in the forbidden bandgap of the undoped active region, resulting in a leakage current pathway.The model implemented for describing tunneling to traps consider the combination of two capture/emission mechanisms: a phonon-assisted inelastic process, whose electron capture rate from the conduction band is and an elastic transition [28], [52], [53] for which the electron capture rate from the conduction band is Here, the capture rates are presented for electrons but, similarly, can be found for holes: E T is the energy level of traps referred to the intrinsic Fermi level (E T = E t − E i ), V T is the interaction volume of the trap, S is the Huang-Rhys factor, ω the energy of the phonons involved in the transition (E phon ), m t is the electron tunneling mass, l is the number of the phonons emitted in the transition, and f is the Bose-Einstein occupation of the phonon state and z = 2S f (f + 1) and χ = √ l 2 + z 2 .The dissipated energy in this process is TAT can be described as an additional SRH-like recombination process, whose rate is given by Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
where N T is the trap density, g n,p are the electron and hole degeneracy factors, and c n,p are the capture rates described above.
The parameter set adopted for our simulation is summarized by Table II: we have set S = 10 and E phon = 91.2meV in agreement with previous literature results [22], [52].In the table there are two different value for the effective tunneling mass: one for the intraband tunneling, previously described in Section IV and one considered in the implementation of TAT.Such difference arises from the different kind of processes, as we will discuss later.In the following, the different tunneling masses will be referred to as effective intraband tunneling mass (m t,e/h ) or tunneling mass for the TAT model (m t,e/h,TAT ).
An initial consideration regards the relatively small value of the tunneling masses adopted for the modeling of TAT, m t,e,TAT = m t,h,TAT = 0.01, which originates from the complexity of the real physical process.As reported in several previous works on simulations regarding the case of single trap-assisted-tunneling [51], [54], [55], the values of the tunneling mass can differ a lot from the theoretical values, due to nonidealities, such as the difference between the real lattice configuration and the one assumed by the simulation model [56], [57], [58].Fig. 5 shows the band diagram at the equilibrium of the simulated device highlighting the traps considered for TAT modeling (see gray area in Fig. 5).Hole trap-assisted-tunneling has been implemented from each well in the active region, starting from the first well near the p-side, to the overall distributed in energy and space traps in the active region of the device, from left to right (see blue arrows in Fig. 5).Similarly, for electrons, trap-assisted-tunneling has been implemented from each well in the active region, starting from the first well near the n-side of the device from right to left (see red arrows in Fig. 5).
The total length of the region along which trap-assisted tunneling can take place correspond to the thickness of the active region, i.e., of 30 pairs of undoped In 0.15 Ga 0.85 N quantum wells (3 nm) and GaN barriers (7 nm).With regard to the traps involved in the tunneling process, donor traps have been placed over the entire active region and coupled through tunneling with the point of origin of the carriers: electrons towards the n-side and holes toward the p-side of the device, respectively.The donor traps have been uniformly distributed in both the energy and spatial domains.In particular, the energy window Fig. 5. TAT model implementation.Blue and red lines represent the conduction and valence band edges, respectively.Electrons (red circles) and holes (blue circles) tunnel into uniformly distributed trap states located close to midgap within the entire MQW region, respectively, from the right and left side of the junction, causing the forward TAT current.Blue and Red lines are representative of the nonlocal meshes defined along the entire active region for each quantum well, for holes and electrons, respectively.ranges from 0.5 eV above to 0.5 eV below the intrinsic Fermi level E i , for a total energy band of E sigma = 1 eV (this energy distribution is delimited in Fig. 5 by the gray lines) with an interaction volume of 6 × 10 −4 μm 3 and a concentration of traps of 5 × 10 15 cm −3 .Despite the energy distribution might seem unrealistic, former studies demonstrated that this kind of trap distribution is required to correctly reproduce the current voltage characteristic of similar devices [59].From a physical standpoint, this distribution can be explained by considering the presence of threading dislocations, i.e., of extended defects that can be associated with a relevant variety of trap states.In fact, while point defects can typically be associated to well discretized levels within the bandgap [59], [60], states related to extended defects, such as threading dislocations, have been found to cover wide portions of the forbidden gap [61] with deep levels located in a wide energy range [62], [63], [64], which can act as efficient nonradiative recombination centers [16], from which the uniformly distributed model could represent the presence of extended defects such as threading dislocations [59].The trap volume, i.e., the interaction volume of the uniformly distributed in space and in energy traps considered for the simulation, V T = 6 × 10 −4 μm 3 , is higher compared to previous works [23].This could be attributed to the presence of extended defects, which have higher interaction regime [18].
Finally, considering a typical concentration of threading dislocations (TD) of 10 8 −10 9 cm 2 for these GaN-based devices grown on sapphire substrate, the concentration of TD-related defect sites has been found to be in the order of 10 15 /10 16 [16], as properly considered in our simulation ( N T = 5 × 10 15 cm 3 , which was also a fitting parameter similar to the values chosen in simulations of similar devices [23], [65]) with typical value of the capture cross section (σ = 5 × 10 −15 cm 2 ) [65].The outcome of the simulations without (dashed red line) and with (solid green line) TAT model enabled, compared with the experimental I-V curve (red diamonds), is shown in Fig. 6(a): as can be noticed, trap-assisted tunneling through traps is the dominant mechanism below the main diode turn ON (at around 2 V-2.5 V) and strongly depends on the energy window defined in the active region, as shown in Fig. 6(b), as a larger energy window increases the probability that carriers will be able to tunnel through traps with lower energy jump.A significant current level before the main diode turn-ON is indicative of an I-V characteristic that is far from the ideality of the device, where trap-mediated conduction processes influence significantly the electrical characteristics of the device resulting in an increase of the subthreshold leakage current and reverse current [23], [65].The matching obtained over ten orders of magnitude of current confirms that intraband tunneling through barriers plays a dominant role in determining the current characteristic at high current level, above the main diode-turn ON, while below the main diode turn-ON, the I-V characteristic is mainly influenced by trap-assisted tunneling through defects located in the active region.Considering the devices under investigation, these defects can be associated to the presence of a high concentration of threading dislocations and related deep levels, which agrees with previous papers demonstrating that states close to mid-gap have a higher TAT probability in similar GaN-based device [19], [51].In our simulation, a model considering TAT mediated by TDs in the active region of the device allows to reproduce the current-voltage characteristic below the main diode turn ON, while tunneling through GaN barrier in the active region, in addition to thermionic emission from quantum barriers [12], [35], [66], represents the main mechanism responsible for current flow at higher bias levels [12], [35], [65].The proposed model is capable of reproducing the I-V characteristic over 10 orders of magnitude of current, thus allowing an effective study of relevant bias-dependent transport mechanisms on these high periodicity InGaN/GaN multiple quantum wells solar cells.

VII. CONCLUSION
Nonlocal tunneling Sentaurus TCAD models have been successfully implemented to reproduce the current-voltage characteristics of GaN-based with high periodicity InGaN/GaN multiple quantum wells solar cells in a current range exceeding ten orders of magnitude.A good matching with the experimental data was obtained through: 1) trap-assisted tunneling through uniformly distributed traps in energy and space located the active region, modeling the presence of high density extended defects such as threading dislocation.This allowed to reproduce the I-V curves below the main diode turn-ON and 2) intraband tunneling through barriers, that allowed to well reproduce the diode turn-ON and the current-voltage characteristics at high current level.This approach can be useful for understanding different transport mechanisms for high-periodicity solar cells, and the key role of defects in determining the electrical properties.

Fig. 3 .
Fig. 3. Blue triangles: experimental I-V curve.(a) Black curve: Simulated I-V curve without intra-barrier tunneling.Other solid curves: Simulated I-V curves with intrabarrier tunneling and effective tunneling masses ranging from 0.003 to 0.3 ×m o .(b) Black curve: Simulated I-V curve with intrabarrier tunneling keeping m t,e constant ( m t,e = 0.15 m o ) and varying m t,h from 0.01 to 0.15 ×m o (no variation is observed by changing m t,h , all curves are represented by the solid black one).Other dashed lines: Simulated I-V curve with intrabarrier tunneling keeping m t,h constant ( m t,h = 0.15 m o ) and varying m t,e from 0.01 to 0.15 ×m o (variation is observed by changing m t,e ).

Fig. 4 .
Fig. 4. Holes and electrons injection efficiency.Holes and electrons density at 0.5 V (a), 2.0 V (b) and 2.5 V (c).In these graphs (a, b, c) band diagram is calculated at 0 V bias and used as indicator of the position in the device.(d) Holes injection with different polarization-related electric field activation ratio at 2.5 V; here the band diagrams is calculated at 2.5 V.

Fig. 6 .
Fig. 6.Current-voltage characteristic.Blue triangles: Experimental I-V curve.Red dashed curve: Simulated I-V curve with intrabarrier tunneling model enabled.(a) Blue curve: Simulated I-V curve with intrabarrier tunneling model and TAT model simultaneously enabled.(b) Sensitivity to the energy window defined in the active region.

ACKNOWLEDGMENT
This study was developed in the framework of the research activities carried out within the Project "Network 4 Energy Sustainable Transition-NEST", Spoke 1., Project code PE0000021, funded under the National Recovery and Resilience Plan (NRRP), Mission 4, Component 2, Investment 1.3-Call for tender No. 1561 of 11.10.2022 of Ministero dell'Universitaà e della Ricerca (MUR); funded by the European Union-NextGenerationEU.This work was also supported by the Arizona State University and Rice University through ULTRA, an Energy Frontier Research Center (EFRC) funded by the U.S. Department of Energy, Office of Science, Basic Energy Sciences under Award DE SC0021230.The views and opinions expressed are, however, only those of the authors and do not necessarily reflect those of the European Union or the European Commission.Neither the European Union nor the European Commission can be held responsible for them.

TABLE I MOST
RELEVANT MATERIAL PARAMETERS ADOPTED FOR THE SIMULATION OF THE DEVICES

TABLE II MOST
SIGNIFICANT SIMULATION PARAMETERS: THE VALUE OF THE PARAMETERS WILL BE COMMENTED LATER IN THIS SECTION