A Closed-Form Model for 2-DEG Charge and Surface Potential in N-Polar Gallium Nitride Heterostructures

A closed-form model to calculate the 2-D electron gas (2-DEG) charge density and surface potential in an N-polar GaN/AlGaN heterostructure is presented. The proposed model is developed from the solution of Schrödinger’s equation for a finite triangular potential well. Asymptotic expressions of Airy functions are used to compute an initial value of subband energy. A new, unified expression to update the surface potential due to Poisson’s equation is presented. The modified Fang–Howard (M.F.H.) function is then employed to develop a closed-form correction to the subband energy using first-order perturbation theory. No fitting parameters are used. The results from our model are validated with numerical data for a wide range of bias conditions, channel layer thickness, and spacer layer mole fraction. Finally, the gate capacitance of the heterostructure is evaluated using automatic differentiation and validated against experimental data.

In a previous work [32], we developed an analytical model for the N-polar heterostructure.We noted that there is a significantly greater spread of the wave function into the spacer layer of the N-polar heterostructure, hence making the infinite confinement assumption inaccurate.Additionally, owing to the reversed polarization direction and resultant electric field, the subband energies are not directly related to the charge density.Thus, an analytical solution of the Schrödinger-Poisson equations in the N-polar orientation is more involved.
Our strategy in [32] was to describe the wave function in the channel, spacer, and barrier layers as linear combinations of Airy functions.The solution of a determinant equation, which enforces continuity conditions on the wavefunctions, yields very good solutions for the wavefunctions and an initial estimate of the subband energies.The final charge density, surface potential, and band diagram are obtained by applying first-order perturbation theory to incorporate the effect of charge density on the shape of the confining potential.The analytical approach presented above was, however, not a closed-form model.A closed-form approach is critical in the development of models for circuit simulation.
In this work, we present a closed-form model for the N-polar heterostructure by making three main modifications to our earlier work.First, we use the asymptotic forms of the Airy functions to convert the determinant equation into an explicit algebraic equation.Next, we develop an explicit formula to self-consistently update the charge and surface potential given the subband energy.Finally, we use the modified Fang-Howard (M.F.H.) wave function [33] to  [18] for the N-polar structure.Note that all energies are shown in volts rather than electron volts, for ease of presentation.Also, x denotes the Al mole fraction in Al x Ga (1−x) N.
simplify the evaluation of integrals involved in the solution of Poisson's equation and perturbation theory.We show that this closed-form model describes the charge density and surface potential very well, for a range of voltages spanning both the OFF and ON states, different mole fractions of the Al x Ga 1−x N spacer layer, and channel layer thickness.We also compute the gate capacitance from the closed-form charge expression using automatic differentiation and benchmark our results with experimental data.
This article is organized as follows.Section II describes the typical N-polar GaN heterostructure and reviews key aspects of our earlier, analytical model [32].The modifications required to convert the analytical model into a closed-form model are presented in Section III.Section IV discusses our results, Section V outlines the scope and limitations of our model, and Section V presents our conclusion.

