Correct implementation of polarization constants in wurtzite materials and impact on III-nitrides

Accurate values for polarization discontinuities between pyroelectric materials are critical for understanding and designing the electronic properties of heterostructures. For wurtzite materials, the zincblende structure has been used in the literature as a reference to determine the effective spontaneous polarization constants. We show that, because the zincblende structure has a nonzero formal polarization, this method results in a spurious contribution to the spontaneous polarization differences between materials. In addition, we address the correct choice of"improper"versus"proper"piezoelectric constants. For the technologically important III-nitride materials GaN, AlN, and InN, we determine polarization discontinuities using a consistent reference based on the layered hexagonal structure and the correct choice of piezoelectric constants, and discuss the results in light of available experimental data.


I. INTRODUCTION
Pyroelectric materials have emerged in a variety of electronic and optoelectronic applications. Because of the symmetry of their crystal structure these materials exhibit spontaneous (SP) and piezoelectric (PZ) dipole moments [1] which manifest themselves as electric fields in heterostructure layers and sheet charges at interfaces. In the technologically important III-nitrides, which have the wurtzite (WZ) structure (space group P 6 3 mc), polarization differences allow for strong carrier confinement and the formation of a two-dimensional electron gas (2DEG) with high density at AlGaN/GaN interfaces, exploited in high electron mobility transistors (HEMTs). The effect of polarization can also be detrimental, for example causing the quantum-confined Stark effect in quantum wells of light-emitting diodes (LEDs), which reduces radiative recombination rates and shifts the emission wavelength. For both HEMTs and LEDs, accurate values of the SP and PZ polarization constants are required for a fundamental understanding as well as for device design.
Since experimental determination of the separate SP and PZ contributions to the total polarization is very difficult, calculated values of SP and PZ polarization constants are widely used in simulations. The PZ polarization constants are, in principle, fairly straightforward to explicitly measure or calculate [2]. However, the reported values exhibit a considerable spread [3]. In addition, the difference between so-called "proper" and "improper" PZ constants [2,4] is often overlooked, even though it can give rise to significant quantitative changes in the resulting polarization fields. This difference is one issue that is elucidated in the present paper.
The definition of SP polarization constants is even more subtle, and they are typically not amenable to explicit experimental determination, except in special cases [5]. The calculation of SP polarization requires the choice of a reference structure, which in the case of WZ semiconductors, has invariably been chosen to be zincblende (ZB) [6][7][8]. In this work we will show that though ZB as a reference structure is intuitively appealing, the SP polarization constants that result have been misinterpreted, introducing a source of error into the predicted values for bound sheet charge densities (and polarization fields). We also demonstrate that a proper choice of reference structure can eliminate these problems, and we provide revised values that can be directly inserted in current simulation tools.
While our theoretical considerations are general, we choose the nitride semiconductors because they provide a suitable example to illustrate the derivations, and because our findings have a significant impact on this materials system of high (and still increasing) technological importance. In Sec. II we review the underlying theory. In Sec. III we address the problems with choosing zincblende as a reference structure, and propose a solution. Sec. IV deals with piezoelectric contributions, specifically the issue of proper versus improper constants. In Sec. V we show that our findings have important consequences for nitride device structures and compare with previous implementations and with experiment. Section VI concludes the paper.

