%===========================%
% DOCUMENT SETTINGS         %
%===========================%
% draftversion = true makes %
%   single column           %
%===========================%
% draftversion = false makes%
%   double column           %
%===========================%
\def\draftversion{false}
%===========================%

% \documentclass[prx,twocolumn,showpacs,preprintnumbers,
%    citeautoscript,amsmath,amssymb]{revtex4-1}

\RequirePackage{ifthen}
\ifthenelse{\equal{\draftversion}{true}}{
  \documentclass[aps,prx,10pt,galley,amsmath,amssymb,showpacs,citeautoscript,
                 superscriptaddress]{revtex4}
}{
  \documentclass[aps,prx,10pt,twocolumn,amsmath,amssymb,showpacs,citeautoscript,
                 superscriptaddress]{revtex4-1}
}

%-------- MARGIN COMMENTS --------------
\def\scr{\scriptsize}
\ifthenelse{\equal{\draftversion}{true}}{
  \marginparwidth 2.7in
  \marginparsep 0.5in
  \newcounter{comm} % counter for commentaries
  % increase counter
  \def\commnext{\stepcounter{comm}}
  % commentary in text
  \def\commtext{{\bf\color{blue}[\arabic{comm}]}}
  % commentary in margin
  \def\commmar{{\bf\color{blue}[\arabic{comm}]}}
  % comment commands for all authors
  \def\dvm#1{\commnext\marginpar{\small DV\commmar: #1}\commtext}
  \def\cdm#1{\commnext\marginpar{\small CD\commmar: #1}\commtext}
  \def\cvm#1{\commnext\marginpar{\small CVdW\commmar: #1}\commtext}
  \def\mlab#1{\marginpar{\small\bf #1}}
  \def\tnewpage{\newpage\marginpar{\small Temporary newpage}}
}{
  \def\dvm#1{}
  \def\cdm#1{}
  \def\cvm#1{}
  \def\mlab#1{}
  \def\tnewpage{}
}
%----------------------------------

%-------- COLORS --------------
\usepackage[usenames]{color} % colors
\def\Red#1{\textcolor{red}{#1}}
\def\Blue#1{\textcolor{blue}{#1}}
\def\Magenta#1{\textcolor{magenta}{#1}}
%----------------------------------

%-------- STRIKEOUT ---------------
\usepackage{ulem}
%----------------------------------



\usepackage{graphicx}% Include figure files
\usepackage{dcolumn}% Align table columns on decimal point
\usepackage{bm}% bold math
\usepackage{epstopdf}
%
% \nofiles
%

\begin{document}

%
% \preprint{APS/123-QED}
%

% \title{Understanding and implementing polarization constants in
% wurtzite materials}
% \title{Spontaneous polarization constants in III-nitride materials
% revisited}
% \title{Correct implementation of spontaneous polarization constants
% in III-nitride materials}
%
\title{Correct implementation of polarization constants in wurtzite
materials and impact on III-nitrides}


\author{Cyrus E. Dreyer}
\affiliation{Materials Department, University of California,
Santa Barbara, CA 93106-5050}
\affiliation{Department of Physics and Astronomy, Rutgers
University, Piscataway, NJ 08845-0849}
\author{Anderson Janotti}
\altaffiliation{Current address: Materials Science and Engineering, University of Delaware, Newark, Delaware 19716-1501, USA}
\affiliation{Materials Department, University of California,
Santa Barbara, CA 93106-5050}
\author{Chris G. Van de Walle}
\affiliation{Materials Department, University of California,
Santa Barbara, CA 93106-5050}
\author{David Vanderbilt}
\affiliation{Department of Physics and Astronomy, Rutgers
University, Piscataway, NJ 08845-0849}

%
%  JOURNAL: prl?
%


\date{December 22, 2015}

\begin{abstract}

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,
%
% (space group $P6_3/mmc$) of these materials
%
and discuss the results in light of available experimental data.



\end{abstract}

\pacs{77.22.Ej,73.40.Lq,85.60.Bt}

%
% 73.40.Lq - Other semiconductor-to-semiconductor contacts, p-n
% junctions, and heterojunctions
% 77.22.Ej - Polarization and depolarization
% 85.60.Bt	-Optoelectronic device characterization, design, and modeling
% 81.05.Ea - III-V semiconductors
%

%
% \keywords{Suggested keywords}%Use showkeys class option if keyword
%
                              %display desired
\maketitle

\section{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 \cite{Nye1957} 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
$P6_3mc$), 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 \cite{Martin1972}.  However, the reported values
exhibit a considerable spread \cite{Feneberg2007}.  In addition,
the difference between so-called ``proper'' and ``improper'' PZ
constants \cite{Martin1972,Vanderbilt2000} 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. 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)
\cite{BFV1997,BFV2001,Bechstedt2000}.  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.~\ref{sec:theory}
we review the underlying theory.  In Sec.~\ref{sec:spont} we address
the problems with choosing zincblende as a reference structure,
and propose a solution.  Sec.~\ref{improper} deals with
piezoelectric contributions, specifically the issue of proper
versus improper constants.  In Sec.~\ref{exp} we show that our
findings have important
consequences for nitride device structures
and compare with previous implementations and with experiment.
Section~\ref{sec:conc} concludes the paper.

\section{Calculating polarization constants in wurtzite}
\label{sec:theory}
%
For WZ films grown in the [0001] direction (\textit{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_\text{SP}$, and the $z$ component of
the PZ polarization \cite{Nye1957}.
That is, for material $m$,
%
\begin{equation}
\label{totp}
\begin{split}
P_3^m=P_\text{SP}^m+(\varepsilon_{1}+\varepsilon_{2})e_{31}^m+\varepsilon_{3}e_{33}^m,
\end{split}
\end{equation}
%
where (in Voigt notation) $\varepsilon_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~\ref{improper}).
Henceforth we drop the subscript ``3''
from $P$ for simplicity; all unbolded quantities pertaining
to wurtzite are assumed to be in the $c$ direction.

\subsection{The modern theory of polarization}

Direct calculation of the polarization constants in
Eq.~(\ref{totp}) 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)
\cite{KSV1993,Resta1994}. For a given structure $\lambda$, the
MTP allows calculation of the so-called ``formal'' polarization
\cite{KSV1993}:
%
\begin{align}
\label{pf}
{\bf P}_\text{f} & ={\bf P}_\text{ion}+{\bf P}_\text{el} \nonumber \\
&=  \frac{e}{\Omega}\sum_s Z_s^\text{ion}{\bf
R}_s^{(\lambda)}+\frac{ief}{8\pi^3}\sum^\text{occ}_j\int_\text{BZ}
d{\bf k}\langle u_{j,\bf k}^{(\lambda)}\vert \nabla_{\bf k}\vert
u_{j,\bf k}^{(\lambda)}\rangle,
\raisetag{3\baselineskip}
\end{align}
%
where $\Omega$ is the cell volume, $Z_s^\text{ion}$ is the
charge of the ion $s$ and ${\bf R}_\text{s}^{(\lambda)}$ is its
position in the $\lambda$ structure, $f$ is the spin degeneracy
of the bands, the sum $j$ runs over occupied bands, and $u_{j,\bf
k}^{(\lambda)}$ are the cell periodic parts of the Bloch wave
functions. The formal polarization is defined only modulo
the ``quantum of polarization'' $e\textbf{R}/\Omega$, where
\textbf{R} is any lattice constant and $e$ is the electron charge
\cite{KSV1993,Resta1994}.