II. STRUCTURE AND REVIEW OF ANALYTICAL MODEL
Consider the typical N-polar GaN heterostructure [17] shown in Fig. 1(b).Here, t i , t ch , t s , t b , t p , and t d refer to the thicknesses of the insulator, channel, spacer, ungraded barrier, graded barrier, and doped GaN layers, respectively.This structure includes both Si doping (N 1 , N 2 ) and polarization doping Q p due to compositional grading.From Fig. 1(c), we identify a charge [32]  Fig. 2 depicts the flow of our earlier analytical model [32].The blocks shaded yellow involve implicit/nonclosed-form expressions.In the present work, we present modifications such that one can obtain closed-form expressions for charge density and surface potential for a given V g .We have made available a Python implementation of our closed-form model in [34] to help avoid any ambiguities.
Our analytical model [32] begins by estimating an initial guess for the 2-DEG charge density n (0) s .We define V off as the gate voltage at which ε f = 0. We set n (0) s = 0 when V g < V off .For the case when V g ≥ V off , we compute n (0)  s by approximating the 2-DEG density as We then estimate an effective electric field (with magnitude E ch ) in the channel region, as a weighted average of the fields at the insulator-channel (E i|ch = σ (0)  m /ϵ ch ) and channel-spacer (E ch|s = σ (0) c /ϵ ch ) interfaces [32].This effective field approximates the confining potential in the channel region to be a triangular quantum well φ (0) (z).Given that the potential in the spacer and the barrier are also linear, we can, in general, represent the subband wavefunctions as linear combinations of Airy functions (Ai(Z ), Bi(Z )) in the channel, spacer, and barrier layers, as shown in Fig. 3(a) for an energy close to ε (0) 0 .Note that Z α (z) denotes a variable transformation to convert Schrödinger's equation into the Airy differential equation [32], [35] in region α ∈ {channel, spacer, barrier}.Continuity of the wave function and its derivative at the channel-spacer and spacer-barrier interfaces result in a 4 × 4 matrix problem.We also showed [32] that working with a reduced basis consisting of {Ai(Z ch ), Ai(Z s )} (and hence a 2 × 2 matrix problem) gave almost identical results for subband energies, due to the large difference in magnitudes of Ai(Z b ) and Bi(Z s ) with respect to Ai(Z ch ).The "2 × 2" method cannot predict the wave function in the barrier layer (since Ai(Z b ) is absent in the reduced basis), but gives identical wavefunctions as the "4 × 4" method in the channel and spacer regions.Note that "2 × 2" does not enforce that the wave function becomes zero at the spacer-barrier interface.Additionally, note a special case when the mole fractions of the spacer and barrier layers are identical, wherein the "2 × 2" method is exactly equivalent to the "4 × 4" method owing to the spacer-barrier interface being homogeneous.
Working with the "2 × 2" method, the wavefunctions in the channel and spacer layers are with where E s = ϵ ch E ch /ϵ s is the magnitude of the electric field in the spacer and m s(ch) and ϵ s(ch) being the effective mass and dielectric constant in the spacer (channel), respectively.a ch and a s are nonzero constants.The solution of Schrödinger's equation results in [A][n] = 0 with the Airy matrix Note that our choice of effective electric field provides a very good estimate for the wave function ξ 0 (z), without any further corrections, as shown in Fig. 3(b).The numerical results are obtained using a finite-difference-method-based Schrödinger-Poisson solver [30], [31].The numerical method enforces infinite-well boundary conditions over a sufficiently wide Schrödinger simulation region.While setting up the calculation, we ensure that the results are converged with (a) Airy Ai, Bi functions that serve as our basis to solve Schrödinger's equation.Note the large differences in magnitudes.The "2 × 2" method uses a reduced basis consisting of {Ai(Z ch ) and Ai(Z s )}.(b) First subband wavefunctions at V g = 0 V for spacer mole fraction x = 1.Note that ξ o (z) matches the numerical results very well, while ξ0 (z) is a good approximation for ξ o (z).respect to the size of this simulation region and the number of mesh points.
However, we need to correct the subband energy.For this, we update the charge density (n (1)  s ) and surface potential (or equivalently Fermi potential ε (1)  f ) by equating the charge density mandated by electrostatics with the charge density supported by the 2-DEG V th (7) setting j = 1.Here, d 0 represents the shift in the charge distribution away from the channel-spacer interface [32] [shown in Fig. 3(b)], given by where D = em ch /π h2 is the density of states in the channel, V th is the thermal voltage, and V go = V g − V off .This also provides a shifted, but still linear potential φ (1) (z).Note again that ( 6) and ( 7) present an implicit system of equations.
Poisson's equation is used to update the potential profile to include the effect of the shape of the charge distribution as with integration constants C ch = (σ c −en (1)  s )/ϵ ch and D ch chosen such that φ (2) (0) = ϕ (1)  s .The subband energy is corrected using first-order perturbation theory as ε (1)  0 = ε (0) 0 + ε 0 with and φ(z) = φ (2) (z) − φ (1) (z) being the perturbing potential.Note that the integrals in ( 9) and ( 10) can be solved analytically, but only in terms of Airy functions.Our closed-form model presented in the next section will eliminate the need for evaluating Airy functions.The final charge density n (2)   s and Fermi potential ε (2)  f are obtained using another self-consistent update ( j = 2).

III. CLOSED-FORM MODEL A. Asymptotic Forms of Airy Functions
We replace the Airy functions and derivatives in (5) using their lowest order asymptotic approximations for Z ch → −∞ and Z s → ∞ [36], [37], [38], evaluated when As an illustration, Fig. 4 compares the asymptotic formulae with the actual Airy functions and derivatives.For typical Npolar heterostructures, we are interested in the energy range ε 0 > 0.1 V. Clearly, within this energy range, the asymptotic formulae match very well with the exact expressions.These asymptotic functions simplify the determinant where 12) is an implicit equation in ε (0) 0 since both Z ch (0) and K are functions of ε (0) 0 .To convert ( 12) into a closed-form expression, we make use of the expressions available for the first energy level ε (0) 0,∞ in an infinite triangular well [29] ch to obtain ( 13) with Fig. 5 compares the initial estimate for energy ε (0) 0 obtained using the exact Airy matrix [see (5)], the asymptotic approximation [see (12)], and the closed-form expression [see (14)].Note that all these methods match quite well, justifying our approximations in deriving (14).