II. CALCULATING POLARIZATION CONSTANTS IN WURTZITE
For WZ films grown in the [0001] direction (i.e., the +c direction), the polarization component P 3 is given by the sum of the SP polarization at the wurtzite material's own lattice parameters, P SP , and the z component of the PZ polarization [1]. That is, for material m, where (in Voigt notation) ǫ i (i=1,2,3) is the strain in the i direction and e 3i are the corresponding piezoelectric constants (specifically, the "improper" ones; see Section IV). Henceforth we drop the subscript "3" from P for simplicity; all unbolded quantities pertaining to wurtzite are assumed to be in the c direction. The goal of this work is to derive the appropriate SP and PZ constants that allow Eq. (1) to be used in accurately determining polarization differences at interfaces between different WZ materials.
A. The modern theory of polarization Direct calculation of the polarization constants in Eq. (1) by first-principles electronic-structure methods was enabled by the formulation of a rigorous theory of bulk polarization, known as the modern theory of polarization (MTP) [9,10]. For a given structure λ, the MTP allows calculation of the so-called "formal" polarization [9]: where Ω is the cell volume, Z ion s is the charge of the ion s and R (λ) s is its position in the λ structure, f is the spin degeneracy of the bands, the sum j runs over occupied bands, and u (λ) j,k are the cell periodic parts of the Bloch wave functions. In Eq. (2), P el is the Berry phase taken over the valence-band manifold [9,10]. The formal polarization is defined only modulo the "quantum of polarization" eR/Ω, where R is any lattice constant and e is the electron charge [9,10].
In the MTP, only differences between formal polarizations of appropriate structures, λ = 0 and λ = 1, are well defined: The choice of the "appropriate" structures λ = 0 and λ = 1 rests on one of two possible considerations to ensure that physical conclusions can be drawn from their formal polarization differences. First, if the two structures are connected by an adiabatic, gap-preserving deformation path [9,10], then their difference in polarization [∆P in Eq. (3)] is given by the expression and corresponds to the zero-field adiabatic displacement current. This quantity can, in principle, be determined experimentally. An obvious application is the calculation of piezoelectric constants, which involves polarization differences between structures with slightly different lattice constants and/or internal structural parameters.

B. Interface theorem
The second consideration, as shown by , is that if an insulating interface can be constructed between two structures, the difference in formal polarization gives the bound charge, σ b , that builds up at the interface as a result of the continuity of the displacement field over an interface with no free charge: This is often referred to as the "interface theorem." Since there is no adiabatic path necessary between the two structures in this consideration, λ = 0 and λ = 1 can be different polymorphs of the same material (such as WZ and ZB structures of GaN) or different materials altogether (such as GaN and AlN); as long as they form an insulating interface, Eq. (5) will give the bound charge accumulation at the interface. From the interface theorem [Eq. (5)] and Eq. (1), the bound polarization charge at the interface between different III-nitride materials (m and n) is SP + e n 31 (ǫ n 1 + ǫ n 2 ) + e n 33 ǫ n 3 ] .
Note that σ b is the charge density of electrons at an interface for which material n has been grown on top of material m in the +c direction.

