ORIGINAL PAPER
www.pss-b.com

The Electronic Properties of Chlorine in GaN: An Ab Initio
Study
Yusuke Fujii, Ruggero Micheletto, and Giovanni Alﬁeri*
Impurities can also be unintentionallyincorporated during processing steps likeetching. Dry etching is a fundamental stepfor the realization of semiconductor devices
and, in the case of GaN, it is typically favored
over wet etching for its selectivity.[13] In particular, Cl-based reactive ion etching (RIE)
and inductive coupled plasma (ICP) RIE
are chosen because they yield higher etchingrates than H-based RIE.[13]Despite Cl-based dry etching being
useful for mesa and termination etching,
this process was found to have an impact on the electronic properties of GaN. Deep levels in E gap , affecting the leakage current
and the breakdown voltage, were observed in both dry-etchedn- and p-GaN.[14–16] A decrease by two orders of magnitude ofcarrier concentration was found in Cl-based treated p-type
GaN.[17] In addition, controversial reports on how Cl-RIE affects
n-GaN ohmic contact resistivity (ρc ) are present in the literature:
while some authors have shown that Cl-RIE decreases ρc ,[18,19]
some other studies have shown that the presence of gallium
vacancies[20,21] and the formation of oxides on the surface[22,23]
increase ρc . Due to these effects, it is straightforward to ask
what the role of Cl, which is used in the RIE plasma, is on
the electronic properties or GaN.
To answer this question, we conducted a theoretical investigation, by means of density functional theory, of the electronic
properties of substitutional chlorine, both isolated and complexed with a vacancy or a dopant, in wurtzite GaN.

Chlorine-based reactive ion etching (RIE) is a fundamental processing step for the
manufacturing of GaN semiconductor devices. As impurities can be unintentionally incorporated in the crystal during processing, the electronic propertiesof chlorine in GaN are investigated. Density functional theory calculations of
substitutional Cl and related complexes (with a vacancy or a dopant) are carriedout. It is found that Cl and its complexes explain the reported effects of Cl
RIE-treated GaN on hole density and ohmic contact resistivity.

1. Introduction
Gallium nitride (GaN) has a wide bandgap (E gap ) and a high electron saturation velocity that make it the ideal candidate for highvoltage and high-frequency applications.[1] In addition, GaN has
also a direct E gap that allows the production of light-emitting
diodes (LEDs). Furthermore, even if homoepitaxial GaN is still
expensive, heteroepitaxial (Si, SiC, or Al2 O3 ) GaN is a valid alternative for the realization of cost-competitive power electronics.[1]Yet, the development of GaN optoelectronic devices cannot
leave out of consideration the study of electrically active point
defects. These are known to alter the functionality of devices;
therefore, it is important to study defects to be able to control
and/or limit their presence. Apart from intrinsic defects likethe nitrogen vacancy (V N ), which hinders the formation of p-type
GaN,[2–4] impurities can also be detrimental for applications.
This is the case of carbon and oxygen which are unintentionallyincorporated during growth. The electronic properties of carbon,
which sits at an N-site (CN), and its related complexes[5–7] havebeen thoroughly studied in the past. CN was found to be responsible for current collapse[8] and for the yellow luminescence(YL).[9] Oxygen (ON), which behaves as a shallow donor, is known
to decorate dislocations[10] and diffuse in the epilayer from a sapphire substrate.[11] In the latter case, it leads to p-type doping passivation due to the formation of neutral complexes with MgGa .[12]Y. Fujii, Prof. R. Micheletto
Department of Materials System Science
Yokohama City University
Yokohama 236-0027, Japan
Dr. G. Alﬁeri
ABB Power Grids Research
Segelhofstrasse 1A, 5405 Baden-Daettwil, Switzerland
E-mail: giovanni.alﬁeri@hitachi-powergrids.com
The ORCID identiﬁcation number(s) for the author(s) of this article
can be found under https://doi.org/10.1002/pssb.202000303.

DOI: 10.1002/pssb.202000303

Phys. Status Solidi B 2020, 2000303

