Electronic and optical properties of Y-doped Si3N4 by density functional theory

Electronic and optical properties of Y-doped Si3N4 by density functional theory

Journal of Alloys and Compounds 637 (2015) 376–381 Contents lists available at ScienceDirect Journal of Alloys and Compounds journal homepage: www.e...

3MB Sizes 2 Downloads 76 Views

Journal of Alloys and Compounds 637 (2015) 376–381

Contents lists available at ScienceDirect

Journal of Alloys and Compounds journal homepage: www.elsevier.com/locate/jalcom

Electronic and optical properties of Y-doped Si3N4 by density functional theory Zhifeng Huang a, Fei Chen a,b,⇑, Rui Su a, Zhihao Wang a, Junyang Li a, Qiang Shen a, Lianmeng Zhang a,b a b

State Key Lab of Advanced Technology for Materials Synthesis and Processing, Wuhan University of Technology, Wuhan 430070, China Key Laboratory of Advanced Technology for Specially Functional Materials, Ministry of Education, Wuhan University of Technology, Wuhan 430070, China

a r t i c l e

i n f o

Article history: Received 5 December 2014 Received in revised form 28 February 2015 Accepted 28 February 2015 Available online 5 March 2015 Keywords: Density functional theory Y-doped Si3N4 Electronic property Optical property

a b s t r a c t Geometry structures, formation energies, electronic and optical properties of Y-doped a-Si3N4 and bSi3N4 are investigated based on the density functional theory (DFT). The low values of formation energies indicate both Y-doped Si3N4 models can be easily synthesized. Besides, the negative formation energies of a-Yi-Si3N4 demonstrate that interstitial Y-doped a-Si3N4 has an excellent stability. The energies of impurity levels are different resulting from the different chemical environment around Y atoms. The impurity levels localized in the band gap reduces the maximum energy gaps, which enhances the optical properties of Si3N4. The static dielectric constants become larger and the optical absorption spectra show the red-shift phenomena for all Y-doped Si3N4 models. Ó 2015 Elsevier B.V. All rights reserved.

1. Introduction Silicon nitride (Si3N4) materials have been reported to show significant electronic and optical properties [1–3]. One of the most effective methods to regulate and control the electronic and optical property of Si3N4 nanomaterial is to change the energy level structure by introducing appropriate doping [4]. Xu et al. [5] prepared Eu2+-doped a-Si3N4 nanowires emitted strong yellow light ranging between 250 and 450 nm, which is related to the 4f6 5d  4f7 energy level transition of Eu2+. Gao et al. [6] found that the light absorption at 2.5 and 4.2 eV of heavily Al-doped Si3N4 nanobelt, which was ascribed to the two impurity bands introduced by Al doping. Wang et al. [7] synthesized a-Si3N4 nanowires showed that an intense violet-blue visible emission starting from 350 nm, which was due to the energy gap decreasing to 3.5 eV with heavily La doping. Simultaneously, first principle calculation based on the density functional theory (DFT) offers an effective method to research electronic structure and explain electron transition mechanism for nanoscale materials [8–10]. Ding et al. [11] found that c-Si3N4 showed the dielectric property at low Al substitution concentration but the metallic property at high Al concentration. Hart et al. [12] showed that the band gap of Ternary silicon germanium

⇑ Corresponding author at: State Key Lab of Advanced Technology for Materials Synthesis and Processing, Wuhan University of Technology, Wuhan 430070, China. Tel.: +86 27 87168606; fax: +86 27 87879468. E-mail address: [email protected] (F. Chen). http://dx.doi.org/10.1016/j.jallcom.2015.02.213 0925-8388/Ó 2015 Elsevier B.V. All rights reserved.