III. REFERENCE STRUCTURE FOR SPONTANEOUS POLARIZATION
A. Effective spontaneous polarization constants We will first address the difference in spontaneous polarization in Eq. (7), ∆P int SP . Strain effects will be taken into account separately in the PZ part, so that ∆P int SP is simply the difference of formal polarizations of the respective zero-strain structures, For purposes of Eq. (1), we would like to define a SP polarization constant that is a property of a single material. Simply taking P m f of Eq. (8) as P m SP is problematic, since formal polarization is multivalued, being only welldefined modulo a quantum of polarization eR/Ω. Therefore, in every situation in which Eq. (1) is applied to determine σ b at an interface, it must be confirmed that formal polarizations of the two materials are taken on the same "branch" of eR/Ω. A better approach is to take P m SP in Eq. (1) as a so-called "effective" SP polarization, P eff , defined by Resta and Vanderbilt [12] to be the ∆P in Eq. (3) that results as the system is taken from a high-symmetry "reference" structure (λ = 0) to the structure of interest (λ = 1). That is, Using P eff to define the SP polarization of the material removes the indeterminacy inherent to the formal polarization.
The reference structure is often chosen to be centrosymmetric, but it is important to recognize that the formal polarization of centrosymmetric crystals is not necessarily zero. This is because, as stated above, P f is a multivalued vector field, so it is possible for a nonzero formal polarization to be unchanged (modulo eR/Ω) under the inversion operator. Nevertheless, high symmetry puts restrictions on the possible values of P (λ=0) f [11]. While in principle effective polarization constants are still differences in formal polarization between λ = 1 and λ = 0 (reference) structures, in practice they can be used to compare spontaneous polarizations of different materials to obtain ∆P SP if such materials share a reference structure with the same formal polarization. Such a comparison then correctly yields the interface charge density according to the interface theorem of Ref. 11. In such cases, ∆P SP is just given by the difference in effective SP polarization of the materials, In the more general case that the reference formal polarizations do not match, the correct change in SP polarization following from Eqs. (8) and (9) is That is, a correction term of the form has to be added to Eq. (10). Unfortunately, this correction term is not typically implemented in device simulation packages (e.g., Ref. 13) or used in the interpretation of experimental data (e.g., Ref. 14, which is considered a standard reference in the field).
When the PZ terms are included as well, the total interface charge given by Eq. (7) becomes σ b = ∆ P int SP + ∆P ref corr − 2ǫ n 1 (e n 31 − e n 33 C n 13 /C n 33 ). (13) Equation (13) is a central result of the present work. B. Correction term for the effective spontaneous polarization with the zincblende reference structure As mentioned before, previous studies [6][7][8] have exclusively used ZB (space group F43m) as a reference structure for calculating the SP polarization of the WZ. This structure is not centrosymmetric, although it has sufficient symmetry to preclude any SP polarization [1]. The fact that an insulating (111) interface can be constructed between the WZ and ZB polytypes [8] makes it an appropriate reference structure. In fact, experimental measurements have deduced the relative polarization between the WZ and ZB phases of GaN [5], which were found to be consistent with the theoretical values in Refs. 6-8. However, there is a subtlety with using ZB as a reference structure: it has a nonzero formal polarization in the [111] direction, P ZB f (modulo eR/Ω). Again, this is consistent with the symmetry considerations because P f is a multivalued vector quantity, and can be nonzero while still remaining unchanged (modulo eR/Ω) under the F43m symmetry operations. These symmetry operations dictate the possible values of P ZB f , and therefore the resulting value depends only on the lattice constant, not on the chemical species of the atoms [11] (see Section S1 of the supplemental material (SM) [15] for more details on the formal polarization of ZB). The ZB reference structures for the reported effective SP polarization values for the III-nitrides were those with lattice constants equal to the in-plane lattice constant of the corresponding wurtzite material [6][7][8] (as confirmed by our calculations), so P ZB f will be different for GaN, AlN, and InN, and do not simply constitute a constant shift of P WZ f for all the materials. Therefore, for the effective SP polarization constants with the ZB reference to be implemented in Eq. (7) to determine the polarization difference between different WZ materials, the correction term of Eq. (12) is required, as in Eq. (13).
Consider the example of the interface charge between InN and GaN. Although, as mentioned above, the formal polarization of zincblende does not necessarily vanish, the symmetry of the structure severely limits the possible values. Specifically there are two possible values of the formal polarization in the [111] direction that are consistent with the symmetry: either P f vanishes, or it is equal to e √ 3/2a 2 n (both modulo e √ 2a n /Ω n ), where a n is the WZ in plane lattice constant of material n and Ω n is the volume of the ZB primitive cell (see Ref. 11 or Section S1 of the SM [15]). For the III-nitrides it is the latter, giving a correction term for GaN/InN: When considering the SP polarization differences between WZ nitrides, this represents a significant correction. In fact, as we will show in Section III C, the correc-tion is an order of magnitude larger than the effective polarizations when they are calculated with the zincblende reference [6,7]. As we shall see later in Section V E, this error is substantially reduced in practice by an approximate error cancellation that occurs in connection with the treatment of the PZ response.
There is nothing intrinsically wrong with using ZB as the reference structure for defining WZ effective SP polarization; however if these values are to be used to obtain polarization differences between different WZ materials, the ∆P ref corr term [Eq. (12), or Eq. (14) for the example of GaN/InN] must be explicitly included in expressions such as Eq. (11) or Eq. (13). To our knowledge, however, this has not been properly implemented in the numerous previous evaluations of SP polarization for nitride interfaces, and it would require changes in the software for the many simulation tools that include modeling of polarization fields in heterostructures.
C. P 63/mmc hexagonal layered structure as an alternative reference In order to avoid extensive changes in the simulation software, and to enhance physical insight, we advocate another approach, namely to determine effective SP polarization constants with respect to a reference structure for which the formal polarization is explicitly zero (so that ∆P int SP = ∆ P int SP ). A straightforward choice for this reference structure is the layered hexagonal (H) structure (space group P 6 3 /mmc), as was used for hexagonal P 6 3 mc ABC materials [16]. This structure is centrosymmetric, and we will show below with explicit firstprinciples calculations that it remains insulating and its formal polarization vanishes. The layered hexagonal structure can be obtained by an adiabatic (gap preserving) increase of the internal structural u parameter from u ≈ 0.37 − 0.38 of the WZ structure to u = 0.5. All that is required to avoid correction terms like Eq. (14) is to replace the effective SP polarization constants currently used in the field (the ones referenced to ZB [6,7]) with those referenced to the H reference structure. We have explicitly verified that this leads to expressions that are identical to those that would be obtained for the ZB reference, provided the second term in Eq. (11) or Eq. (13) is included.
The first-principles calculations of P f for the H, WZ, and ZB structures of the III-nitrides were performed using density functional theory with the screened hybrid functional of Heyd, Scuseria, and Ernzerhof (HSE) [17] as implemented in the Vasp code [18]. Hartree-Fock mixing parameters of 31 % for AlN and GaN, and 25 % for InN were used to correctly describe the band gaps and structural parameters of each material. Conventional functionals based on the local density approximation (LDA) or generalized gradient approximation (GGA) predict InN to be a metal, precluding the calculation of the polarization constants if the Γ point is included in the k- point mesh (which is required in Vasp). Projector augmented wave potentials (PAW) [19], with the In and Ga d electrons frozen in the core, were used. All calculations were performed on bulk primitive cells, with a 6 × 6 × 8 Monkhorst-Pack [20] k-point mesh to sample the Brillouin zone, and a large energy cutoff of 600 eV for the plane-wave basis set, chosen to ensure convergence of the internal structural parameter u. The calculated lattice parameters and band gaps, listed in Section S2 of the SM [15], show good agreement with experimental data.
We have calculated the electronic structure for structures with increasing u, ranging from u ≈ 0.37 to u = 0.5 ( Fig. 1), and confirmed that this path between WZ and H is gap preserving. These calculations also show that the formal polarization of the H structure is zero (modulo eR/Ω) for the III-nitrides (Fig. 1). We remind the reader that this was not guaranteed, since P f can be nonzero and still consistent with inversion symmetry, if the inversion operator changes P f by a multiple of eR/Ω. We have therefore verified that the hexagonal phase is a reference structure for which there is no spurious term in Eq. (11).
In addition, by correcting for any discontinuities (in the amount of a multiple of eR/Ω) that may occur in the calculations of formal polarizations along the path between WZ and H, we have insured that we are comparing formal polarizations of WZ GaN, AlN, InN on the same branch of eR/Ω [12].
The calculated spontaneous polarization coefficients for the WZ structure using either H or ZB as a reference are given in Table I. The results obtained by Bernardini et al. [6,7] are listed for comparison. The GGA functional used in that work provides results that are very close to those we obtained with HSE; the discrepancy is the largest for InN, which is probably related to the fact, mentioned above, that GGA predicts InN to be a metal. Table I also shows, however, that the choice of reference structure makes a significant difference. The magnitudes of the coefficients are much larger, and their signs are different when H is used as the reference. We observe that it is not just the absolute values, but also the relative differences between the calculated polarization constants of the three materials that differ from the previously reported values [6,7]. The difference in sign of P demonstrates that the conventional wisdom that the SP polarization in WZ points in the -c direction is misleading. The formal polarization of WZ has no definite sign as this would depend on the chosen branch. The effective SP depends on the polarization difference, and therefore the sign will depend on the sign and magnitude of the formal polarization of the reference structure. Even though the values reported in Table I  clearly differ in sign, absolute magnitude, and relative differences between materials, we will show in Section V that the final predictions based on both formulations are actually rather similar, because of the way the PZ contributions have been included in the previous work (cf. Section IV).