2. Theoretical Methodology
The siesta code[24] was used for conducting density functional
theory calculations, using Troullier-Martins’ norm-conserving
pseudopotentials[25] to account for the effect of core electrons
and the Ceperley-Alder form of local density approximation(LDA) for the exchange correlation potential.[26] A double-ζplus polarization was used as atomic orbitals basis set throughout the calculations. For the Ga atoms the d-orbital electrons
were treated as valence electrons. The equilibrium conﬁgurations
pﬃﬃﬃ of the pristine and defected 96-atom GaN supercells(2 3a  3a  2c, a ¼ 3.184 Å, c ¼ 5.184 Å) were obtained byrelaxing the atomic coordinates with a conjugate gradient algorithm, until the maximum atomic forces were less than
0.025 eV Å1 and the stress tensor less than 0.5 GPa. The charge
density was projected onto a real space grid with an equivalent
cutoff of 375 Ry and eight Monkhorst-Pack special k-points[27]were used. The a ¼ b and c lattice constants were both 1% larger
than the experimental value with a direct bandgap of 2 eV.

2000303 (1 of 6)© 2020 Wiley-VCH GmbH

www.advancedsciencenews.com

www.pss-b.com

The formation energy (E form ) of a defect with charge state q
was calculated using the followng equationXqqni μi þ qðE F þ E V þ Δv0=b Þ þ E corr ðqÞ (1)
E form ¼ E tot  E bulk
tot 
i

q

E tot , ni , μi , E F , E V , Δv0=b , and E corr are the total energy (for thedefective or pristine supercell), the number of added (ni > 0) or
removed (ni < 0) i-atoms, the chemical potential of the i-atom,
the Fermi energy, the top of the valence band, the electrostatic
potential alignment term (between the bulk and defective supercells), and the ﬁnite-size correction term, respectively. This lastone (E corr ) was calculated following both the Makov–Payne
method[28] (hereafter MP) and the one reported by Freysoldt
et al.[29] (hereafter FNV), as implemented in a previous study.[30]As it is known that variations of E form (and thus of charge transition levels) depend on the choice of the functional, we applied acorrection scheme based on the alignment of electrons andchemical potential of the atomic species, as in a previous
study.[31] This correction is applied to the FNV results and will
be labeled as FNV*.
Calculation of μi was conducted as described in a previous
study.[32] μGa and μN were calculated using a Ga supercell
and N2 molecule, respectively. For the FNV* scheme, μi wascorrected using Equation (14) of a previous study[31] and using
the hybrid functional (HSE) value of the formation enthalpy(1.34 eV). For μCl , we used the LDA value. However, the optimized geometry of the Cl2 molecule, as obtained by siesta, wasalso used to calculate μCl with a B3LYP functional (kinetic energy
cutoff of 180 Ry).[33]3. Validation of the Methodology
To begin with, we conducted the calculation of the electronic
properties of two intrinsic defects, the V N and V Ga . This wasdone to evaluate the corrections to μN , μGa , and E V . Once these
corrections were obtained, we tested the FNV* scheme on two
defect complexes, one related to V N and another to V Ga .
Figure 1 shows the results of the LDA calculations carried out
on V N . Our LDA calculations predict the negative-U behavior of
V N .[3,4] The negative-U behavior (with U being the electron–
electron repulsion) occurs when the (q  1/q) charge state ofa defect lies closer to E V than the (q/q þ 1) one.[34,35] As a consequence, only a (q  1/q þ 1) charge state is found in E gap . The
MP correction (Figure 1a) does not show V N (þ/3þ) in the
Kohn-Sham E gap . As a matter of fact, this level is resonant in
the valence band (E V –0.30 eV) and so are the singly positive
and negative charge states, which are resonant in the conduction
band (shadowed in Figure 1a). When the FNV correction wasused, we obtained similar results to MP with V N (þ/3þ) predicted at E V  0.52 eV. Figure 1b shows the results of E form
calculatedusing the FNV* correction. By lowering E V relative
to VN(þ/3þ), we adjust the charge transition level position in
E gap . μN is corrected to match the E form of VN(þ/3þ) reported
by Diallo et al.[36] The corrected values of μN are 270.8 ≤ μN ≤
269.5 eV atom1 (from 271.2 ≤ μN ≤ 269.8 eV atom1 ). We
point out that, by using a bigger supercell (144 atoms), thevariation of E form is ≤ 0.1 eV.
Phys. Status Solidi B 2020, 2000303

Figure 1. Formation energy plot of the VN calculated with the a) MP andb) FNV* correction schemes. Red (black) line is for the N (Ga)-rich case.
The shadowed area represents the limit of the Kohn-Sham Egap.