nitride (Si1xGex)3N4 was indirect for the silicon-rich compounds but became direct and decreased as the germanium content increases, due to greater mixing of s and p states in the conduction band. Yttrium has been extensively added into synthesizing Si3N4 based materials to stabilize crystal structures and improve optical properties [13–15]. Xu et al. [16] found that the calculated static dielectric constant increased and the energy gap decreased as Y content increased in c-Si3N4. The calculated density of states of interstitial Y-doped a-SiAlONs showed a new energy level appearing in the band gap originated from Y 4d orbit electrons [17,18]. However, researches have not been reported currently about the impact of local structure and bond characteristic on electronic and optical property of Y-doped Si3N4 based materials. Therefore, in this paper, Y-doped a-Si3N4 and b-Si3N4 models were discussed, which are marked as a-YSi-Si3N4, a-Yi-Si3N4 and b-YSi-Si3N4, respectively. Geometry optimizations and local deformations were firstly analyzed based on the density functional theory (DFT) calculation [19]. Secondly, formation energies were calculated to determine the stabilities of Y-doped models. Finally, combined with the results of local structures and bond characteristics, the electronic structures and optical properties caused by Y doping were discussed. 2. Methodology Fig. 1 shows the primitive cell of a-Si3N4 and b-Si3N4 models. Both a-Si3N4 and bSi3N4 crystallize in hexagonal structure [20]. For b-Si3N4, the periodicity is . . .ABAB. . . on the c axis. For a-Si3N4, the periodic stacking structure is ...ABCDABCD... on the c

377

Z. Huang et al. / Journal of Alloys and Compounds 637 (2015) 376–381

Fig. 1. Crystal structures of primitive (a) a-Si3N4 and (b) b-Si3N4 cell perpendicular to the c axis and (c) a-Si3N4 and (d) b-Si3N4 cell parallel to the c axis.

Fig. 2. Optimized geometry structures of Y-doped Si3N4 models: (a) a-YSi-Si3N4, (b) a-Yi-Si3N4 and (c) b-YSi-Si3N4.

Table 1 Optimized cell parameters and bond lengths for pure and Y-doped Si3N4 models. Supercells

a

c

Bond length (Å) b (Å)

c (Å)

a (o)

b (o)

c (o)

Volume (Å3)

Si–N

b-Si3N4 a-Si3N4 [17] a-Si3N4 [24] b-Si3N4 [20] b-Si3N4 [24] a-YSi-Si3N4

15.564 15.267 15.620 15.562 15.152 15.254 15.709

15.564 15.267 15.620 15.562 15.152 15.254 15.727

5.640 5.832 5.655 5.641 5.748 5.836 5.643

90.0 90.0

90.0 90.0

120.0 120.0

1.743a 1.740a 1.747a

90.0

90.0

120.0

90.0 90.2

90.0 90.2

120.0 120.1

1183.07 1177.24 1194.7 1182.80 1150.02 1176.13 1205.96

1.748a 1.673b

2.203 2.223 2.242 2.326

a-Yi-Si3N4

15.460

15.460

5.619

90.0

90.0

120.0

1163.01

1.741a 1.767b

2.388 2.573(3)c 2.576(3)c

b-YSi-Si3N4

15.339

15.302

5.863

90.0

90.0

119.4

1198.70

1.745a 1.664b

2.199(2)c 2.236 2.256

a-Si3N4

b

Cell parameters a (Å)

Y–N

The average value of all Si–N bonds. The N atoms of the Si–N bonds are adjacent to Y atom. The numbers in the parentheses show the number of bonds with same length.

axis. The CD layer is similar to AB layer except that it is rotated by 180° on the c axis. Therefore, the lattice parameter of a-Si3N4 in the c axis is approximately two times larger than that of b-Si3N4. Besides, there is a large interstice between AB layer and CD layer, which means that a-Si3N4 can accommodate a large interstitial atom. In our work, we considered 2  2  1 a-Si3N4 and 2  2  2 b-Si3N4 supercells, both of which contain 112 atoms and are of similar dimensions. The substitutional Y-doped a-Si3N4 and b-Si3N4 models are marked as a-YSi-Si3N4 and b-YSi-Si3N4, namely a Y atom substituting a Si atom in the supercells. The interstitial Y-doped a-Si3N4 is marked as a-Yi-Si3N4 that is a Y atom occupying the interstice of a-Si3N4 supercell.

All DFT calculations were performed using CASTEP package based on the plane wave pseudo potential approach [21]. The exchange–correlation function is Perdew Burke Ernzerhof (PEB) of the generalized gradient approximation (GGA) [22]. The interaction between valence electrons and the ionic core was described by the ultrasoft pseudo potential [23]. A plane-wave cutoff energy of 400 eV and a k-point mesh of 3  3  5 were used for the geometry optimizations. Reference configurations for the valence electrons are Si 3s 3p, N 2s 2p and Y 4p 4d 5s. The threshold for self-consistent field iterations is 5  105 eV per atom. The convergence tolerance parameters of optimized calculation were

378

Z. Huang et al. / Journal of Alloys and Compounds 637 (2015) 376–381

Table 2 Formation energies (eV) of Y-doped Si3N4 models. Models