IV. IMPROPER VERSUS PROPER PIEZOELECTRIC CONSTANTS
We now address the specifics of the PZ terms in Eqs. (7) and (13). A complication that must be addressed is the choice between improper and proper e 31 (e 33 has no such complication) [2,4].
As we have done above, consider a thin layer of a WZ material grown in the c direction. If the layer is strained perpendicular to the c direction, the total bound charge on the +c and −c surfaces will change as a result of the polarization current, or redistribution of charge, in the layer. If metallic contacts on the +c and −c surfaces are short-circuited when the strain occurs, the current flow can be measured directly and will give the proper PZ constant, denoted e prop 31 [2,4].
If the +c and −c faces are in open-circuit boundary conditions, the layer will have a field across it due to the SP polarization, which will be modified by the strain via two mechanisms. The first is the same as in the proper case, as the strain will cause a flow of polarization current. But in addition, since the field depends on the charge density, the change in the area of the c-plane as a result of ǫ 1 will dilute or concentrate the pre-strain bound charge. For small strains the latter is given by the zero-strain formal polarization [2,4]. Taking both of these mechanisms into account gives the improper PZ constant, e imp 31 . In the case of, e.g., Eq. (13), the PZ constants correspond to the improper case, since their role in the equation is to take into account the change in formal polarization of material n with strain, so that σ b corresponds to the bound charge at the coherent interface with the in-plane lattice constant of material m. The change in formal polarization with strain is an alternative definition of the improper PZ constants [4].
where P n f ǫ=0 is the zero-strain formal polarization of material n. There is no change to the e 33 PZ constant. The proper PZ constant is a well-defined bulk quantity, as it is related to the polarization current; however the improper PZ constant is branch dependent [4]. Here also, defining polarization with respect to the H reference proves useful. Since the formal polarization of the H structure vanishes ( Fig. 1), P n,(H ref) eff = P n,WZ f ǫ=0 ; this also ensures that improper PZ constants for the different materials are taken on the same branch, in the same way as this is confirmed for the SP polarization constants. Therefore, consistent use of the H reference structure allows us to write Eq. (7) as where P