We repeat analogous calculations for V Ga (Figure 2). The MP
correction (Figure 2a) does not show any donor level in the
Kohn-Sham E gap (E V  0.14 eV) and leads to the single, double,
and triple negative charge state transition levels at E V þ 0.33 eV,
E V þ 1.06 eV, and E V þ 2.08 eV, respectively. Using the FNV correction, the single donor level of VGa is located at E V  0.31 eV
and the single, double, and triple acceptors at E V þ 0.25 eV,
E V þ 0.88 eV, and E V þ 0.186 eV, respectively. By lowering
EV relative to VGa(0/þ),[36] we found that the acceptor levels
agree well with HSE results[3,36] to within 0.2 eV (Figure 2b).
The same occurs for E form (1808.5≤ μGa ≤1807.2 eV atom1from 1809.3≤ μGa ≤1807.8 eV atom1).[36]Having obtained the corrected values of μN , μGa , and E V , wethen proceed to validate the FNV* scheme used in this study.
To do this, we calculated the charge transition levels of the
MgGa V N and ONVGa complexes.

Figure 2. Formation energy plot of the VGa calculated with the a) MP andb) FNV* correction schemes. Red (black) line is for Ga (N)-rich case.

2000303 (2 of 6)© 2020 Wiley-VCH GmbH

www.advancedsciencenews.com

www.pss-b.com

The (0/2þ) charge transition level of MgGa V N [2] is foundresonant in the valence band at E V  0.19 and E V  0.39
for the MP and FNV correction scheme, respectively. To correct
the charge state transition level position in E gap , we employed
the value of EV used to adjust the VN(þ/3þ) transition level.
μMg is calculated using a Mg3 N2 supercell[2] and used to computethe E form of Mg(/0). μMg was corrected so that E form of Mg(/0)matched the HSE value reported in the literature.[2,3] With E V andμMg , we ﬁnd that MgGa V N (0/2þ) is located at EV þ 0.64 eV, 0.2 eVlower than the value reported in a previous study[2], and the E form
is in good agreement with a previous study.[2]The FNV (MP) correction scheme shows that ON V Ga (0/þ),
ON V Ga (/0), and ON V Ga (2/) are found at E V  0.30(20) eV,
EV þ 0.48(0.52) eV, and EV þ 1.32(1.68) eV, respectively. The
adjusted electron chemical potential for V Ga and the adjusted
value of μO [37] were used to calculate the charge state transitionlevel positions and their E form . The single donor, single, anddouble acceptors are found at EV þ 0.96 eV, EV þ 1.74 eV, andEV þ 2.5 eV, respectively. These levels are in agreement up to<0.20 eV with those obtained by HSEs, with the exception of
Table 1. Charge state transition levels of the VN, VGa, and ONVGa, calculatedusing the MP, FNV, and FNV* schemes. For ONVGa, the correction to EV is
the one taken from VGa. Numerical values indicate EV þ Et and units are in
eV. In brackets, the difference between the values obtained by HSE
calculations and those of the present study is shown.
Defect

MP

FNV

FNV*

V N (þ/3þ)

0.30

0.52

0.48 (0)3.13 (0.03)

V N (0/þ)

2.06

2.42

V Ga (0/þ)

0.14

0.31

0.97 (0)V Ga (/0)

0.32

0.25

1.51 (0.17)

V Ga (2/)

1.07

0.88

2.14 (0.16)

V Ga (3/2)

2.08

1.86

2.95 (0.15)

MgGa V N (0/2þ)

0.19

0.39

0.64 (0.20)

ON V Ga (0/þ)

0.20

0.30

0.96 (0.14)

ON V Ga (/0)

0.52

0.48

1.74 (0.06)

ON V Ga (2/)

1.68

1.32

2.58 (0.28)

the double donor level, which is 0.28 eV higher than the hybrid
counterpart.[37,38] Their E form also agrees reasonably well withthe literature.[37]A summary of the above results is shown in Table 1.