N-rich

N-poor

a-YSi-Si3N4 a-Yi-Si3N4

5.83 25.94 5.87

6.66 28.42 6.70

b-YSi-Si3N4

the energy of 2.0  105 eV per atom, the maximum force of 0.05 eV Å, the maximum inner stress of 0.1 GPa and the maximum displacement of 2.0  104 nm. Higher plane-wave cutoff energy and more k-point meshes had also been used for the calculation in order to examine the accuracy of the calculations, the results showed almost no change in the total energy and geometry structure. After finishing geometry optimizations, formation energies, electronic and optical properties of Y-doped Si3N4 were calculated.

3. Results and discussion 3.1. Geometry optimization results Fig. 2 shows the optimized geometry structures of Y-doped Si3N4 supercells models. Table 1 summarizes the optimized cell parameters and bond lengths for all models. The calculated cell parameters of 2  2  1 a-Si3N4 and 2  2  2 b-Si3N4 supercells are in agreement with previous DFT results [17,20] and experimental values [24]. The minor changes mean the crystal structures of Ydoped a-Si3N4 and b-Si3N4 supercells are still hexagonal. The variations of cell parameters of Y-doped models are considered coming from the changes in bond lengths. Since the covalent radii of Y (180 Å) is larger than both that of Si (1 1 1 Å) and N (70 Å), local deformations of the Y-doped models are expected. As seen in Table 1, for a-YSi-Si3N4 and b-YSi-Si3N4, the different Y–N bonds mean that the Y atom are bonded to four adjacent N atoms forming an irregular tetrahedron, respectively. The average values of the Si–N bonds adjacent to Y atom are shorter by 4.0% and 4.4% for

Fig. 3. Band structures and projected density of states plots for (a) a-Si3N4, (b) aYSi-Si3N4 and (c) a-Yi-Si3N4.

Fig. 4. Band structures and projected density of states plots for (a) b-Si3N4 and (b) b-YSi-Si3N4.

Z. Huang et al. / Journal of Alloys and Compounds 637 (2015) 376–381

379

Fig. 5. Partial optimized geometry structures and corresponding electron density difference plots: Si–N–Si–N in (a) a-Si3N4 and (b) b-Si3N4, Y–N–Si–N in (c) a-YSi-Si3N4, (d) b-YSi-Si3N4 and (e) a-Yi-Si3N4.

3.2. Formation energies

Table 3 Mulliken populations of pure and Y-doped Si3N4 models. Supercells

Y population

N population

Si population

Adjacent to Y

Adjacent to Y

Si

N

a-YSi-Si3N4

1.84 1.85 1.82

1.38 1.39 1.37

1.83

1.28(2)a 1.21 1.18

a-Yi-Si3N4

1.81

1.37

1.43

1.37 1.34(3)a 1.31(3)a

1.77(3)a 1.68(3)a 1.63(6)a

b-YSi-Si3N4

1.82

1.37

1.79

1.30 1.27 1.23(2)a

1.74 1.64

a-Si3N4 b-Si3N4

a

Average population

The numbers in the parentheses show the same Mulliken populations.

a-YSi-Si3N4 and b-YSi-Si3N4, respectively. For a-Yi-Si3N4, there are eight N atoms around the Y atom, however, the axial distance (5.619 Å) of the top and bottom two N atoms is too large to form two equivalent Y–N bonds. Therefore, the Y atom is off-centered of the eight N atoms and forms only one strong axial Y–N bond (2.388 Å), considerably shorter than the six lateral Y–N bonds (>2.57 Å), which is consistent with the previous result [17]. The strong forces between Y and N atoms lead to the regional narrowing. These small changes of Si–N bonds near to the Y atom eliminate the lattice destruction of Si3N4 brought by introducing the big Y atom.

To investigate the stabilities of Y-doped Si3N4 models, the formation energies were calculated by the following formula [25]:

Eform ¼ Edoped  Epure  mlY þ nlSi

ð1Þ

where Eform and Epure are the total energy of Y-doped and pure Si3N4. lY, lSi and lN are the chemical potentials of Y, Si and N atoms, respectively. The coefficient m is the number of the Y atom in Y-doped model and n is the number of the Si atom substituted by the Y atom. It is important to note that the formation energy is not fixed but depends on growing conditions, which can be changed from N-rich to N-poor (Si-rich) [26]. The relations of lN, lY and lSi as follows:

lSi þ 4lN ¼ laSi3 N4