n,(H ref) eff
can be taken from Table I. Calculated proper PZ constants are given in the "proper" column of Table II. Since the HSE hybrid functional was used (and therefore density functional perturbation theory was not implemented), finite differences were used to calculate the derivatives with strain, following the procedure outlined in Eqs. (4)- (6) in Ref. 6. Specifically, improper PZ constants were calculated, and converted to proper constants by adding P n f ǫ=0 [see Eq. (15)] as determined in the calculation. This removes any dependence on the branch choice used in finite-difference calculations [4]. We then convert back to improper constants using P n,(H ref) eff as discussed above, in order to ensure that the constants are reported for the same branch for each material ("improper" column in Table II).
The WZ structure does have another nonzero piezoelectric constant e 15 , which couples a shear deformation in a plane perpendicular to the c plane (ǫ 13 or ǫ 31 ) to the polarization in the c plane. In this case, there are two improper PZ constants, since a correction must be included in the case of ǫ 13 but not for ǫ 31 [4]. These elements do not enter in the situation we consider in this work (plane stress conditions with the c plane as the growth plane), but may be important for growth on nonpolar or semipolar planes [21].
It is important to comment on the PZ constants reported in the literature [3]. When PZ polarization constants have been implemented in simulations (e.g. Refs. 3,13,and 14) it has never been specified which PZ constants are used for WZ III-nitrides. However, by comparing our calculations of proper and improper PZ constants ("proper" and "improper" columns of Table II) with the reported PZ constants in the literature ("prev. reported" of Table II) we have found that the reported constants are more likely to be the proper PZ constants.
From an experimental perspective, most of the experimental techniques have measured total polarization, and then deduced the PZ constants in the nitrides using the SP constants from, e.g., Ref. 6 using Eq. 1. As we will show in Section V, observation of the effects of total polarization can be misleading with regards to the differentiation between proper and improper PZ constants, due to the error cancellation from the use of the ZB reference in defining the SP polarization (discussed in Section III B).
There have been direct measurements of the PZ constants, either by probing the electromechanical coupling constants via surface acoustic waves [22,23], or by using interferometry to determine the strain caused by the application of a voltage [24][25][26][27]. Both of these techniques measure the proper constants, since neither is sensitive to the change in surface charge density resulting from the deformation. These reported values indeed agree well with our calculated values for the proper PZ constants, both in sign and in magnitude.
In previous work [7], P n,(ZB ref) eff was used instead of P n f ǫ=0 in Eq. (15) to convert improper to proper e 31 PZ constants (cf. where a n is the equilibrium, in-plane lattice constant of the WZ material n. To our knowledge the inclusion of the last term in Eq. (17) (17)], the distinction between proper and improper PZ constants is actually very significant.