4. Results and Discussion
To ﬁnd out whether Cl sits at a N–(ClN) or at a Ga-site (ClGa), we
calculated the E form of ClN and ClGa. First, we focused on ClN, as
shown in Figure 3. After geometry optimization, the Cl─N bondlength ranges between 2.3 and 2.4 Å and ClN has a tetrahedral
symmetry. The ten electrons (three from VN and seven fromCl) occupy two singlets and one triplet state. Inspection of thehighest occupied molecular orbital (HOMO) reveals a fully occupied a1 state (Figure 3a), meaning that ClN is a double donor.
Analysis of E form , with the MP correction (Figure 3b), showsthat ClN is a double donor for each value of E F in the Kohn-Sham
E gap . Two transition levels (doubly and singly positive charge
states) are found resonant in the conduction band. The sameoccurs after application of the FNV correction.[39] Using theE V corrected for V N (þ/3þ), the E form of ClN shows only one
charge state transition level in E gap (Figure 3c), at EV þ 3.0 eV.
The use of μCl calculated using a HSE[33] results in a higher
E form of 0.4 eV. This applies also to the results shown in theremainder of this study.
Analogous calculations were also conducted on ClGa.
Geometry optimization results in a severe change of the Cl─N
bond length, with the shortest being 1.6 Å and the longest
2.6 Å. As a result, ClGa has a C s symmetry (Figure 4a), meaningthat the 12 electrons found in the defect molecule model areaccommodated in singlet levels. The analysis of the E form with
the MP correction reveals the presence of two charge states in
the Kohn-Sham E gap (Figure 4b). The singly and doubly negative
are found at EV þ 0.60 eV and EV þ 1.22 eV, respectively. To
conduct FNV* correction (Figure 4c), we use the value of EV
corrected for VGa(0/þ). This leads to the presence of four chargetransition levels in E gap at EVþ0.76 eV, EVþ1.56 eV, EVþ2.01 eV,
and EVþ2.70 eV for ClGa(þ/2þ), ClGa(0/þ), ClGa(/0), andClGa(2/), respectively.

Figure 3. a) Isosurfaces of the calculated wave function for ClN (0:09 e bohr3 ). The red (green) isosurface represents the positive (negative) valuesof the wave function. Formation energy plot of ClN calculated with the b) MP and c) FNV* correction schemes. Red (black) line is for N (Ga)-rich case.
Ga (N) atoms are yellow (blue), whereas Cl is light green.

Phys. Status Solidi B 2020, 2000303

2000303 (3 of 6)© 2020 Wiley-VCH GmbH

www.advancedsciencenews.com

www.pss-b.com

Figure 4. a) Isosurfaces of the calculated wave function for ClGa (0:03 e bohr3 ). The red (green) isosurface represents the positive (negative) valuesof the wave function. Ga (N) atoms are yellow (blue), whereas Cl is light green. Formation energy plot of ClGa calculated with the b) MP and c) FNV*
correction schemes. Red (black) line is for N (Ga)-rich case.

As shown in Figure 3c and 4c, the most stable conﬁguration,
for most of the E F -level values, is ClN. ClGa has a lower E form than
ClN (0.5 eV) only when 3.2 eV< E V þ E F ≤3.4 eV, that is, when
GaN is heavily n-type doped.
As Cl preferably occupies N-site for a wider range of E F values,
it is straightforward to analyze the ClN V Ga complex (Figure 5).
We considered the parallel conﬁguration, that is, the defect isparallel to the c-axis. This is because it was seen that for the
CN V Ga , no signiﬁcant difference of E form was present between
the parallel and perpendicular conﬁgurations.[6] Geometry optimization results in a C 3v symmetry, meaning that the sevenelectrons (ﬁve of V Ga and two of ClN ) can be found in two
singlets (a1 ) and one doublet (e). We ﬁnd that the spare electron
occupies one doublet as the HOMO is an e state (Figure 5a). As
Lyons et al.[38] have shown, the higher the number of donors
binding to V Ga , the lower the number of acceptors of the resulting complex. This means that if the analogy with the ON V Ga H
complex is correct,[38] the ClN V Ga complex will give rise to one
donor or one acceptor level (or both) in E gap .
Figure 5b shows the results of our calculations for the ClN V Ga
complex, with the MP scheme. We ﬁnd a donor level at
EV þ 0.21 eV and an acceptor level at EV þ 1.28 eV. Anotheracceptor level, doubly negative, is resonant in the conduction
band, at EV þ 2.8 eV. The FNV shows the singly and doubly
acceptor levels at EV þ 0.58 eV and EV þ 1.16 eV, respectively.
We also ﬁnd a donor level at EV þ 0.11 eV. To apply the FNV*
scheme, similarly to what we have done for the ON V Ga complex,
we consider the corrected EV of V Ga (0/þ). By doing so, we obtain
ClN V Ga (0/þ) at EV þ 1.61 eV, whereas the negative one is at
EV þ 2.70 eV (Figure 5c).
As the ClGa is more stable than ClN , in heavily doped n-GaN,
we also examined the ClGa V N complex. Similarly to ClN V Ga , thiscomplex was also studied in the parallel conﬁguration. Geometry
optimization, of the neutral and charged defects, shows that the
Cl atom leaves the Ga-site and moves into the N-site, that is.
ClGa V N transforms into ClN V Ga . An analogous transformation
was reported also in rare-earth-doped GaN[40] and demonstrates
that, although ClGa might be stable for a limited range of E F
values, it does not form complexes with native defects.
Next, we evaluate the effects of Cl on GaN dopant species.
To do this, we consider how ClN affects Mg, which is the p-type
dopant of choice for GaN and sits at a Ga-site.[41] The same cannot be done for ClGa , as the n-type dopant is Si and also sits at a
Ga-site. We investigated the ClN MgGa complex (Figure 6) in the

