*2.1. Leakage Modeling*

Since individual leakage mechanisms have distinct temperature dependencies, temperature dependent I-V behavior is obtained. Reverse biased diode characteristics over a

range of temperatures (T) from 50 ◦C to 130 ◦C are displayed in Figure 3a. The maximum cathode voltage (VCathode ) was limited to *VBR*/2 to avoid degrading the samples, and obtain clean trends with T for medium voltages. The presence of two different natures of variation with T is found, hence two regions were identified to be modelled separately.

**Figure 3.** Modeling of the reverse-biased characteristics of the p + -n diodes under test [55]. (**a**) Reverse diode characteristics from T = 50 ◦C to 130 ◦C. The two distinct regions identified in (**a**) are fitted using the Coulombic potential well model in (**c**) for VCathode from 0.5 V to 30 V (in direction of arrow), and using the variable range hopping model in (**d**) for VCathode from 70 V to 75 V (in direction of arrow). (**b**) Displays the good conformity of the fits with adjusted R <sup>2</sup> <sup>≈</sup>1 using the statistical parameter of adjusted R-square (coefficient of determination).

The first region, from VCathode = 0 V to 30 V, with a strong increase in current with temperature, was found to best represent conduction from Coulombic traps through thermionic emission [54,72]. The corresponding fit data is presented in Figure 3c. This mechanism is based on the assumption that the potential around traps at low electric fields can be considered Coulombic, while at higher fields, according to the Poole-Frenkel effect, a lowering of the potential barrier is expected with a square root dependency on field, strengthening the emission process of the trap [59,72–74]. This is expressed in the following formula, and the parameters are defined in [55]:

$$I\_{TE} = AT^2 \exp\left(\frac{-E\_A}{kT}\right) \tag{1}$$

$$e\_{\rm ll} \propto \exp\left(-\frac{E\_T - \beta F^{\frac{1}{2}}}{k\_B T}\right),\tag{2}$$