A. Correct expressions for total polarization for wurtzite materials
Before discussing specific quantitative results for nitride semiconductors, we briefly summarize the main points of the previous sections and rigorously express the polarization of a given WZ material [Eq. (1)]. Spontaneous polarization constants must be defined with respect to a reference structure, and this choice of reference structure must be taken into account when evaluating polarization discontinuities at interfaces. We determined a correction term [Eq. (12)] that is necessary when effective SP polarization constants are used to determine the SP polarization difference between materials at an interface. This correction term is significant when the ZB reference structure is used [e.g., Eq. (14)], but is zero for the H reference. Using H as a reference is therefore more straightforward and is the approach we advocate, with the SP constants P  (Table II). Therefore, in the notation of this paper, Eq. (1) is written rigorously as P = P

B. Calculation of sheet charges for III-nitrides
Because of the important impact of polarization on device performance and design, a plethora of experimental studies have been aimed at determining the effects of polarization at GaN/InGaN and GaN/AlGaN heterostructures. We have plotted these reported results in Fig. 2, expressed as the magnitude of polarization bound charge at the interface, as a function of alloy content (a full list of references is provided in Section S3 of the SM [15]).
For GaN grown in the +c direction with the InGaN (AlGaN) grown on top, the sign of the bound charge at the interface will be negative (positive) [28].
In Fig. 2, the black dashed curves correspond to the current practice in the field: sheet charges are predicted based on (i) SP constants referenced to the ZB structure (P (ZB ref) eff in Table I) and Eq. (13) without the correction term ∆P

(ZB ref) corr
; and (ii) proper PZ constants ("proper" column in Table II). Quantities for alloys were obtained using linear interpolation. For an explicit expression in terms of alloy content, see Eq. (3) in Section S4 of the SM [15]. Elastic constants were taken from Ref. 29.
The red solid line in Fig. 2 corresponds to the implementation recommended in this work, i.e., using the H reference structure and the improper PZ constants, as in Eqs. (16) and (18) [and Eq. (5) in Section S4 of the SM [15]].
In view of the arguments given above, it may seem surprising that the dashed black and solid red curves agree as well as they do; we will return to this point in Section V E.

C. InGaN/GaN interfaces
For the InGaN/GaN system, most experimental studies have applied optical techniques to determine the polarization fields in GaN/InGaN/GaN quantum wells (QWs). This field can be probed by varying QW width [30] or external biases [31] and measuring the change in the optical properties of the QW (labeled "optical" in Fig. 2). In addition, there have also been studies using time-resolved PL to measure shifts due to screening of the polarization field by photoexcited or electrically injected carriers [32]. Other studies have been based on electron holography [33], where cross-sectional transmission electron microscopy is conducted on InGaN/GaN heterostructures to determine the depth-resolved electrostatic potential in the growth direction, and capacitancevoltage (CV) profiling of the fields [34]. When fields are reported, we convert to bound charge for the purposes of Fig. 2(a), assuming GaN/InGaN/GaN quantum well (with thick barriers such that the electric field in the barriers is presumed zero) using a simple parallel-plate capacitor model (E = σ/ε 0 ε r , using a relative dielectric constant for GaN of 10 [35] and for InN of 15 [36] and a linear interpolation for the dielectric constant of InGaN).
For specific values of the points in Fig. 2, see Section S3 of the SM [15].
The red curve in Fig. 2(a) is indeed in reasonable agreement with the experimental observations, appearing to be an upper bound of the data. The optical experiments usually rely on Schrödinger-Poisson simulations to determine the field magnitude from the measured optical properties. Uncertainties in input parameters to these models such as well widths, compositions, and composition profiles can result in quantitative differences. It has been shown recently that taking into account the deviations from an ideal QW structures when interpreting experimental observations can account for the apparent discrepancy between the measurements and theoretical prediction of polarization fields [37,38]. Such deviations are expected to be significant for InGaN/GaN because of the large lattice mismatch and the large difference in optimal growth temperatures for GaN and InGaN.