Figure 5. a) Isosurfaces of the calculated wave function for ClNVGa (0:09 e bohr3 ). The red (green) isosurface represents the positive (negative) valuesof the wave function. Red translucent sphere represents VGa. Formation energy plot of ClN calculated with the b) MP and c) FNV* correction schemes.
Red (black) line is for the N (Ga)-rich case.

Phys. Status Solidi B 2020, 2000303

2000303 (4 of 6)© 2020 Wiley-VCH GmbH

www.advancedsciencenews.com

www.pss-b.com

Figure 6. a) Isosurfaces of the calculated wave function for ClNMgGa (0:09 e bohr3 ). The red (green) isosurface represents the positive (negative)
values of the wave function. Mg atom is red. Formation energy plot of ClN calculated with the b) MP and c) FNV* correction schemes. Red (black)line is for the N (Ga)-rich case.

perpendicular conﬁguration (the Cl─Mg bond is perpendicular tothe c-axis), because the SiGa CN defect has a slighlty lower E form insuch a conﬁguration.[7] Based on the fact that ON is a single donor,
it was found that the ON MgGa complex is neutral.[12] This suggeststhat ClN MgGa could be a single donor. Geometry optimization
reveals that the ClN MgGa complex has a C 1 symmetry, meaningthat of the seven electrons of this complex (ﬁve from MgGa and two
from ClN ), there will be an unpaired electron in a singlet level(Figure 6a).
Using the MP correction scheme, we see that ClN MgGa is asingle donor throughout the Kohn-Sham E gap (Figure 6b). The
same occurs after applying the FNV correction. To apply the
FNV* scheme, we point out that as V Ga was used as a reference
for ClN VGa complex, in the present case we will use the MgGa
defect as a reference. We adjust E V and μMg to match the chargetransition level and E form of Mg(/0) calculated by HSE.[2,3] Then,these values are used for the ClN MgGa complex (Figure 6c) which
is found to be a single donor throughout the Kohn-Sham E gap .
A summary of the Cl-related charge state levels calculated inthis study is shown in Table 2.
As our results show that the electronic properties of Cl can havea signiﬁcant impact on GaN electronic devices. The reported
decrease by two orders of magnitude of the carrier concentrationin Cl-based-treated p-type GaN was attributed to the formation of
MgGa V N complexes.[17] However, this defect might already bepresent in the as-grown p-type material acting as a compensating
center. As our results show, ClN , ClN V Ga , and ClN MgGa behave as
donors in p-type GaN. Due to their low E form , ClN and ClN MgGa
might also be quite abundant, meaning that their presence cansigniﬁcantly reduce the hole density. ClN V Ga possesses thehighest E form possibly due to the fact that V Ga has a rather highE form in as-grown p-type.[3] However, high ion ﬂuxes andenergies, used during RIE,[13] might contribute to the formationof V Ga which, binding to Cl, can reduce the hole concentration,
possibly causing type conversion.[13]Regarding the controversy on the ohmic contact resistivity,
its improvement was attributed to the formation of donors, likeV N or divacancies (V N V Ga , V Ga ðV N Þ2 ).[18,19] However, in n-type
GaN, V N is neutral and divacancies are acceptors.[2,42] In addition, our results support the reported increase of ρc , following
Cl-RIE. As a matter of fact, if Cl is incorporated in GaN,
ClGa and ClN V Ga are acceptors in n-GaN and can lead to donor
compensation, thus worsening the ohmic behavior.
Last, we comment on the electron and chemical reservoir
method: as described in a previous study,[31] this is derived fromthe “marker method.” Typically, the charge state transition level
position in E gap of a well-localized defect A is corrected byshifting E V . This E V is then used to calculate the charge statetransition levels of another defect B. No constraints are givenon the geometry and/or symmetry of the defects. As the marker
method performs best when the wavefunctions of A and B
defects have the same symmetry,[43] the same should apply toa previous study.[31] The main obstacle to this approach is the
difﬁculty in ﬁnding suitable defects calculated by HSE, whose
symmetry matches the one under investigation. For instance,even if both ON V Ga [44] and ClN V Ga have C3v symmetry aftergeometry optimization, the respective highest occupied orbitals
do not, as it is a1 for ON V Ga [44] and e for ClN V Ga .