B. Explicit Update of Surface Potential
We now look at updating the surface potential (or equivalently, the Fermi potential) based on our estimate of subband energy.Since charge density depends on the relative positions of the subband energy level and Fermi level, this step also leads to an update in the value of charge density.From Fig. 2, the update of surface potential is performed twice in the flow of our model.Hence, for ease of presentation, we drop the superscripts on ε f , ε 0 , and n s in this section.
We first develop expressions for the weak charge, moderate charge, and strong charge regions.These regions are defined based on the relative values of the voltage V go with respect to the subband energy ε 0 and quantities α d = eDV th /C s and α n = α d • exp(1).Then, these expressions are combined to give a unified model for the Fermi potential which is valid across all the regions.Some broad aspects of our approach are similar to those available in the literature for Ga-polar heterostructures [19], [20].However, our approximations and resultant simplifications for the N-polar structure are very different from those for the Ga-polar structure.The results from our model are presented in Fig. 6.
1) Weak Charge Region (V go < ε 0 − 3V th ): Here, we assume that the charge in the channel is negligible.Eliminating ε f from ( 6) and ( 7), using en s /C s ≪ V go and approximating ln (1 + x) ≈ x gives Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
0 .Trace C presents true values, while traces A and B have been shifted to improve the clarity of presentation.Trace A demonstrates the validity of the region-wise approximations [see (17), (19), and (20)].Traces B and C demonstrate the interpolation forms ( 21) and (22).Note that the relative error with respect to an implicit solution [see ( 6) and ( 7)] is within 1%.(b) Auxiliary functions Γ a , Γ b that help to smoothly perform the interpolation in (22).
Substituting the above in (6) yields the following explicit function: In this region, we do not neglect the charge n s .Note, however, that the Fermi level lies below the lowest subband energy level (ε f < ε 0 ), as shown, for example, in Fig. 8(a) and (b).Hence, after eliminating the charge n s from ( 6) and ( 7), we make the same logarithmic approximation as in the weak charge region We further approximate ln(1 − x) ≈ −x−x 2 /2 and use ε f /V 2 go ≈ ε 0 /V 2 go to obtain the explicit expression 3) Strong Charge Region (V go > α n ): Here, the Fermi level lies significantly higher than the lowest subband energy level, as shown, for example, in Fig. 8(a) and (b).Thus, ( 7) simplifies into n s ≈ D(ε f − ε 0 ), which when combined with (6) yields The trace A in Fig. 6(a) compares our explicit expressions for the weak charge [see (17)], moderate charge [see (19)], and strong charge [see (20)] regions with a numerical solution of ( 6) and ( 7), evaluated with ε (1)  0 as the input.The good match justifies our approximations.Similar to [19], we unify the moderate charge and strong charge regions using interpolation In the moderate charge region, V go ≪ α n,d yields V gon,d → V go .Similarly, in the strong charge region, V go ≫ α n,d gives V gon → exp (1)/β and V god → 1/β.Furthermore, V th ε 0 /2V 2 go can be neglected in this region.Hence, ( 21) reduces smoothly to (19) or (20), as shown by trace B in Fig. 6(a).
We finally propose a single unified and explicit expression for ε f valid over the entire range of voltage (above and below V off ) as where δ = 0.01 V is a smoothing factor similar to that used in [19].When V go > 0, V gop ≈ V go given the small value of δ.However, when V go < 0, V gop is a voltage in the range of V th .Hence, the auxiliary functions a and b smoothly transition between b ≫ a in the weak charge region to a ≫ b in the moderate and strong charge regions, as shown in Fig. 6(b).
In the weak charge region, with the approximation ln(1+ x) ≈ x, we see that the unified expression (22) reduces to (17).Furthermore, in the moderate and strong charge regions, the approximation ln(1 + e x ) ≈ x reduces ( 22) to (21).Trace C in Fig. 6(a) shows that our unified model to obtain ε f (given the subband energy ε 0 ) matches the implicit numerical solution of ( 6) and ( 7) very well, with a relative error in the range of 1%.The charge density can be derived from ( 22) and (6) as