D. AlGaN/GaN interfaces
For the AlGaN/GaN system, there are two basic strategies for experimentally determining polarization effects. The first is to directly measure the polarization field in an AlGaN/GaN/AlGaN (QW) structure with the same methods as used in the InGaN/GaN case [39,40]. For the purposes of Fig. 2(b) we have converted these fields to bound sheet charge densities in an Al-GaN/GaN/AlGaN quantum well (using a relative dielectric constant for GaN of 10 [35]).
The other strategy is to measure the density of the 2DEG at the AlGaN/GaN interface in a HEMT structure (GaN channel, AlGaN barrier); from this, the bound interface charge, σ b , can be derived [41]. The 2DEG density can be determined either by Hall effect [14] or CV [41] measurements.
The significant scatter in the experimental data in Fig. 2(b) may have several origins. There are experimental uncertainties that can influence fields, such as incomplete strain relaxation in buffer layers [42], and differences in background doping [43][44][45].
As in the case of InGaN/GaN, the predicted sheet charges appear to be an overestimation compared to the experimental observations. In the case of optical measurements, Schrödinger-Poisson modeling is again typically used to interpret the measured properties, and the same uncertainties and systematic errors may arise as discussed in the case of InGaN/GaN [46,47]. For the cases where the compensating 2DEG density is measured, the thickness of the AlGaN layer and Schottky barrier height at the AlGaN surface will determine whether the entire bound charge is compensated, which could be a reason that the observations are slightly lower than predicted theoretically [48]. Also, interface roughness and electron traps due to dislocations and/or surface states have been proposed to explain the reduced 2DEG density [14].