In the MTP, only {\it differences} between formal polarizations
of appropriate structures, $\lambda=0$ and $\lambda=1$, are
well defined:
%
\begin{equation}
\begin{split}
\label{delp}
\Delta{\bf P}={\bf P}_\text{f}^{(\lambda = 1)}-{\bf P}_\text{f}^{(\lambda = 0)}.
\end{split}
\end{equation}
%
The choice of the ``appropriate'' structures $\lambda = 0$
and $\lambda = 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
\cite{KSV1993,Resta1994}, then their difference in polarization
[$\Delta{\bf P}$ in Eq.~(\ref{delp})] is given by the expression
%
\begin{equation}
\begin{split}
\label{polcur}
\Delta{\bf P}=\int^1_0 d\lambda \frac{\partial {\bf P}}{\partial \lambda},
\end{split}
\end{equation}
%
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.

\subsection{Interface theorem}

The second consideration, as shown by Vanderbilt and King-Smith
\cite{VKS1993}, is that {\it if} an insulating interface can be
constructed between two structures, the difference in formal
polarization gives the bound charge, $\sigma_\text{b}$, that
builds up at the interface as a result of the continuity of the
displacement field over an interface with no free charge:
%
\begin{equation}
\label{disp}
\begin{split}
\sigma_\text{b}=({\bf P}_\text{f}^{\lambda=1} -{\bf
P}_\text{f}^{\lambda=0})\cdot {\bf \hat{n}}. %\equiv\Delta
%\textbf{P}_\text{int}\cdot{\bf \hat{n}}.
\end{split}
\end{equation}
This is often referred to as the ``interface theorem.'' Since
there is no adiabatic path necessary between the two structures in
this consideration, $\lambda=0$ and $\lambda=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.~(\ref{disp})
will give the bound charge accumulation at the interface.

From the interface theorem [Eq.~(\ref{disp})] and Eq.~(\ref{totp}),
the bound polarization charge at the interface between different
III-nitride materials ($m$ and $n$) is
%
\begin{equation}
\begin{split}
\sigma_\text{b} =\left[P_\text{SP}^m + e^m_{31}
 (\varepsilon_1^m+\varepsilon_2^m)+e^m_{33}\varepsilon^m_3\right]\\
- \left[P_\text{SP}^n + e^n_{31}(\varepsilon_1^n+\varepsilon^n_2)
 +e^n_{33}\varepsilon^n_3\right].