C. Fang-Howard Wave Function for Solution of Poisson's Equation and Perturbative Correction
To avoid using Airy functions in the integrals related to Poisson's equation ( 9) and perturbative correction (10),  (27) we approximate the wave function ξ 0 (z) using the M.F.H. function [39] with , and z 0 = 3/(κ + b) capturing an effective shift in the wave function.N is a dimensionless normalization constant given by The expression for b listed is approximate [40]; the exact variational problem with ξ 0 (z) cannot be solved analytically.However, since we already have estimates for E ch and ε (0) 0 , we can construct ξ 0 (z).As shown in Fig. 3(b), ξ 0 (z) approximates the Airy wave function ξ 0 (z) very well.The main advantage of using the M.F.H. function is in computing integrals involving wavefunctions, which can now be written in closed form.
Poisson's equation (9) for the updated potential φ (2) (z) is simplified using The perturbative contribution ε 0 | s from the spacer is very small compared to that from the channel ( ε 0 | ch ).Hence, we replace the upper limit of the integral in (10) by z = 0 and use the M.F.H. wave function ξ 0 (z) to obtain a closed-form expression for ε 0 = (D ch − ϕ (1)  s ) en (1)   s with θ k and k as listed in Table I.We have made use of the normalization condition on ξ 0 (z) while deriving the above expression.
f ) and numerical [31] data for different spacer mole fractions x.Note that ε 1 lies above the Fermi level.(c) and (d) Relative contributions of perturbative corrections from the channel and the spacer, justifying our approximation ignoring the latter.

IV. RESULTS AND DISCUSSION
Section III described the process of converting our analytical model [32] into a closed-form model.Both these models take advantage of the fact that our choice of the effective electric field E ch directly yields the correct wave function, without having to make any further corrections, as shown in Fig. 3(b).Hence, we can analytically estimate the shift d 0 in the charge distribution (away from the channel-spacer interface) using (8).A correct estimate of this shift helps to capture the effect of the shape of the wave function while describing the charge density mandated by the electrostatics in our structure [see (6)].Fig. 7 verifies that (8) correctly predicts the shift in 2-DEG charge distribution for different bias voltages, channel thickness, and spacer mole fractions.
Fig. 8(a) and (b) presents the final Fermi potential ε (2)  f and subband energy ε (1)  0 as a function of gate voltage, for different mole fractions of the spacer.Our closed-form results match well with numerical Schrödinger-Poisson [31] simulations.
Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.Fig. 9. Carrier density as a function of (a) gate voltage and (b) channel layer thickness showing a good match between the numerical data [31] and our model (23).Fig. 10.Gate capacitance as a function of (a) channel layer thickness and (b) gate voltage for different spacer layer mole fraction x.The experimental data (x = 0.3) is taken from [17].Results from our model are obtained using automatic differentiation.Also shown are results from numerical Schrödinger-Poisson simulations [31].
Note that the range ε (1)   0 > 0.1 V used in the asymptotic approximations [see (11)] of the Airy functions (see Fig. 4) is reasonable, allowing us to use the asymptotic forms to simplify the determinant equation.Also, the second subband ε 1 lies above the Fermi level justifying our approach of computing charge only from the first subband.The final value of the 2-DEG charge density n (2)   s obtained using our model is compared against numerical simulations [31] in Fig. 9.We find a good fit with numerical data for different values of channel thickness and spacer mole fraction.
Finally, the carrier density and Fermi potential models help us to determine the gate capacitance C g = edn (2)  s /d V g .Since we have a closed-form expression for charge density n (2)  s , the capacitance can be calculated using automatic differentiation [41].We use the Autograd library [42] to evaluate C g .We benchmark our results with experimental data for an N-polar stack with different channel thicknesses (but without a spacer) from [17], as shown in Fig. 10.We also include the results from numerical Schrödinger-Poisson [31] simulation for heterostructures with (x = 1) and without a spacer (x = 0.3).Our model matches exactly with the numerical results and also captures the experimental data very well.