E. Comparison of theoretical implementations
The degree of agreement between results obtained based on the current practice in the field (ZB reference, no correction term, proper PZ) and our revised implementation (H reference, improper PZ constants) merits some discussion. It is clear from Table I that the SP polarization values determined with the H reference structure are very different from those determined with the ZB reference; also, from Table II, the improper and proper e 31 coefficients are very different. However, the similarity between the red solid and black dashed curves in Fig. 2 demonstrates that these two large corrections cancel when we evaluate Eq. (18). I.e., the correct implementation (H reference, improper PZ constants) results in values that are only slightly different from the current (incorrect) practice in the field (ZB reference, no correction term, proper PZ). We now demonstrate analytically how this accidental agreement comes about.
If we take the difference between our revised implementation (solid red curve in Fig. 2) and the current practice in the field (dashed black curve in Fig. 2) for a given x, we obtain the total error in using the current practice in the field: A derivation of this expression is given in Section S5 of the SM [15]. For both AlGaN/GaN and InGaN/GaN, the two terms on the right-hand side of Eq. (19) have opposite signs and a tendency to cancel. To understand why, first note that ǫ 1 (x) is approximately linear in x, so that both terms can be regarded as being roughly proportional to ǫ 1 (x). In particular, linearizing Eq. (14) in ǫ 1 , we find x∆P [see also Eq. (8) of the SM [15]]. Thus where in the second step we used the fact that the formal polarization of the H structure vanishes. For the III-nitrides, |ǫ 1 | < 0.1, and we see from Table I  , is related to the difference in in-plane lattice constants between material n and material m (Section III B). The largest value it will take for the materials considered in this study is 0.28 C/m 2 for InN on GaN, as calculated in Eq. (14), and it will be significantly smaller for lower alloy content and for the case of AlGaN on GaN. The error is therefore the product of small factors, and thus small in practice, significantly smaller than the errors in the SP and PZ parts individually.
The small magnitude of P n,ZB eff (Table I) demonstrates that P WZ f ∼ P ZB f in these materials (see Section S1 of the SM [15]). The similarity between P ZB f and P WZ f is not unexpected. Although ZB has sufficient symmetry to preclude SP polarization, in the [111] direction the structure only differs from the WZ c direction by the stacking of the cation/anion planes and a small deviation from the ideal WZ u parameter and c/a lattice constants. In WZ materials other than the III-nitrides, such deviations could in principle be larger, and in that case P ZB ref eff will be larger, resulting in ∆P error being more significant.
For AlGaN/GaN, the relatively modest difference in lattice constants between AlN and GaN (and therefore modest ǫ 1 values for coherently strained alloy layers), and an almost exact cancellation between P GaN,ZB in a significant deviation at higher In content; this will be important for the prediction of polarization fields in applications such as tunnel field-effect transistors based on thin, high In-content interlayers [49].
The two implementations differ significantly in the relative contributions of SP and PZ polarization. An effect of this is illustrated by the case where there is strain relaxation in the alloy layer. In Fig. 3 the predicted polarization bound charges for In 0.2 Ga 0.8 N/GaN (blue curves) and Al 0.2 Ga 0.8 N/GaN (green curves) are shown as a function of strain relaxation of the layer, modeled by simply scaling ǫ 1 . For both InGaN/GaN and Al-GaN/GaN, the revised implementation of this work predicts a much faster decrease in bound charge at the interface than the current practice in the field. Of course, strain relaxation is associated with the presence of edge dislocations at the interface, which may themselves influ-ence the interface bound charge; this effect has not been taken into account in either version of the implementation.

VI. CONCLUSIONS
We have derived a rigorously correct implementation of polarization constants in wurtzite materials, focusing on the example of the III-nitrides. Our derivation has demonstrated the impact of the choice of reference structure when calculating spontaneous polarization constants using the modern theory of polarization. Insufficient care in using the values can result in spurious contributions to the polarization discontinuities at heterostructure interfaces. We have provided new values calculated with a consistent hexagonal (rather than zincblende) reference structure. In addition, we have demonstrated the importance of choosing the correct piezoelectric constants (improper), and provided values for these improper constants. These revised values of the spontaneous and piezoelectric constants can be directly used in simulations and to interpret experimental observations. The revised implementation predicts a more rapid decrease of polarization charge with strain relaxation for an alloy layer on GaN.

S1. POLARIZATION IN ZINCBLENDE
There is an additional subtlety when calculating the formal polarization of the zincblende (ZB) structure that is related to the choice of unit cell. If the primitive unit cell is used for the calculation [with the origin chosen such that the Ga atom is at (-1/8,-1/8,-1/8) and the N at (1/8,1/8,1/8)], then the result for III-nitirides will be what was determined in Ref. S1, Section III E. Specifically, the electronic part of the polarization vanishes, so the contribution simply comes from P ion of Eq. (2) of the main text: where the a WZ is the in plane lattice parameter of the wurzite ( 2a WZ (1/2,0,1/2)] that will result in a quantum of polarization eR/Ω that will take P ZB f to zero, and therefore ZB truly has a nonvanishing formal polarization. For the III-nitrides, we list the values of formal polarization for ZB along with WZ and the layered hexagonal (H) structure in Table I below.

S4. BOUND CHARGES AT NITRIDE INTERFACES
Here we present the specific equations used to generate Fig. 2 in the main text. As in the main text, we assume a coherent c plane interface of GaN and the alloy (InGaN or AlGaN), with the alloy layer under biaxial stress. The current practice in the field (black dashed curve in Fig. 2 of the main text) is to use the effective spontaneous (SP) polarization constants with respect to the zincblende (ZB) reference, without the correction term [∆P ref corr introduced in Eq. (11)  x + e GaN,prop where ∆ P where ∆P (ZB ref) and similarly for InGaN/GaN.