ð2Þ

lY þ lN ¼ lYN

ð3Þ

Under N-rich growing condition, lN is gained from the groundstate energy of the N2 molecule (lN = lN/2), while lY and lSi are calculated by the formulae (2) and (3), respectively. Under N-poor growing condition, lSi is determined from the energy of one Si atom in bulk Si (lSi = lbulk-Si), then lN and lY obtained by the formulae (2) and (3), respectively. Table 2 exhibits formation energies of Y-doped Si3N4 models under N-rich and N-poor growing conditions. The small values of formation energies indicate that Y-doped Si3N4 models can be easily synthesized. The negative values of formation energies for a-Yi-Si3N4 demonstrate that interstitial

380

Z. Huang et al. / Journal of Alloys and Compounds 637 (2015) 376–381

Y-doped a-Si3N4 has an excellent stability, which is the reason yttrium has been extensively used as stabilizer in synthesizing Si3N4 based materials [13–15]. 3.3. Electronic properties Figs. 3 and 4 show the calculated band structures and projected density of states (PDOS) plots for pure and Y-doped Si3N4 models. The dashed lines in the band structures represent the Fermi level (EF). The energy gaps (Eg) are 4.78 and 4.31 eV for pure a-Si3N4 and b-Si3N4, respectively, which are consistent with the previous DFT result [27]. The bottom of conduction bands is mainly constituted of Si 3s 3p states, and the top of valence bands is occupied by N 2p states for both pure a-Si3N4 and b-Si3N4. For a-YSi-Si3N4 and b-YSi-Si3N4, the bottom of conduction bands is consists of Si 3s 3p and Y 4p 4d states, while the top of the valence bands is N 2p states. Besides, for a-YSi-Si3N4, the conduction band extends forward the valence band since the energy of Y 4p states are lower than Si 3s 3p states, resulting to the Eg decreased to 4.09 eV. For b-YSi-Si3N4, the hybrid energy levels of Y 4d and Si 3p states below the conduction band reduce the Eg to 3.65 eV. For a-Yi-Si3N4, Si 3s 3p and Y 4p 4d 5s states occupy the bottom of conduction band and N 2p states occupy the top of valence band. Besides, Y 4d, Si 3p and N 2p states appear at the same energy region below the conduction band, which indicates the hybridization occurs between Y 4d with Si 3p and N 2p orbits. As a result, the Eg is 4.51 eV but the maximum energy gap is reduced to 2.47 eV. To investigate the effects of bond characteristics on electronic properties, electron density difference plots and Mulliken populations are discussed. Fig. 5 shows the partial optimized geometry structures and corresponding electron density difference plots of all models. The electron density of Si–N bond gathers toward from the Si atom to the N atom obviously. However, in Y-doped Si3N4 models, the electron density of Y–N bond extend not visibly from the Y atom to the N atom, which indicates the Y–N bond strengths are weaker than Si–N bond. Table 3 shows the Mulliken populations of all models. The average populations of N and Si atoms in pure a-Si3N4 and b-Si3N4 are different to their formal charges 3 and +4, indicating the ioniccovalent character for Si–N bond. The changes of populations of N and Si atoms in Y-doped models suggested the changes of electron transfers. The populations of N atoms adjacent to Y atom for a-YSi-Si3N4 and b-YSi-Si3N4 are larger than that of N atoms bonded to Si atoms, which indicates less electrons transferring from Y atom to adjacent N atoms compared with Si atoms to N atoms. Consequently, the quantity of vacancy bands in N 2p states increased and the top of valence bands moves up to the Fermi level, which is corresponding with the PDOS showed in Figs. 3(b) and 4(b). For a-Yi-Si3N4, the populations of Si atoms adjacent to Y atom are smaller than 1.84, which means less electrons transferring from Si atoms to N atoms due to electrons transferring from Y atom to N atoms. The population 1.37 of just one N atom indicates only one strong axial Y–N bond between the Y atom and axial N atom is formed, and the lateral Y–N bonds are relatively weaker in a-YiSi3N4. In a word, the different bond characteristics and electron transfers are originated from the different chemical environments of Y atoms, which generates the difference of the impurity levels positions for Y-doped Si3N4 models.

12 pffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi e21 ðxÞ þ e22 ðxÞ  e1 ðxÞ

gðxÞ ¼ 2x

ð4Þ