V. SCOPE AND LIMITATION
In this section, we discuss the range of applicability of our closed-form model.With reference to Fig. 8  we see that the second subband comes closer to the Fermi level with an increase in gate voltage V g .We estimate that the charge from the second subband will grow once V g > V * g (where ε f = ε 1 − V th ) and become significant when V g > V # g (where ε f = ε 1 ).However, there is a voltage V gm0 (with V * g < V gm0 < V # g ) as shown in Fig. 11(a) where the charge on the gate σ m = 0.The direction of the electric field at the oxide-channel interface reverses around V gm0 , leading to a second quantum well at this interface, as shown in Fig. 11(b).When biased with V g > V gm0 , the device will suffer from enhanced gate leakage current [15] due to the second quantum well.Indeed, this could be the reason that experimental results for N-polar HEMTs (see [4], [7], [8], [9], [10], [11], [12], [13], [14]) rarely go beyond V g = 0 V. Finally, Fig. 11(c) shows the match between our closed-form model and numerical data for V g = 0.4 V. Hence, our model can comfortably be used for describing most of the experimentally measured data.

VI. CONCLUSION
We have presented a closed-form model for 2-DEG charge and surface potential in an N-polar GaN MIS-HEMT.This model is based on our earlier analytical model [31] obtained by solving the Schrödinger-Poisson equations in a finite triangular well and corrected using first-order perturbation theory.We convert this analytical model into a closed-form model by using asymptotic formulae for Airy functions, a unified interpolation formula for updating surface potential, and the Fang-Howard expression for wavefunctions while evaluating integrals.We have demonstrated that our surface potential model fits very well with the numerical results in all operating regimes.Finally, we have evaluated the gate capacitance using automatic differentiation and obtained very good agreement with experimental data.We believe that this work will be useful in the development of closed-form current-voltage models for circuit-level analysis involving N-polar MIS-HEMTs.

Fig. 1 .
Fig. 1.(a) Shortcomings in using models derived from Ga-polar heterostructures for a typical [17] N-polar GaN/AlGaN MIS-HEMT heterostructure as shown in (b).(c) Charge distribution, (d) conduction band profile, and (e) parameter values[18] for the N-polar structure.Note that all energies are shown in volts rather than electron volts, for ease of presentation.Also, x denotes the Al mole fraction in Al x Ga (1−x) N.
) such that σ m + en s − σ c = 0, with −σ m being the charge on the gate metal, n s the 2-D electron gas (2-DEG) charge density, and σ x being the polarization charge corresponding to the mole fraction x of the Al x Ga 1−x N layers.Fig.1(d) shows a portion of the band diagram for the structure in Fig. 1(b).ϕ s denotes the surface potential.All parameter values are listed in Fig. 1(e).

Fig. 2 .
Fig.2.Flowchart of our earlier, analytical model[32].The blocks shaded yellow include implicit/nonclosed-form expressions.Details are presented in Section II.The closed-from model presented in this work follows the same flow.
) and [n] = [a ch , a s ] T .The subband energy ε (0) 0 is obtained by solving the implicit equation setting the determinant |[A]| = 0.The constants a ch and a s are obtained as the null-space of [A].

Fig. 3 .
Fig. 3.(a) Airy Ai, Bi functions that serve as our basis to solve Schrödinger's equation.Note the large differences in magnitudes.The "2 × 2" method uses a reduced basis consisting of {Ai(Z ch ) and Ai(Z s )}.(b) First subband wavefunctions at V g = 0 V for spacer mole fraction x = 1.Note that ξ o (z) matches the numerical results very well, while ξ0 (z) is a good approximation for ξ o (z).

Fig. 4 .
Fig. 4. Airy functions, derivatives, and asymptotic approximations (evaluated at z = 0) in the (a) channel region and (b) spacer region.The asymptotic forms approximate the exact functions well in the observed range [see Fig. 8(a) and (b)] of subband energy ε (0) 0 > 0.1 V.

Fig. 5 .
Fig. 5. First estimate of the subband energy ε (0) 0 for spacer mole fraction (a) x = 0.3 and (b) x = 1 obtained using the Airy matrix, asymptotic formulae, infinite well and the closed-form expressions.Also shown are results from the numerical method [31], emphasizing the need for perturbative corrections to ε (0) 0 .

Fig. 7 .
Fig. 7. Shift in 2-DEG charge density away from the channel-spacer interface as a function of (a) gate voltage and (b) thickness of the channel layer for different mole fractions x of the spacer.Our model [see(8)] captures the numerical Schrödinger-Poisson results[31] well.

Fig. 11 .
Fig. 11.(a) Gate voltages V * g , V gm0 , and V # g for different channel layer thickness.(b) Band diagrams showing changing nature of the electric field at the insulator-channel interface.(c) Validity of our closed-form model at V g = 0.4 V. Numerical results are computed using [31].