5. Conclusions
Table 2. Charge state transition levels detected in the present study
for the two substitutional (ClN and ClGa ) and complex ClN V Ga defects,
calculated using the FNV* correction scheme.
(2/)

(/0)

(0/þ)

(þ/2þ)

ClN

–

–

3.0

–

ClGa

0.76

1.56

2.01

2.70

ClN V Ga

–

1.61

2.70

–

Phys. Status Solidi B 2020, 2000303

Density functional theory calculations revealed that Cl sits at a
N-site and ClN behaves as a double donor, in agreement withthe fact that O is a single donor. ClN -related complexes, with
V Ga or MgGa , were also found to behave as donors. The presenceof these defects, due to possible unintentional incorporation
during processing, might explain the observed increase (decrease)of the electron (hole) concentration, reported in dry etched n-(p-)
type GaN.

2000303 (5 of 6)© 2020 Wiley-VCH GmbH

www.advancedsciencenews.com

www.pss-b.com

Conﬂict of Interest
The authors declare no conﬂict of interest.

Keywords
chlorine, density functional theory, dry etching, point defects
Received: June 5, 2020
Revised: August 14, 2020
Published online:[1] Y. Zhang, A. Dadgar, T. Palacios, J. Phys. D: Appl. Phys. 2018, 51,
273001.
[2] Q. Yan, A. Janotti, M. Schefﬂer, C. G. van de Walle, Appl. Phys. Lett.
2012, 100, 142110.
[3] G. Miceli, A. Pasquarello, Microelectron. Eng. 2015, 147, 51.
[4] J. L. Lyons, C. G. van de Walle, npj Comput. Mater. 2017, 3, 12.
[5] J. L. Lyons, A. Janotti, C. G. van de Walle, Appl. Phys. Lett. 2010, 97,
152108.
[6] M. Matsubara, E. Bellotti, Phys. Rev. B 2017, 121, 195701.
[7] M. Matsubara, E. Bellotti, Phys. Rev. B 2017, 121, 195702.
[8] M. J. Uren, M. Caesar, S. Karboyan, P. Moens, P. Vanmeerbeek,
M. Kuball, IEEE Electron Device Lett. 2015, 36, 826.
[9] S. O. Kucheyev, M. Toth, M. R. Phillips, J. S. Williams, C. Jagadish,G. Li, J. Appl. Phys. 2002, 91, 5867.
[10] P. P. Michalowski, S. Zlotnik, M. Rudzinski, Chem. Commun. 2019,55, 11539.
[11] R. Jakiela, E. Dumiszewska, P. Caban, A. Stonert, A. Turos, A. Barcz,
Phys. Status Solidi C 2011, 8, 1513.
[12] A. F. Wright, J. Appl. Phys. 2005, 98, 103531.
[13] S. J. Pearton, R. J. Shul, F. Ren, MRS Internet J. Nitride Semicond. Res.
2000, 5, 11.
[14] M. Kato, K. Mikamo, M. Ichimura, M. Kanechika, O. Ishiguro,
T. Kachi, J. Appl. Phys. 2008, 103, 093701.
[15] K. J. Choi, H. W. Jang, J. L. Lee, Appl. Phys. Lett. 2003, 82, 1233.
[16] C. Y. Lee, H. Sekiguchi, H. Okada, A. Wakahara, Jpn. J. Appl. Phys.
2012, 51, 076503.
[17] J. M. Lee, K. S. Lee, S. J. Park, J. Vac. Sci. Technol. 2004, 22, 479.
[18] A. T. Ping, Q. Chen, J. W. Yang, M. Asif Khan, I. Adesida, J. Electron.
Mater. 1998, 27, 261.
[19] T. Fujishima, S. Joglekar, D. Piedra, H. S. Lee, Y. Zhang, A. Uedono,
T. Palacios, Appl. Phys. Lett. 2013, 103, 083508.