\end{split}
\end{equation}
%
As an example, we will take a realistic situation that occurs
in heterostructures, by assuming that material $n$ is strained
coherently to $m$ ($\epsilon_1^m=\epsilon_2^m=\epsilon_3^m=0$),
\textit{i.e.}, under plane stress ($\epsilon^n_1=\epsilon^n_2,
\epsilon^n_3= -2C^n_{13}/C^n_{33}\epsilon^n_1$, where $C_{ij}$
are the elastic constants). Therefore we have
%
\begin{equation}
\begin{split}
\label{bound}
\sigma_\text{b}&=(P_\text{SP}^m-P_\text{SP}^n )-
2\epsilon^n_1(e^n_{31}-e^n_{33}C^n_{13}/C^n_{33}) \\
&=\Delta P_\text{SP}^\text{int}-
2\epsilon^n_1(e^n_{31}-e^n_{33}C^n_{13}/C^n_{33}).
\end{split}
\end{equation}
%
Note that $\sigma_\text{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.

\section{Reference structure for spontaneous polarization}
\label{sec:spont}

\subsection{Effective spontaneous polarization constants}

We will first address the difference in spontaneous polarization
in Eq.~(\ref{bound}), $\Delta
P_\text{SP}^\text{int}$.
Strain effects will be taken into account separately in the
PZ part, so that
%
$\Delta P_\text{SP}^\text{int}$ is simply the
difference of formal polarizations of the respective
zero-strain structures,
%
\begin{equation}
\label{Pint}
\Delta P_\text{SP}^\text{int} =P^{m}_\text{f}\big |_{\varepsilon=0}
  -P^{n}_\text{f}\big |_{\varepsilon=0}.
\end{equation}

%CEDrev
For purposes of Eq.~(\ref{totp}), we would like to define a SP
polarization constant that is a property of a single material. Simply
taking $P^{m}_\text{f}$ of Eq.~(\ref{Pint}) as $P^m_{\text{SP}}$ is
problematic, since formal polarization is multivalued, being only
well-defined modulo a quantum of polarization
$e\textbf{R}/\Omega$. Therefore, in every situation in which
Eq.~(\ref{totp}) is applied to determine $\sigma_b$ at an interface,
it must be confirmed that formal polarizations of the two materials
are taken on the same ``branch'' of $e\textbf{R}/\Omega$. A better approach is to take
$P^m_{\text{SP}}$ in Eq.~(\ref{totp}) as a so-called
``effective'' SP polarization, ${\bf P}_\text{eff}$, defined by Resta
and Vanderbilt \cite{Resta2007} to be the $\Delta{\bf P}$ in
Eq.~(\ref{delp}) that results as the system is taken from a
high-symmetry ``reference'' structure ($\lambda=0$) to the structure
of interest ($\lambda=1$).
%
That is,
%
\begin{equation}
\label{Peff} \textbf{P}_\textrm{eff}
=\textbf{P}_\textrm{f}^{(\lambda=1)}
-\textbf{P}_\textrm{f}^{(\lambda=0)} =\textbf{P}_\textrm{f}
-\textbf{P}_\textrm{f}^\textrm{ref} .
\end{equation}
%
%(The term ``effective %polarization'' was actually introduced only
%later, by Resta and %Vanderbilt \cite{Resta2007}.)
Using
\textbf{P}$_{\text{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,
$\textbf{P}_\text{f}$ is a multivalued vector field,
% (as a result of the indeterminacy  modulo $e\textbf{R}/\Omega$),
%
so it is
possible for a nonzero formal polarization to be unchanged (modulo
$e\textbf{R}/\Omega$) under the inversion operator. Nevertheless,
high symmetry puts restrictions on the possible values of
${\bf P}^{(\lambda=0)}_\text{f}$ \cite{VKS1993}.

While in principle effective polarization constants are still
{\it differences} in formal polarization between $\lambda=1$
and $\lambda=0$ (reference) structures, in practice they
can be used to compare spontaneous polarizations of different
materials to obtain $\Delta P_\text{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.~\onlinecite{VKS1993}.
In such cases, $\Delta P_\text{SP}$ is just given by the difference
in \textit{effective} SP polarization of the materials,
%
\begin{equation}
\label{Pinttilde}
\Delta\widetilde{P}^\text{int}_\text{SP} =P^{m}_\text{eff}-P^{n}_\text{eff} .
\end{equation}
%

In the more general case that the reference formal polarizations do not match,
the correct change in SP polarization following from
Eqs.~(\ref{Pint}) and (\ref{Peff}) is
%
\begin{equation}
\label{Pint2}
\Delta P_\text{SP}^\text{int}=\Delta\widetilde{P}^\text{int}_\text{SP}
  +(P^{m,\text{ref}}_\text{f}-P^{n,\text{ref}}_\text{f}).
\end{equation}
%
That is, a correction term of the form
%
\begin{equation}
\label{Pintcorr}
\Delta P_\text{corr}^\text{ref}\equiv P^{m,\text{ref}}_\text{f}-P^{n,\text{ref}}_\text{f}
\end{equation}
%
has to be added to Eq.~(\ref{Pinttilde}). Unfortunately, this
correction term is not typically implemented in device simulation
packages (\textit{e.g.}, Ref.~\onlinecite{silense})
or used in the interpretation of experimental data (\textit{e.g.},
Ref.~\onlinecite{Ambacher2000}, which is considered a standard
reference in the field).

When the PZ terms are included as well, the total interface charge
given by Eq.~(\ref{bound}) becomes
%
\begin{equation}
\begin{split}
\label{boundcorr}
\sigma_\text{b}&= \Delta\widetilde{P}^\text{int}_\text{SP}
 +\Delta P_\text{corr}^\text{ref}
 - 2\epsilon^n_1(e^n_{31}-e^n_{33}C^n_{13}/C^n_{33}).
\end{split}
\end{equation}
%
Equation (\ref{boundcorr}) is a central result of the present work.




\subsection{Correction term for the effective spontaneous
polarization with the zincblende reference structure
\label{zbcorr}}

As mentioned before, previous studies
\cite{BFV1997,BFV2001,Bechstedt2000} have exclusively used
ZB (space group $F\bar{4}3m$) 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 \cite{Nye1957}. The fact that an insulating
(111) interface can be constructed between the WZ and ZB
polytypes \cite{Bechstedt2000} makes it
an appropriate reference structure.

However, there is a subtlety with using ZB as a reference
structure: it has a nonzero formal polarization in
the [111] direction, $P_\text{f}^\text{ZB}$ (modulo
$e\textbf{R}/\Omega$). Again, this is consistent with the
symmetry considerations because $P_\text{f}$ is a multivalued
vector quantity, and can be nonzero while still remaining
unchanged (modulo $e\textbf{R}/\Omega$) under the $F\bar{4}3m$
symmetry operations.  These symmetry operations dictate the
possible values of $P_\text{f}^\text{ZB}$, and therefore the
resulting value depends only on the lattice constant, not on the
chemical species of the atoms \cite{VKS1993}. 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 \cite{BFV1997,BFV2001,Bechstedt2000} (as confirmed by our
calculations), so $P_\text{f}^\text{ZB}$ will be different for each
material.  Therefore, for the effective SP polarization constants
with the ZB reference to be implemented in Eq.~(\ref{bound}),
the correction term of Eq.~(\ref{Pintcorr}) is required, as in
Eq.~(\ref{boundcorr}).

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_\text{f}$ vanishes, or it is equal to $e\sqrt{3}a_n/2\Omega_n$
(both modulo $e a_n/\Omega_n$), where $a_n$ is the lattice constant
of material $n$ and $\Omega_n$ is the volume \cite{VKS1993}. For
the III-nitrides our calculations result in the latter, giving
a correction term for GaN/InN:
%
\begin{equation}
\label{corr}
\begin{split}
\Delta P_\text{corr}^\text{(ZB ref)}&=P_\text{f}^\text{GaN,ZB}
  -P_\text{f}^\text{InN,ZB} \\
&=\frac{2e\sqrt{3}}{4}\left(\frac{1}{(a_\text{GaN})^2}
 -\frac{1}{(a_\text{InN})^2}\right) \\
&=0.28 \text{ C/m}^2.
\end{split}
\end{equation}
When considering the SP polarization differences between WZ nitrides,
this represents a significant correction. In fact, as we will show in
Section~\ref{hex}, the correction is an order of magnitude larger than
the effective polarizations when they are calculated with the
zincblende reference \cite{BFV1997,BFV2001}.  As we shall see later in
Section~\ref{sec:compar}, 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, as long as the $\Delta P_\text{corr}^\text{ref}$
term [Eq.~(\ref{Pintcorr}), or Eq.~(\ref{corr}) for the example
of GaN/InN] is explicitly included in expressions such as
Eq.~(\ref{Pint2}) or Eq.~(\ref{boundcorr}) . 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.

\subsection{$P6_3/mmc$ hexagonal layered structure as an
alternative reference \label{hex}}

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 $\Delta P^\text{int}_\text{SP}
=\Delta\widetilde{P}^\text{int}_\text{SP}$). A straightforward
choice for this reference structure is the layered hexagonal (H)
structure (space group $P6_3/mmc$), as was used for hexagonal
$P6_3mc$ $ABC$ materials \cite{Bennett2012}. This structure
is centrosymmetric, and we will show below with explicit
first-principles 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\approx 0.37-0.38$ of the WZ structure to
$u=0.5$. All that is required to avoid correction terms like
Eq.~(\ref{corr}) is to replace the effective SP polarization
constants currently used in the field (the ones referenced
to ZB \cite{BFV1997,BFV2001}) 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.~(\ref{Pint2})
or Eq.~(\ref{boundcorr}) is included.

The first-principles calculations of $P_\text{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) \cite{HSE06} as implemented in the
{\sc Vasp} code \cite{vasp1}.  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 $\Gamma$ point is included in the $k$-point mesh  (which is
required in {\sc Vasp}). Projector augmented wave potentials (PAW)
\cite{paw}, with the In and Ga $d$ electrons frozen in the core,
were used.  All calculations were performed on bulk primitive
cells, with a $6\times6\times8$ Monkhorst-Pack \cite{mp1976}
$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~S1
of the supplemental material (SM) \cite{Supp}, show good agreement
with experimental data.

We have calculated the electronic structure for structures
with increasing $u$, ranging from $u\approx 0.37$ to $u=0.5$
(Fig.~\ref{uint}), and confirmed that this path between WZ
and H is gap preserving by checking the electronic structure
for each calculation. These calculations also show that
the formal polarization of the H structure is zero (modulo
$e\textbf{R}/\Omega$) for the III-nitrides (Fig.~\ref{uint}). We
remind the reader that this was not guaranteed, since ${\bf
P}_\text{f}$ can be nonzero and still consistent with inversion
symmetry, if the inversion operator changes ${\bf P}_\text{f}$
by a multiple of $e\textbf{R}/\Omega$. We have therefore verified
that the hexagonal phase is a reference structure for which there
is no spurious term in Eq.~(\ref{Pint2}).

In addition, by correcting for any discontinuities (in the amount
of a multiple of $e\textbf{R}/\Omega$) 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 $e\textbf{R}/\Omega$
\cite{Resta2007}.

\begin{figure}
\includegraphics[width=\columnwidth]{Fig1.eps}
\caption{\label{uint}
Formal polarization of InN, GaN, and AlN for structures as a
function of the internal structural parameter {\it u}, varying
between fully relaxed WZ (circled symbols, labelled with relaxed
$u$ value) and H ($u=0.5$), as shown schematically by ball and
stick models where smaller balls represent N atoms and larger
ones represent the cations. All other lattice parameters were
fixed at their relaxed WZ values.}
\end{figure}

\begin{table}
\caption{\label{sptab}
Effective spontaneous polarization constants in units of C/m$^2$
of wurtzite (WZ) GaN, AlN, and InN calculated using either the
hexagonal (H, space group $P6_3/mmc$) or zincblende (ZB, space
group $F\bar{4}3m$) reference structures. The lattice constant of
the ZB structure is chosen to match the in-plane lattice constant
of the WZ structure for the same material.   Results from previous
calculations \cite{BFV2001} that used the ZB reference are listed
for comparison.}

\begin{ruledtabular}
\begin{tabular}{cccc}
   & $P_\text{eff}^{\text{(H ref)}}$   &
     $P_\text{eff}^{\text{(ZB ref)}}$  &
     $P_\text{eff}^{\text{(ZB ref)}}$, prev.\footnotemark[1] \\
\hline
     GaN & $ 1.312$   & $-0.035$     & $-0.034$ \\
     AlN & $1.351$     & $-0.090$       & $-0.090$ \\
     InN & $1.026$    & $-0.053$       & $-0.042$ \\
\end{tabular}
\end{ruledtabular}
\footnotetext[1]{From Ref.~\onlinecite{BFV2001}.}
\end{table}

The calculated spontaneous polarization coefficients for the
WZ structure using either H or ZB as a reference are given in
Table~\ref{sptab}.  The results obtained by Bernardini {\it et
al.} \cite{BFV2001,BFV1997} 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~\ref{sptab} also
shows, however, that the choice of reference structure makes a
significant difference.
%
% While the sign of all the spontaneous polarization constants remains the same,
%
%CEDrev
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
\cite{BFV2001,BFV1997}.

The difference in sign of $P_\text{eff}^{\text{(H ref)}}$ compared to
$P_\text{eff}^{\text{(ZB ref)}}$ 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. 
%When the PZ contribution is correctly included (Section~\ref{improper}), the physically relevant quantity in the case of the III-nitrides, the bound charge at the interface between different materials, has the same sign when either reference structure is used.


\section{Improper versus proper piezoelectric constants
\label{improper}}

%
% ************Imp only rewrite***************
%

We now address the specifics of the PZ terms in
Eqs.~(\ref{bound}) and (\ref{boundcorr}). A complication that must be addressed is the
choice between improper and proper $e_{31}$ ($e_{33}$ has no such
complication) \cite{Vanderbilt2000,Martin1972}.

%
% This is a result of the fact that strain perpendicular to the
% direction of SP polarization alters the bound charge at the
% surface/interface of a WZ layer via two mechanisms.
%

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_{31}^\text{prop}$
\cite{Vanderbilt2000,Martin1972}.

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 \textit{density}, the change in the area
of the $c$-plane as a result of $\varepsilon_1$
will dilute or concentrate the pre-strain bound
charge. For small strains the latter is given by the zero-strain
formal polarization \cite{Vanderbilt2000,Martin1972}. Taking both
of these mechanisms into account gives the improper PZ constant,
$e_{31}^\text{imp}$.

In the case of, \textit{e.g.},~Eq.~(\ref{boundcorr}), 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 $\sigma_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 \cite{Vanderbilt2000}.


%
% The first is that the total bound charge on the surface changes as
% the result of the strain, due to the polarization current that is
% caused by the strain; if the two sides of the layer are connected,
% this current can be measured directly and will give the proper PZ
% constant \cite{Vanderbilt2000,Martin1972}. However, the change
% in the field over the layer with strain (assuming the opposing
% surfaces/interfaces of the layer are not connected) depends not
% only on the total charge, but on the charge \textit{density}
% on the surface, which is also changing because the area of the
% layer is changing with strain. At small strains, this corresponds
% to the concentration or dilution of the bound charge resulting
% from the existing SP polarization \cite{Vanderbilt2000,Martin1972}.
%
%
% In the case of, i.e.~Eq.~\ref{boundcorr}, the PZ constants take
% into account the change in formal polarization of material $n$
% as a result of the requirement that it grow coherently on material
% $m$. The polarization difference between $m$ and $n$
%
%
% ***********************************************
%
%
%  **********Possible rewrite**********
%
%
% We now address the specifics of the PZ terms in
% Eq.~(\ref{bound}). A complication that must be addressed is the
% choice between improper and proper $e_{31}$ ($e_{33}$ has no
% such complication) \cite{Vanderbilt2000,Martin1972}. This is
% a result of the fact that for a finite layer of a WZ material
% (grown in the $c$ direction), the bound surface charge on the
% surface/interface terminated the layer is affected by in-plane
% strain (given by $e_{31}$) via two possible mechanisms.
%
% In the case where there is no free charge present and
% the opposite surfaces of the layer are not coupled (so no
% polarization current may flow between the surfaces), in-plane
% strain simply dilutes or concentrates the existing bound
% polarization charge, and the so-called {\it improper} constants
% \cite{Vanderbilt2000,Martin1972,Szabo1998} are appropriate to use:
%
%
% \begin{equation}
%
%
% \begin{split}
%
%
% \label{boundimp0}
%
%
% \sigma_\text{b,imp}=&\Delta\widetilde{P}^\text{int,(H
% ref)}_\text{SP}- 2\epsilon^n_1\Big({e}^{n,\text{imp}}_{31}
% -e^{n,\text{prop}}_{33}C^n_{13}/C^n_{33}\Big),
%
%
% \end{split}
%
%
% \end{equation}
%
%
% If new charge is allowed to arrive at the surfaces as the
% strain occurs, either from the opposite surface through a flow
% of polarization current, or an external source, the {\it proper}
% PZ constants should be used \cite{Vanderbilt2000}. In this case,
% the bound charge at an interface in is given by
%
%
% \begin{equation}
%
%
% \begin{split}
%
%
% \label{boundprop}
%
%
% \sigma_\text{b,prop}=&\Delta\widetilde{P}^\text{int,(H
% ref)}_\text{SP}- 2\epsilon^n_1\Big({e}^{n,\text{prop}}_{31}
% -e^{n,\text{prop}}_{33}C^n_{13}/C^n_{33}\Big),
%
%
% \end{split}
%
%
% \end{equation}
%
%
% *****************************************
%
%
% We now address the specifics of the PZ terms in
% Eq.~(\ref{bound}). A complication that must be addressed is the
% choice between improper and proper $e_{31}$ ($e_{33}$ has no such
% complication) \cite{Vanderbilt2000,Martin1972}. Specifically, in
% situations where the boundary conditions require a null electric
% field (as is the case for an infinite bulk crystal), the {\it
% proper} PZ constants (which are determined in a computational or
% experimental situation with null electric fields) should be used
% \cite{Vanderbilt2000}. In more practical terms, this corresponds
% to an experimental situation where any macroscopic fields are
% screened by free charge \cite{Martin1972}. The bound charge at
% an interface in this case would be given by
%
%
% \begin{equation}
%
%
% \begin{split}
%
%
% \label{boundprop}
%
%
% \sigma_\text{b,prop}=&\Delta\widetilde{P}^\text{int,(H ref)}_\text{SP}
% - 2\epsilon^n_1\Big({e}^{n,\text{prop}}_{31}
% -e^{n,\text{prop}}_{33}C^n_{13}/C^n_{33}\Big),
%
%
% \end{split}
%
%
% \end{equation}
%
%
% In contrast, situations where the depolarizing field is allowed
% to exist, such as in a thin quantum well with two insulating
% interfaces (no free charge present), require an extra term in
% the PZ contribution that comes from rotations and dilations of
% the zero-strain polarization moment \cite{Martin1972,Lax1973},
% resulting in the so-called {\it improper} constants
% \cite{Vanderbilt2000,Martin1972,Szabo1998}.
%

From Refs.~\onlinecite{Martin1972} and \onlinecite{Vanderbilt2000},
the improper
PZ constant $e^{n,\text{imp}}_{31}$ is related to the proper
constant by
%
\begin{equation}
\label{imp}
e^{n,\text{imp}}_{31}={e}^{n,\text{prop}}_{31}
  -P_\text{f}^n\big|_{\varepsilon=0},
\end{equation}
%
where $P_\text{f}^n\big|_{\varepsilon=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 \cite{Vanderbilt2000}. Here
also, defining polarization with respect to the H reference
proves useful. Since the formal polarization of the H
structure vanishes (Fig.~\ref{uint}), $P_\text{eff}^{n,\text{(H
ref)}}=P_\text{f}^{n,\text{WZ}}\big|_{\varepsilon=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.~(\ref{bound}) as
%
\begin{equation}
\begin{split}
\label{boundimp}
\sigma_\text{b}&=\Delta\widetilde{P}^\text{int,(H ref)}_\text{SP}\\
&- 2\epsilon^n_1\Big({e}^{n,\text{prop}}_{31}
-P_\text{eff}^{n,\text{(H ref)}}-e^{n,\text{prop}}_{33}C^n_{13}/C^n_{33}\Big),
\end{split}
\end{equation}
%
where $P_\text{eff}^{n,\text{(H ref)}}$ can be taken from Table \ref{sptab}.
\begin{table}[h]
\caption{\label{pztab}
Calculated piezoelectric polarization constants in units of C/m$^2$
compared with reported
%
%\Red{theoretical} \Blue{[Right?]}
values from the literature.
}
\begin{ruledtabular}
\begin{tabular}{c c c c c}
&	& proper	& improper 	& prev. reported\footnotemark[1] \\ \hline
GaN & e$_{31}$ 	& --0.551	& --1.863		&--0.22 to --0.55  \\
	& e$_{33}$ 	& 1.020 	& 1.020		&0.43 to 1.12 \\\hline
AlN 	& e$_{31}$ 	& --0.676 	& --2.027		&--0.38 to --0.81  \\
	& e$_{33}$ 	& 1.569 	& 1.569		&1.29 to 1.94 \\\hline
InN 	& e$_{31}$ 	&--0.604 	& --1.63		&--0.23 to --0.59  \\
	& e$_{33}$ 	&1.238 	& 1.238		& 0.39 to 1.09 \\
\end{tabular}
\end{ruledtabular}
\footnotetext[1]{From Ref.~\onlinecite{Feneberg2007} and references therein.}
\end{table}

Calculated proper PZ constants are given in the ``proper''
column of Table \ref{pztab}. 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.~\onlinecite{BFV1997}. Specifically,
improper PZ constants were calculated, and converted to proper
constants by adding $P_\text{f}^n\big|_{\varepsilon=0}$ [see
Eq.~(\ref{imp})] as determined in the calculation. This removes
any dependence on the branch choice used in finite-difference
calculations \cite{Vanderbilt2000}. We then convert \textit{back}
to improper constants using $P_\text{eff}^{n,\text{(H ref)}}$
as discussed above, in order to ensure that the constants are
reported for the same branch for each material (``improper''
column in Table \ref{pztab}).

It is important to comment on the PZ constants reported in the
literature \cite{Feneberg2007}.
When PZ polarization
constants have been implemented in simulations (\textit{e.g}.
Refs.~\onlinecite{silense,Ambacher2000,Feneberg2007}) it has
never been specified \textit{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
\ref{pztab}) with the  reported PZ constants in the literature
(``prev. reported'' of Table \ref{pztab}) we have found that the
reported constants are more likely to be the proper PZ constants.

In previous work \cite{BFV2001},
$P_\text{eff}^{n,\text{(ZB ref)}}$ was used instead
of $P_\text{f}^n\big|_{\varepsilon=0}$ in Eq.~(\ref{imp}) to
convert improper to proper $e_{31}$ PZ constants (cf. Table
VI and V of Ref.~\onlinecite{BFV2001}).  Because of the
nonvanishing formal polarization of the ZB reference
structure, $P_\text{eff}^{n,\text{(ZB ref)}} \neq
P_\text{f}^{n,\text{WZ}}\big|_{\varepsilon=0}$; instead, we
see from the discussion resulting in Eq.~(\ref{corr}) that
Eq.~(\ref{imp}) can be expressed as
%
\begin{equation}
\label{zbimp}
e^{n,\text{imp}}_{31}={e}^{n,\text{prop}}_{31}
 -\left(P_\text{eff}^{n,\text{ZB ref}}
 +\frac{e\sqrt{3}}{2a_n^2}\right),
%
% P_\text{eff}^{n,\text{(ZB ref)}}=P_\text{f}^{n,\text{WZ}}
%  \big|_{\varepsilon=0}-P_\text{f}^{n,\text{ZB}}\big|_{\varepsilon=0}
% = P_\text{f}^{n,\text{WZ}}\big|_{\varepsilon=0}
%
\end{equation}
%
where $a_n$ is the equilibrium lattice constant of material $n$. To
our knowledge the inclusion of the last term in Eq.~(\ref{zbimp})
has not been discussed in the literature.
%
% and therefore introduces a spurious dependence of the WZ PZ
% constants on the formal polarization of the ZB polytype of the
% material.

Because of the small magnitude of $P_\text{eff}^{n,\text{(ZB
ref)}}$, neglecting the last term in Eq.~(\ref{zbimp}) led to
the conclusion in Ref.~\onlinecite{BFV2001} that the difference
between the proper and improper PZ constants is small, seemingly
rendering the distinction of no consequence.
%
Instead,
because of the large magnitude of $P_\text{eff}^{n,\text{(H ref)}}$
[and $e\sqrt{3}/2a_n^2$ in Eq.~(\ref{zbimp})],
the distinction between proper and improper PZ constants is actually
very significant.

\section{Comparison with reported experimental results \label{exp}}

\subsection{Correct expressions for total polarization for wurtzite materials}
\label{review}

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.~(\ref{totp})]. 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.~(\ref{Pintcorr})] 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 [\textit{e.g.},
Eq.~(\ref{corr})], 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_\text{eff}^{\text{(H ref)}}$ listed
in Table \ref{pztab}.

In addition, the improper PZ constants should be used to determine
interface bound charge and fields in heterostructure layers. These
can be obtained from the proper constant $e^\text{prop}_{31}$
by subtracting $P_\text{eff}^{\text{(H ref)}}$ (Table
\ref{pztab}). Therefore, in the notation of this paper,
Eq.~(\ref{totp}) is written rigorously as
%
\begin{equation}
\label{Ptotfin}
P= P_\text{eff}^{\text{(H ref)}}+(\varepsilon_{1}+\varepsilon_{2})
  \left(e^\text{prop}_{31}-P_\text{eff}^\text{(H ref)}\right)
  +\varepsilon_{3}e^\text{prop}_{33}.
\end{equation}

%
% where we write $P_\text{imp}$ to indicate the situation with
% ideal insulating interfaces. Bound sheet charges are obtained
% using Eq.~(\ref{boundimp}).
%

%
% On the other hand, if the strain occurs in a situation where
% external free charge can flow to the surfaces/interfaces that
% terminate a layer, or polarization current can flow between the
% surfaces, the proper PZ constants are the relevant ones. Therefore,
% if we call the polarization in this situation $P_\text{prop}$,
% we have
%
%
% \begin{equation}
%
%
% \label{Ptotfin2}
%
%
% P_\text{prop}= P_\text{eff}^{\text{(H ref)}}+
%  (\varepsilon_{1}+\varepsilon_{2})e^\text{prop}_{31}
%  +\varepsilon_{3}e^\text{prop}_{33}.
%
%
% \end{equation}
%
%
% Bound sheet charges are obtained using Eq.~(\ref{boundprop}).
%

\subsection{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.~\ref{Both}, 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~S2 of
the SM \cite{Supp}).

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) \cite{Ambacher1999}.

In Fig.~\ref{Both}, 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_\text{eff}^{\text{(ZB ref)}}$ in Table \ref{sptab}) and
Eq.~(\ref{boundcorr}) \textit{without} the correction term $\Delta
P_\text{corr}^{\text{(ZB ref)}}$; and (ii) proper PZ constants
(``proper'' column in Table \ref{pztab}). Quantities for alloys
were obtained using linear interpolation (Vegard's law). For
an explicit expression in terms of alloy content, see Eq.~(1)
in Section~S3 of the SM \cite{Supp}. Elastic constants were
taken from Ref.~\onlinecite{Vurgaftman2003}.

The red solid line in Fig.~\ref{Both} corresponds to
the implementation recommended in this work, \textit{i.e.,} using the
H reference structure and the improper PZ constants, as in
Eqs.~(\ref{boundimp}) and (\ref{Ptotfin}) [and Eq.~(3) in Section~S3
of the SM \cite{Supp}].
%
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~\ref{sec:compar}.

%
% The green dot-dashed curve corresponds to an implementation
% relevant for situations where free charge screens the polarization
% field, \textit{i.e.}, using the H reference structure and
% the proper PZ constants, as in Eqs.~(\ref{boundprop})
% and (\ref{Ptotfin2}) [and Eq.~(3) in Section~S3 of
% Ref.~\onlinecite{Supp})].
%

\begin{figure}
\includegraphics[width=\columnwidth]{Fig2.eps}
\caption{\label{Both}
Absolute values for polarization sheet charges at the (a)
InGaN/GaN and (b) AlGaN/GaN interface as a function of alloy
content predicted from the spontaneous polarization constants
calculated using either the zincblende (ZB) reference structure
[without correction term, Eq.~(\ref{corr})] and the proper
piezoelectric constants (black dashed curve), or the hexagonal (H)
reference structure and improper piezoelectric constants (red solid
curve). Points are experimental values from the literature (see
section S2 of the SM \cite{Supp} for references and values). }
\end{figure}



\subsection{InGaN/GaN interfaces\label{ingan}}

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 \cite{Hangleiter2003} or external biases
\cite{Renner2002} and measuring the change in the optical
properties of the QW (labeled ``optical'' in Fig.~\ref{Both}).
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
\cite{Turchinovich2003}. Other studies have been based on
electron holography \cite{Cherns1999}, where cross-sectional
transmission electron microscopy is conducted on InGaN/GaN
heterostructures to determine the depth-resolved electrostatic
potential in the growth direction, and capacitance-voltage
(CV) profiling of the fields \cite{Zhang2004}. When fields
are reported, we convert to bound charge for the purposes of
Fig.~\ref{Both}(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=\sigma/\epsilon_0\epsilon_r$, using a relative dielectric
constant for GaN of 10 \cite{Barker1973} and for InN of 15
\cite{Madelung} and a linear interpolation for the dielectric
constant of InGaN). For specific values of the points in
Fig.~\ref{Both}, see Section~S2 of the SM \cite{Supp}.

%
% Since the measurements in Fig.~\ref{Both}(a)  were made on thin
% with nominally no free charge present, quantum wells, the improper
% PZ constants should be used.
%

The red curve in Fig.~\ref{Both}(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\"{o}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 \cite{Wu2012,Mcbride2014} can result in
quantitative differences.  Deviations from an ideal QW structure
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.

\subsection{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 \cite{Grandjean1999,Leroux1999}. For the purposes of
Fig.~\ref{Both}(b) we have converted these fields to bound sheet
charge densities in an AlGaN/GaN/AlGaN quantum well (using a
relative dielectric constant for GaN of 10 \cite{Barker1973}).
%
% We see that these optical measurements tend to be on the lower
% end of the experimental scatter.
%

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, $\sigma_\text{b}$,
can be derived \cite{Yu1997}.  The 2DEG density can be determined
either by Hall effect \cite{Ambacher2000} or CV \cite{Yu1997}
measurements.

Again, the predicted sheet charges are in reasonable agreement
with experimental observations. The significant scatter in
the experimental data in Fig.~\ref{Both}(b) may have several
origins. There are experimental uncertainties that can influence
fields, such as  incomplete strain relaxation in buffer layers
\cite{Zhou2002}, and differences in background doping. For
optical measurements, Schr\"{o}dinger-Poisson modeling is again
typically used to interpret the measured properties, and the
same uncertainties arise as discussed in the case of InGaN/GaN
\cite{Arulkumaran2003,Mkhoyan2004}.

%
% We see from Fig.~\ref{Both}(b) that all of the largest sheet
% charge densities were derived from Hall measurements; this may
% be an overestimation of the bound charge due to the fact that
% Hall measurements are affected by parallel conduction pathways
% in underlying layers \cite{Kohler2011,Davidsson2004,Simon2001},
% especially as background $n$-type doping in bulk GaN and AlGaN
% can be as high as $\sim$10$^{18}$cm$^{-3}$.
%

%
% Also, We find that for the case of AlGaN/GaN, the predicted
% polarization sheet charges underestimate the experimental values
% reported in the literature. Though the source of this descrepancy
% is not clear may be a result of the sources of scatter just
% mentioned. For example, strain relaxation will result in larger
% predicted sheet charges for a given Al content; the green,
% dot-dash curve in Fig.~\ref{Both} represents the case of full
% strain relaxation of the alloy layer. The disagreement between
% the reported experimental and predicted sheet charges merits
% further investigation.
%

\subsection{Comparison of theoretical implementations}
\label{sec:compar}

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.  The success of current practice in the field
has hinged on a fortuitous cancellation of errors between the spurious
term from the use of the ZB reference for the SP polarization and the
proper PZ constants instead of the improper ones. Comparing the
revised implementation [Eq.~(3) of the SM \cite{Supp}] to the current
practice in the field [Eq.~(1) of the SM \cite{Supp}] gives the
difference between the solid red and dashed black curves for a given $x$ in
Fig.~\ref{Both} as $ x \Delta P_{\text{corr}}^{\text{(ZB
ref)}}+2\varepsilon_1(x)P_{\text{eff}}^{n, \text{(H ref)}}(x)$ (see Section~S4
of the SM \cite{Supp} for a detailed derivation).  For both AlGaN/GaN
and InGaN/GaN, $\Delta P_{\text{corr}}^{\text{(ZB ref)}}$
and $2\varepsilon_1P_{\text{eff}}^{\text{(H ref)}}$ have opposite
signs and a tendency to cancel.

To understand why, note that for small lattice-constant differences,
the correction term can be approximated as $x\Delta
P_{\text{corr}}^{\text{(ZB ref)}} \simeq -2\epsilon_1(x) P_\text{f}^{m,
\text{ZB}}$ [see Eq.~(6) of the SM \cite{Supp}].  Therefore, the
difference between the implementations depends on the in-plane strain
and the difference between $P_\text{f}^{m, \text{ZB}}$ and
$P_{\text{eff}}^{n, \text{(H ref)}}$ [see Eq.~(7) of the SM
\cite{Supp}]. We see from Table \ref{sptab} that
$P_\text{eff}^{\text{ZB
ref}}=P_\text{f}^{\text{WZ}}-P_\text{f}^{\text{ZB}}$ is one to two
orders of magnitude less than $P_\text{eff}^{\text{(H
ref)}}=P_\text{f}^{\text{WZ}}$, therefore
$P_\text{eff}^{\text{(H ref)}}\sim P_\text{f}^{\text{ZB}}$.

The similarity between $P_\text{f}^{\text{ZB}}$ and
$P_\text{f}^{\text{WZ}}$ 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. 
%Therefore, considering a projection of the
%ion charges and bond/atom centered Wannier centers (which act rigorously as
%classical point charges, the movement of which gives the change in
%polarization \cite{KSV1993,VKS1993,Resta1994}) on the $c$ and [111]
%axes, the WZ and ZB structures would look similar except for small
%differences due to different lattice constants and $u$ parameters.


For AlGaN/GaN, the relatively modest difference in lattice constants
between AlN and GaN (and therefore modest $\varepsilon_1$ values for
coherently strained alloy layers), and an almost exact cancellation
between $P_\text{f}^{\text{GaN,ZB}}$ and
$P_{\text{eff}}^{\text{AlGaN,(H ref)}}$ means that the difference
between implementations is small over the whole composition range
(cf. Section~S4 of the SM \cite{Supp}). For InGaN/GaN, the large
lattice mismatch of InN and GaN and a less complete cancellation of
$P_\text{f}^{\text{GaN,ZB}}$ and $P_{\text{eff}}^{\text{InGaN,(H ref)}}$
results 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 \cite{Li2015}.

The two implementations differ significantly in the
\textit{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.~\ref{rel} the predicted
polarization bound charges for In$_{0.2}$Ga$_{0.4}$N/GaN (blue
curves) and Al$_{0.2}$Ga$_{0.4}$N/GaN (green curves) are shown
as a function
of strain relaxation of the layer, modeled by simply scaling
$\epsilon_1$. For both InGaN/GaN and AlGaN/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 influence the interface bound charge;
this effect has not been taken into account
in either version of the implementation.

\begin{figure}
\includegraphics[width=\columnwidth]{Fig3.eps}
\caption{\label{rel}
Absolute values for polarization sheet charges at the
In$_{0.2}$Ga$_{0.4}$N/GaN (blue) and Al$_{0.2}$Ga$_{0.4}$N/GaN
(green) interface as a function of percent strain relaxation.
0 \% relaxation corresponds to perfectly strained layers, 100 \% relaxation to
an unstrained overlayer at its bulk lattice constant.  Solid
curves are the correct implementation described in this work (H reference
for SP and improper PZ constants); dashed curves are the current
practice in the field (ZB reference for SP, without the correction
term, and proper PZ constants). }
\end{figure}

%
% For the case of InGaN/GaN, Fig.~\ref{Both}(a) shows that the
% predicted polarization charges are very similar for the two
% implementations.This in spite of the fact that, as seen from
% Table \ref{sptab}, there is a significant difference in the SP
% polarization between InN and GaN when a consistent reference
% structure is used (H, or with inclusion of a correction term for
% ZB), contrary to use of uncorrected ZB-referenced values. It is
% the use of improper PZ constants, along with the H reference,
% that brings the polarization values to a magnitude that is not
% very different from the values obtained using the current practice
% in the field (black dashed line).
%
%
% For InGaN/GaN heterostructures where the polarization is screened
% by free charge (analogous to the situation of AlGaN/GaN HEMTS),
% the proper PZ constants would be the correct ones to use, and
% the polarization bound charge would be significantly larger.
%
% For the case of AlGaN/GaN, it is critical to take into account the
% structure of the device on which the polarization measurements
% are performed. In general, the polarization will lie between
% that predicted using the proper and improper PZ constants (green
% dot-dashed and black dashed lines), depending on the extent
% to which the polarization field is screened. For the case of a
% HEMT-like structure, the results are expected to be closer to
% the predictions based on use of the proper PZ constants.  Again,
% these results (green dot-dashed line) are not very different
% from the current practice in the field (red line); here this
% can be attributed to the fact that the difference between the
% effective SP polarizations of AlN and GaN (Table~\ref{sptab})
% happens to result in similar predicted polarization sheet charges,
% regardless of whether the ZB or H reference structure is used. In
% the case of an AlGaN/GaN heterostructure with insulating interfaces
% and no free charge, however, the bound charge (calculated using
% improper PZ constants, black dashed line) is expected to be
% significantly lower.
%
% For the case of AlGaN/GaN, there is a large difference between the
% predicted sheet charges with the revised implementation of this
% work, and the current practice in the field. In this case we find
% that the differences in effective SP polarization between GaN and
% AlN are similar for the ZB and H reference structure (\textit{i.e.}
% the correction term is small), so it is the inclusion of the
% improper PZ constants which significantly decreases the predicted
% sheet charges in the revised implementation. The current practice
% in the field predicts higher sheet charges than reported in
% the literature, while the revised implementation predicts lower
% sheet charges.

\section{Conclusions}
\label{sec:conc}

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.

%
% Comparing the predicted sheet charges for AlGaN/GaN with reported
% experimental measurements suggests that the determination of 2DEG
% density may need to be reexamined.
%

\acknowledgements{

Work by C. E. D. was supported by the U.S. Department of Energy
(DOE), Office of Science, Basic Energy Sciences (BES) under Award
DE-SC0010689.  Work by A. J. was supported by the Center for Low
Energy Systems Technology (LEAST), one of the six SRC STARnet
Centers, sponsored by MARCO and DARPA.
Work by D. V. was supported by ONR Grant N00014-12-1-1035.
Computational resources were provided by the CSC at the CNSI and
MRL (an NSF MRSEC, DMR-1121053) (NSF CNS-0960316) and by the National
Energy Research Scientific Computing Center, a DOE Office of
Science User Facility supported by the Office of Science of the
U.S. Department of Energy under Contract No. DE-AC02-05CH11231.
}

%
% \bibliography{GaN}
%


\begin{thebibliography}{20}
\bibitem{Nye1957}
J. F. Nye (1957). {\it Physical Properties of Crystals}
(pp. 300). London: Oxford University Press.


\bibitem{Martin1972}
R. M. Martin, Phys. Rev. B {\bf 5}, 1607 (1972).


\bibitem{Feneberg2007}
M. Feneberg and K. Thonke, J. Phys.: Condens. Matter {\bf 19},
403201 (2007).


\bibitem{Vanderbilt2000}
D. Vanderbilt, J. Phys. Chem. Sol. {\bf 61}, 147 (2000).


\bibitem{BFV1997}
F. Bernardini, V. Fiorentini, D. Vanderbilt, Phys. Rev. B {\bf 56},
R10024 (1997).


\bibitem{BFV2001}
F. Bernardini, V. Fiorentini, D. Vanderbilt, Phys. Rev. B {\bf 63},
193201 (2001).


\bibitem{Bechstedt2000}
F. Bechstedt, U. Grossner, and J. Furthm\"{u}ller, Phys. Rev. B
{\bf 62}, 8003 (2000).


\bibitem{KSV1993}
R. D. King-Smith and D. Vanderbilt, Phys. Rev. B {\bf 47}, 1651
(1993).


\bibitem{Resta1994}
R. Resta, Rev. Mod. Phys. {\bf 66}, 899 (1994).


\bibitem{VKS1993}
D. Vanderbilt and R. D. King-Smith, Phys. Rev. B {\bf 48}, 4442 (1993).


\bibitem{Resta2007}
R. Resta and D. Vanderbilt, in \textit{Physics of Ferroelectrics:
A Modern Perspective}, edited by K. Rabe, Ch. H. Ahn, and
J.-M. Triscone (Springer-Verlag, Berlin, 2007).


\bibitem{silense}
SiLENSe Physics Summary v5.2 (2011).


\bibitem{Ambacher2000}
O. Ambacher, B. Foutz, J. Smart, J. R. Shealy, N. G. Weimann,
K. Chu, M. Murphy, R. Dimitrov, A. Mitchell, and M. Stutzmann,
J. Appl. Phys. {\bf 87}, 334 (2000).


\bibitem{Bennett2012}
J. W. Bennett, K. F. Garrity, K. M. Rabe, and D. Vanderbilt,
Phys. Rev. Lett. {\bf 109}, 167602 (2012).


\bibitem{HSE06}
J. Heyd, G. E. Scuseria, M. Ernzerhof, J. Chem. Phys. {\bf 118},
8207 (2003); {\it ibid.} {\bf 124}, 219906 (2006).


\bibitem{vasp1}
G. Kresse and J. Furthm$\ddot{\rm u}$ller, Phys. Rev. B. {\bf 54},
11169 (1996).


\bibitem{paw}
P. E. Bl\"{o}chl, Phys. Rev. B 50, 17953 (1994).


\bibitem{mp1976}
H.J. Monkhorst, J. D. Pack, Phys. Rev. B {\bf 13}, 5188 (1976).


\bibitem{Supp}
See Supplemental Material at [URL will be inserted by publisher] for a
table of structural properties and band gaps, a list of experimental
references, equations used to generate curves in Fig.~\ref{Both}, and
equations demonstrating the difference between implementations.


\bibitem{Ambacher1999}
O. Ambacher, R. Dimitrov, M. Stutzmann, B. E. Foutz, M. J. Murphy,
J. A. Smart, J. R. Shealy, N. G. Weimann, K. Chu, M. Chumbes,
B. Green, A. J. Sierakowski, W. J. Schaff, and L. F. Eastman,
Phys. status solidi (b) \textbf{216}, 381 (1999).


\bibitem{Vurgaftman2003}
I. Vurgaftman and J. R. Meyer, J. Appl. Phys. \textbf{94}, 3675 (2003).


\bibitem{Hangleiter2003}
A. Hangleiter, F. Hitzel, S. Lahmann, and U. Rossow,
Appl. Phys. Lett. \textbf{83}, 1169 (2003).


\bibitem{Renner2002}
F. Renner, P. Kiesel, G. H. D\"{o}hler, M. Kneissl, C. G. Van de
Walle, and N. M. Johnson, Appl. Phys. Lett. {\bf 81}, 490 (2002).


\bibitem{Turchinovich2003}
D. Turchinovich, P. Uhd Jepsen, B. S. Monozon, M. Koch, S. Lahmann,
U. Rossow, and A. Hangleiter, Phys. Rev. B {\bf 68}, 241307 (2003).


\bibitem{Cherns1999}
D. Cherns, J. Barnard, and F. A. Ponce, Solid State Commun. {\bf
111}, 281 (1999).


\bibitem{Zhang2004}
H. Zhang, E. J. Miller, E. T. Yu, C. Poblenz, and J. S. Speck,
J. Vac. Sci. Technol. B \textbf{22}, 2169 (2004).


\bibitem{Barker1973}
A. S. Barker, M. Ilegems, Phys. Rev. B {\bf 7}, 743 (1973).


\bibitem{Madelung}
O. Madelung, \textit{Semiconductors: Data Handbook} (Springer
Berlin Heidelberg, 2004).


\bibitem{Wu2012}
Y.-R. Wu, R. Shivaraman, K.-C. Wang, and J. S. Speck,
Appl. Phys. Lett. {\bf 101}, 083505 (2012).


\bibitem{Mcbride2014}
P. M. McBride, Q. Yan, and C. G. Van de Walle,
Appl. Phys. Lett. \textbf{105}, 083507 (2014).


\bibitem{Grandjean1999}
N. Grandjean, B. Damilano, S. Dalmasso, M. Leroux, M. La\:{u}gt,
and J. Massies, J. Appl. Phys. {\bf 86}, 3714 (1999).


\bibitem{Leroux1999}
M. Leroux, N. Grandjean, J. Massies, B. Gil, P. Lefebvre, and
P. Bigenwald, Phys. Rev. B \textbf{60}, 1496 (1999).


\bibitem{Yu1997}
E. T. Yu, G. J. Sullivan, P. M. Asbeck, C. D. Wang, D. Qiao,
and S. S. Lau, Appl. Phys. Lett. {\bf 71},2794 (1997).


\bibitem{Zhou2002}
Y. Zhou, B. Shen, T. Someya, H. Yu, J. Liu, H. Zhou, R. Zhang,
Y. Shi, Y. Zheng, and Y. Arakawa, Jpn. J. Appl. Phys. {\bf 41},
2531 (2002).


\bibitem{Arulkumaran2003}
S. Arulkumaran, T. Egawa, H. Ishikawa, and T. Jimbo, J. of
Vac. Sci. Technol. B {\bf 21}, 888 (2003).


\bibitem{Mkhoyan2004}
K. A. Mkhoyan, J. Silcox, Z. Yu, W. J. Schaff, L. F. Eastman,
J. Appl. Phys. {\bf 95}, 1843 (2004).


\bibitem{Li2015}
W. Li, S. Sharmin, H. Ilatikhameneh, R. Rahman, Y. Lu, J. Wang,
X. Yan, A. Seabaugh, G. Klimeck, D. Jena, and P. Fay, IEEE
J. Explor. Solid-State Comput. Devices Circuits \textbf{1}, 28
(2015).





% Surface/interfaces

%

% \bibitem{Yan2009} Q. Yan, P. Rinke, M. Scheffler, and C. G. Van

% de Walle, Appl. Phys. Lett. {\bf 95}, 121111 (2009).

%

% \bibitem{Takayama2001} T. Takayama, M. Yuri, K. Itoh, T. Baba,

% and J. S. Harris Jr, J. Cryst. Growth {\bf 222}, 29 (2001).

%

% \bibitem{Chen1992} W. Chen and T. G. Andersson,

% Semicond. Sci. Technol. \textbf{7}, 828 (1992).

%

% \bibitem{Zhang2000} Y. Zhang, I. P. Smorchkova, C. R. Elsass,

% S. Keller, J. P. Ibbetson, S. P. DenBaars, U. K. Mishra, J. Singh,

% J. Appl. Phys. {\bf 87}, 7981 (2000).

%

% \bibitem{BF2002} F. Bernardini and V. Fiorentini,

% Phys. Stat. Sol. (a) {\bf 190}, 65 (2002).

%

% \bibitem{Schulz1977} H. Schulz and K. H. Themann, Solid State

% Commun. {\bf 23}, 815 (1977).

\end{thebibliography}





% Unused bibitems
%
% \bibitem{Lax1973}
% D. Nelson and M. Lax, Phys. Rev. Lett. \textbf{31}, 763 (1973).
%
% \bibitem{Davydov1999}
% V. Yu. Davydov, V. V. Emtsev, I. N. Goncharuk, A. N. Smirnov,
% V. D. Petrikov, V. V. Mamutin, V. A. Vekshin, S. V. Ivanov,
% M. B. Smirnov, and T. Inus
% hima, Appl. Phys. Lett. {\bf 75}, 3297 (1999).
%
% \bibitem{Ibbetson2000}
% J. P. Ibbetson, P. T. Fini, K. D. Ness, S. P. DenBaars,
% J. S. Speck, and U. K. Mishra, Appl. Phys. Lett. {\bf 77}, 250
% (2000).
%
% \bibitem{Simon2001}
% J. Simon, R. Langer, A. Barski, M. Zervos, and N. T. Pelekanos,
% Phys. Stat. Sol. (a) {\bf 188}, 867 (2001).
%
% \bibitem{Kohler2011}
% K. K\"{o}hler, S. M\"{u}ller, P. Waltereit, W. Pletschen,
% V. Polyakov, T. Lim, L. Kirste, H. P. Menner, P. Br\"{u}ckner,
% O. Ambacher, C. Buchheim, a nd R. Goldhahn, J. Appl. Phys. {\bf
% 109}, 053705 (2011).
%
% \bibitem{Szabo1998}
% G. S\'{a}ghi-Szab\'{o}, R. Cohen, and H. Krakauer,
% Phys. Rev. Lett. \textbf{80}, 4321 (1998).
%
% \bibitem{Davidsson2004}
% S. K. Davidsson, M. Gurusinghe, T. G. Andersson, H. Zirath,
% J. Electron. Mater. {\bf 33}, 440 (2004).

\end{document}