Fig. 6 illustrates the real parts e1(x) and imaginary parts e2(x) of the dielectric function e(x) of all models. The calculated e1(x) and e2(x) of pure a-Si3N4 and b-Si3N4 are similar to the results obtained with the ab initio wave functions calculation by Xu and Ching [29]. The main peaks at about 7.0 and 10.0 eV of e2(x) are originated from the electron transitions between the occupied states of N 2p states in the top of valence bands and empty states in the conduction bands for both a-Si3N4 and b-Si3N4. For a-YSi-Si3N4, The new peak at 0.2 eV originates from the electron transition from the occupied states to the empty states in N 2p orbitals at the top of valence band. For a-Yi-Si3N4, the new peak starting at 2.9 eV is derived from the electron transition from the occupied states of N 2p states in the top of valence band to hybrid energy levels of Y 4d with Si 3p and N 2p. For b-YSi-Si3N4, the new strong peak at about 0–2.0 eV stems from the occupied states to the empty states in N 2p orbitals at the top of valence band and the hybrid energy levels of Y 4d and Si 3p states to the conduction band. The zero frequency of e1(x) reflects the static dielectric constant e1(0). The static dielectric constants are 5.23, 4.23 and 13.56 for a-YSi-Si3N4, a-Yi-Si3N4 and b-YSi-Si3N4, which are larger than 4.19 and 4.05 for pure a-Si3N4 and b-Si3N4. The increases of the static dielectric constants are caused by the reductions of the Eg. The large static dielectric constant is corresponding to the large refractive index, which means that Y-doped Si3N4 materials may have special application in antireflection coating for solar cells [30].

3.4. Optical properties Evaluation of optical properties is based on the dielectric function, eðxÞ ¼ e1 ðxÞ þ ie2 ðxÞ, where e1(x) and e2(x) are the real parts and imaginary part, respectively [28]. The absorption coefficient g(x) can be obtained directly from the following formula:

Fig. 6. Dielectric function plots for pure and Y-doped (a) a-Si3N4 and (b) b-Si3N4 models.

Z. Huang et al. / Journal of Alloys and Compounds 637 (2015) 376–381

381

can be easily synthesized and the negative formation energies demonstrate the stability of a-Si3N4 is improved by interstitial Y doping. The different bond characteristics and electron transfers are originated from the different chemical environments of Y atoms, which generates the difference of the impurity levels positions for Y-doped models. The impurity levels localized in the band gap reduces the maximum energy gaps for all Y-doped models, which enhances the optical properties of Si3N4. With the decreases of the Eg, the static dielectric constants for all Y-doped Si3N4 models become larger. The optical absorption spectra show the red-shift phenomena for all Y-doped Si3N4 models, especially for a-Yi-Si3N4, of which the absorption edge energy of 2.9 eV is corresponding to purple light. This study will provide the theoretical foundation for Si3N4 materials used as novel optical materials. Acknowledgement The project is supported by the National Natural Science Foundation of China (Nos. 51202171 and 51472188), the Specialized Research Fund for the Doctoral Program of Higher Education of China (No. 20120143120004) and the ‘‘111’’ project (No. B13035). References

Fig. 7. Absorption spectra plots for pure and Y-doped (a) a-Si3N4 and (b) b-Si3N4 models.

Fig. 7 shows the absorption spectra of all models. The strong absorption bands of all models locate at about 4.5–15.5 eV corresponding to electron transitions in N 2p states to the empty states in the conduction band. Compared with pure a-Si3N4 and b-Si3N4, the absorption energies of starting positions and host peaks for Y-doped models become smaller, resulting from the maximum energy gap reductions. The lower absorption energies signify the red-shift phenomenon occur after Y doping, and the magnitudes of the red-shift are shown as the light red rectangles in Fig. 7. The energies of the absorption host peaks are almost no difference. However, the starting positions of absorption peaks are red-shift by 0.7, 1.9 and 0.6 eV for a-YSi-Si3N4, a-Yi-Si3N4 and b-YSi-Si3N4, respectively. The absorption edge energy of 2.9 eV is corresponding to the wavelength of 427 nm, which means a-YiSi3N4 has the potential as purple light absorber. 4. Conclusions In this paper, geometry structures, formation energies, electronic and optical properties of Y-doped a-Si3N4 and b-Si3N4 are systematically investigated based on DFT. The low values of formation energies indicate that both the Y-doped a-Si3N4 and b-Si3N4