$$
\beta = \sqrt{\frac{q^3}{\pi \varepsilon'}} \tag{3}
$$

The slope extracted from the fitting (not shown) revealed an activation energy *E<sup>A</sup>* of ≈0.85 eV, usually associated with the presence of carbon acceptors [80,81], with an effective lowering in *E<sup>A</sup>* (∆*EA*) = 70 meV, the corresponding Poole-Frenkel coefficient *β* (=1.77 × 10 −5 eV V <sup>−</sup>1/2 m1/2 ) was found to be close to the theoretical value [55]

The second region, from VCathode = 70 V to 75 V, was modelled using variable range hopping (VRH), the leakage evolution fit to the VRH model is presented in Figure 3d. The corresponding equation is written as in Equation (4), and the parameters are described in [55]:

$$I\_{VRH} = I\_0 \exp\left[-1.76 \left(\frac{T\_0}{T}\right)^{\frac{1}{4}} + \mathcal{C}\_{VRH} \left(\frac{T\_0}{T}\right)^{\frac{2}{4}} F^2\right] \tag{4}$$

VRH describes the conduction of electrons across multiple trap states distributed within the bandgap. With the high occurrence of substantial defect densities in GaN epitaxial layers, VRH is commonly observed in GaN diodes [59,62–70], ascribed to the hopping of charged carriers through localized defect states in depletion regions.

For both the fits in Figure 3c,d, the adjusted R-Square (Adj. R-Square) [82] is found to be close to 1, as presented in Figure 3b, attesting to the good conformity of the fits. The R-square, also referred to as the coefficient of determination, always lies between 0 to 1, corresponding to whether the fit line is able to describe 0% or 100% of the variability of the data around the mean. Adj. R-Square is a modification which takes the number of predictors (within the fitted line) into account.

#### *2.2. Simulation of Doping Constraints in Diode Breakdown*

The investigation of breakdown issues is especially suited to using TCAD simulations, which provide versatile, non-destructive and rapid optimization solutions. A representative and simplified (fully vertical) model of the test devices was built using the Sentaurus tool from Synopsys in order to investigate the nature of breakdown, relative to the chosen concentration of p-doping in GaN diodes [55]. The drift diffusion transport model is used, along with appropriate polarization, mobility and recombination models. The n<sup>+</sup> layers are doped with N<sup>D</sup> = 5 <sup>×</sup> <sup>10</sup><sup>18</sup> cm−<sup>3</sup> , and the n− drift layer doping is fixed at <sup>N</sup><sup>D</sup> = 4 <sup>×</sup> <sup>10</sup><sup>16</sup> cm−<sup>3</sup> . For the p-body doping, Mg is defined as the dopant species. As discussed earlier, the Mg acceptors are not expected to be completely ionized at room temperature. Hence, to correctly estimate the effects of p-doping, using the incomplete ionization model is more physical. This model takes the parameters of the individual acceptor species into account, in particular, the ionization energy. Based on this, the simulator internally computes the effective doping concentration under different conditions. For example, a defined Mg concentration of N<sup>A</sup> = 6 <sup>×</sup> <sup>10</sup><sup>19</sup> cm−<sup>3</sup> with an ionization energy of 0.16 eV, leads to an effective base doping within the p-GaN region of <sup>≈</sup> <sup>4</sup> <sup>×</sup> <sup>10</sup><sup>18</sup> cm−<sup>3</sup> (6%), except within the depletion regions around the p-n junctions, where the defined N<sup>A</sup> is almost completely ionized.

Since the measured breakdown voltage of the test diodes is 170 V, the electric field evolution within the vertical diode is visualized at 160 V with different N<sup>A</sup> values in Figure 4. In Figure 4a,b, the chosen N<sup>A</sup> values are relatively low = 4 <sup>×</sup> <sup>10</sup><sup>17</sup> cm−<sup>3</sup> (see Figure 4a), 6 <sup>×</sup> <sup>10</sup><sup>17</sup> cm−<sup>3</sup> and 1 <sup>×</sup> <sup>10</sup><sup>18</sup> cm−<sup>3</sup> . In this scenario, the p-GaN region is observed to be severely depleted, with reach through occurring for the N<sup>A</sup> = 4 <sup>×</sup> <sup>10</sup><sup>17</sup> cm−<sup>3</sup> case, once the depletion regions from the n<sup>+</sup> -p and p-n− junctions intersect. Thus, a lower bound for setting the p-doping is identified owing to this constraint. In a real growth scenario, this constraint could be considerably tighter. If the reduction in Mg concentrations due to hydrogen passivation or other impurities were considered, the breakdown could occur faster (at lower voltages) for equivalent N<sup>A</sup> settings.

In Figure 4c,d the higher N<sup>A</sup> values are considered, including the representative value for the structures under test with N<sup>A</sup> = 6 <sup>×</sup> <sup>10</sup><sup>19</sup> cm−<sup>3</sup> (see Figure 4c). For these cases, the applied voltage drops almost entirely across the lightly doped n− GaN region, leading to smaller depletion of the p<sup>+</sup> GaN layer. On the other hand, the peak electric field at the p + to n− interface is significantly higher. In this scenario, breakdown in expected to be field-triggered, in fact, for N<sup>A</sup> = 6 <sup>×</sup> <sup>10</sup><sup>19</sup> cm−<sup>3</sup> , we are approaching critical field for GaN (≈ 3 MV/cm [83]) at the 160 V condition, which is found to agree reasonably well with the measured breakdown voltage of 170 V. Thus, the higher bound for N<sup>A</sup> settings is identified.

Based on the results in Section 2, we infer that the density of defects within the drift region need to be optimized to control the leakage current and its temperature sensitivity. The contribution of the residual carbon concentration is found to be relevant to the low voltage regimes, and needs to be optimized to improve the leakage performance. Regarding p-doping-induced constraints on the breakdown voltage, for a lightly doped drift layer, keeping the p-doping low can reduce the peak electric field, pushing *VBR* to higher voltages. However, the trade-off dictates that the value still needs to be high enough to avoid complete depletion of the p GaN layer unexpectedly at low voltages.

**Figure 4.** TCAD modeling of vertical p + -n diodes under different p-doping conditions describing the expected breakdown processes (**a**) TCAD structure visualized at N<sup>A</sup> = 4 × 10 17 cm−<sup>3</sup> ; (**b**) Electric field evolution for low p doping values illustrates complete depletion (punch-through) of the p-GaN region; (**c**) TCAD structure visualized at N<sup>A</sup> = 6 × 10 19 cm−<sup>3</sup> ; (**d**) Electric field evolution for high p doping values illustrates high electric fields (approaching critical field for GaN) at the p + -n − interface.