Phys. Status Solidi B 2020, 2000303[20] S. J. Chua, H. W. Choi, J. Zhang, P. Li, Phys. Rev. B 2001, 64, 205302.
[21] Y. J. Han, S. Xue, W. P. Guo, C. Z. Sun, Z. B. Hao, Y. Luo, Jpn. J. Appl.
Phys. 2003, 42, 6409.
[22] A. A. Kobelev, Y. V. Barsukov, N. A. Andrianov, A. S. Smirnov,
J. Phys.: Conf. Ser. 2015, 586, 012013.
[23] L. K. Li, L. S. Tan, E. F. Chor, J. Cryst. Growth 2004, 268, 499.
[24] J. M. Soler, E. Artacho, J. D. Gale, A. Garcia, J. Junquera, P. Ordejon,
D. Sanchez-Portal, J. Phys.: Condens. Matter 2002, 14, 2745.
[25] N. Troullier, J. L. Martins, Phys. Rev. B 1991, 43, 1993.
[26] D. M. Ceperley, B. J. Alder, Phys. Rev. Lett. 1980, 45, 566.
[27] H. Monkhorst, J. D. Pack, Phys. Rev. B 1976, 13, 5188.
[28] G. Makov, M. C. Payne, Phys. Rev. B 1995, 51, 4014.
[29] C. Freysoldt, J. Neugebauer, C. G. van de Walle, Phys. Rev. Lett. 2009,
102, 016402.
[30] M. H. Naik, M. Jain, Comput. Phys. Commun. 2018, 226, 114.
[31] C. Freysoldt, B. Lange, J. Neugebauer, Q. Yan, J. L. Lyons, A. Janotti,
C. G. van de Walle, Phys. Rev. B 2016, 93, 165206.
[32] D. F. Hevia, C. Stampﬂ, F. Viñes, F. Illas, Phys. Rev. B 2013, 88,
085202.
[33] P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni,
D. Ceresoli, G. L. Chiarotti, M. Cococcioni, I. Dabo, A. Dal Corso,
S. Fabris, G. Fratesi, S. de Gironcoli, R. Gebauer, U. Gerstmann,
C. Gougoussis, A. Kokalj, M. Lazzeri, L. Martin-Samos, N. Marzari,
F. Mauri, R. Mazzarello, S. Paolini, A. Pasquarello, L. Paulatto,
C. Sbraccia, S. Scandolo, G. Sclauzero, A. P. Seitsonen, A. Smogunov,P. Umari, R. M. Wentzcovitch, J. Phys.: Condens. Matter 2009, 21,395502.
[34] G. D. Watkins, J. R. Troxell, Phys. Rev. Lett. 1980, 44, 593.
[35] J. Coutinho, V. P. Markevich, A. R. Peaker, J. Phys.: Condens. Matter
2020, 32, 323001.
[36] I. C. Diallo, D. O. Demchenko, Phys. Rev. Appl. 2016, 6, 064002.
[37] Z. Xie, Y. Sui, J. Buckeridge, A. A. Sokol, T. W. Keal, A. Walsh,Appl. Phys. Lett. 2018, 112, 262104.
[38] J. L. Lyons, A. Janotti, C. G. van de Walle, Phys. Status Solidi B 2015,
252, 900.
[39] C. Freysoldt, J. Neugebauer, C. G. van de Walle, Phys. Status Solidi B
2011, 248, 1067.
[40] S. Sanna, W. G. Schmidt, T. Frauenheim, U. Gerstman, Phys. Rev. B
2009, 80, 104120.
[41] G. Miceli, A. Pasquarello, Phys. Rev. B 2016, 93, 165207.
[42] H. Li, M. Huang, S. Chen, J. Semicond. 2020, 41, 032104.
[43] N. Fujita, Ph.D. dissertation, University of Exeter, 2009.
[44] A. Alkauskas, C. E. Dreyer, J. L. Lyons, C. G. Van de Walle, Phys. Rev. B
2016, 93, 201304(R).

2000303 (6 of 6)© 2020 Wiley-VCH GmbH