[1] M. Zeuner, S. Pagano, W. Schnick, Angew. Chem. Int. Ed. 50 (2011) 7754–7775. [2] L.G. Zhang, H. Jin, W.Y. Yang, Z.P. Xie, H.Z. Miao, L.N. An, Appl. Phys. Lett. 86 (2005) 061908. [3] L.W. Yin, Y. Bando, Y.C. Zhu, Y.B. Li, Appl. Phys. Lett. 83 (2003) 3584–3586. [4] H. Wang, Y. Chen, Y. Kaneta, S. Iwata, J. Alloys Comp. 491 (2010) 550–559. [5] X. Xu, T. Nishimura, Q. Huang, R.J. Xie, N. Hirosaki, H. Tanaka, J. Am. Ceram. Soc. 90 (2007) 4047–4049. [6] F. Gao, Y. Wang, L. Zhang, W. Yang, L. An, J. Am. Ceram. Soc. 93 (2010) 1364– 1367. [7] Z. Wang, Z. Huang, F. Chen, Q. Shen, L. Zhang, J. Lumin. 151 (2014) 66–70. [8] H.T. Hintzen, M.R.M.M. Hendrix, H. Wondergem, C.M. Fang, T. Sekine, G. de With, J. Alloys Comp. 351 (2003) 40–42. [9] C.D. Valentin, G. Palma, G. Pacchioni, J. Phys. Chem. C 115 (2011) 561–569. [10] A. Haddou, H. Khachai, R. Khenata, F. Litimein, A. Bouhemadou, G. Murtaza, Z.A. Alahmed, S. Bin-Omran, B. Abbar, J. Mater. Sci. 48 (2013) 8235–8243. [11] Y.C. Ding, A.P. Xiang, J. Luo, X.J. He, Q. Cai, X.F. Hu, Physica B 405 (2010) 828– 833. [12] J.N. Hart, N.L. Allan, F. Claeyssens, Phys. Rev. B 84 (2011) 245209. [13] B. Joshi, G. Gyawali, H. Wang, T. Sekino, S.W. Lee, J. Alloys Comp. 557 (2013) 112–119. [14] C.R. Zhou, Z.B. Yu, V.D. Krstic, J. Eur. Ceram. Soc. 27 (2007) 437–443. [15] L. Ceja-Cárdenas, J. Lemus-Ruíz, D. Jaramillo-Vigueras, S.D. De la Torre, J. Alloys Comp. 501 (2010) 345–351. [16] M. Xu, Y.C. Ding, G. Xiong, W.J. Zhu, H.L. He, Physica B 403 (2008) 2515–2520. [17] L. Benco, J. Hafner, Z. Lences, P. Sajgalik, J. Eur. Ceram. Soc. 28 (2008) 995– 1002. [18] T. Suehiro, H. Onuma, N. Hirosaki, R.J. Xie, T. Sato, A. Miyamoto, J. Phys. Chem. C 114 (2009) 1337–1342. [19] W. Kohn, L.J. Sham, Phys. Rev. A 140 (1965) 1133–1138. [20] A. Kuwabara, K. Matsunaga, I. Tanaka, Phys. Rev. B 78 (2008) 064104. [21] M.D. Segall, P.J.D. Lindan, M.J. Probert, C.J. Pickard, P.J. Hasnip, S.J. Clark, M.C. Payne, J. Phys.-Condens. Matter 14 (2002) 2717–2744. [22] J.P. Perdew, K. Burke, M. Ernzerhof, Phys. Rev. Lett. 77 (1996) 3865–3868. [23] D. Vanderbilt, Phys. Rev. B 41 (1990) 7892–7895. [24] M. Billy, J.C. Labbe, A. Selvaraj, G. Roult, Mater. Res. Bull. 18 (1983) 921–934. [25] R. Long, N.J. English, J. Phys. Chem. C 113 (2009) 8373–8377. [26] F. Oba, K. Tatsumi, I. Tanaka, H. Adachi, J. Am. Ceram. Soc. 85 (2002) 97–100. [27] K. Ulman, R. Sathiyanarayanan, R.K. Pandey, K.V.R.M. Murali, S. Narasimhan, J. Appl. Phys. 113 (2013) 234102. [28] H. Ehrenreich, H.R. Philipp, Phys. Rev. 128 (1962) 1622–1629. [29] Y.N. Xu, W.Y. Ching, Phys. Rev. B 51 (1995) 17379–17389. [30] R.J. Xie, H.T. Hintzen, J. Am. Ceram. Soc. 96 (2013) 665–687.