Piezoelectricity in K1− xNa xNbO3: First-principles calculation*
Li Qiang, Zhang Rui†, Lv Tian-Quan‡, Zheng Li-Mei
Condensed Matter Science and Technology Institute, Harbin Institute of Technology, Harbin 150080, China

Corresponding author. E-mail: ruizhang ccmst@hit.edu.cn

Corresponding author. E-mail: ltq@hit.edu.cn

*Project supported by the National Basic Research Program of China (Grant No. 2013CB632900).

Abstract

The piezoelectric properties of K1− xNa xNbO3 are studied by using first-principles calculations within virtual crystal approximation. To understand the critical factors for the high piezoelectric response in K1− xNa xNbO3, the total energy, piezoelectric coefficient, elastic property, density of state, Born effective charge, and energy barrier on polarization rotation paths are systematically investigated. The morphotropic phase boundary in K1− xNa xNbO3 is predicted to occur at x = 0.521, which is in good agreement with the available experimental data. At the morphotropic phase boundary, the longitudinal piezoelectric coefficient d33 of orthorhombic K0.5Na0.5NbO3 reaches a maximum value. The rotated maximum of is found to be along the 50° direction away from the spontaneous polarization (close to the [001] direction). The moderate bulk and shear modulus are conducive to improving the piezoelectric response. By analyzing the energy barrier on polarization rotation paths, it is found that the polarization rotation of orthorhombic K0.5Na0.5NbO3 becomes easier compared with orthorhombic KNbO3, which proves that the high piezoelectric response is attributed to the flattening of the free energy at compositions close to the morphotropic phase boundary.

Keyword: 31.15.A–; 77.65.–j; 77.84.Cg; first principles; morphotropic phase boundary; piezoelectricity; virtual crystal approximation
1. Introduction

Lead-based piezoelectric materials have been widely used in sensors, actuators, and accelerometers.[1] However, due to Pb’ s toxicity, there is an increasing attention in the development of high-performance lead-free piezoelectric materials. Environmentally friendly (K, Na)NbO3 (KNN) and its solid solutions show the superior piezoelectric properties at a K/Na ratio of around 50/50, and are the potential substitutes for lead-based piezoelectrics.[27] It is well known that the piezoelectric response is enhanced for piezoelectrics at the morphotropic phase boundary (MPB).[1, 8] Therefore, a great deal of experimental effort has been devoted to search the MPB composition in the K1− xNaxNbO3 system. Studies[912] on the phase transition process in the KNN system suggested that an MPB lies near the composition with 50  mol.% NaNbO3. In those studies, the Na content intervals are too large to verify the accurate K/Na ratio at the MPB. The small intervals in K1− xNaxNbO3 ceramics have been derived by Dai et al.[13] Their results suggested that a typical morphotropic phase boundary exists at x = 0.52– 0.525. Although the accurate MPB in K1− xNaxNbO3 has been confirmed experimentally, some fundamental problems still need to be solved, such as the dominant factor for the high piezoelectric response in K1− xNaxNbO3, the effect of K– O (or Na– O) and Nb– O bonding on the piezoelectric performance, and the effect of Na substitution on the polarization rotation. First-principles calculations can predict various properties well, [1418] such as the phase transition and thermodynamic properties, which have been adopted for the purpose of this report.

In this paper, we use first-principles density functional theory within the virtual crystal approximation to predict the MPB in K1− xNaxNbO3 and evaluate its piezoelectric response behaviors. Moreover, the mechanism of a high piezoelectric response and the effect of the domain-engineering method on polarization rotation in K1− xNaxNbO3 near the MPB are also investigated. To the best of our best knowledge, it is the first time that the MPB of KNN has been obtained theoretically. This paper is organized as follows. Theoretical methods are described in Section 2. The results for the structural parameter, bulk modulus, band gap, Born dynamical charge, and piezoelectric property are discussed in Section 3. Finally, the conclusions are given in Section 4.

2. Computational methods

The first-principles calculations are conducted using the density functional theory (DFT) and plane wave pseudopotential method implemented in the ABINIT software package.[1922] In order to quantitatively investigate the MPB and related phenomena, the virtual crystal approximation (VCA)[23] is adopted because it has been successfully applied to some perovskite solutions (e.g., PZT, [24, 25] NBT-BT, [26] PSN, [27] and LNTO[28]). The exchange correlation potential is described by the local-density approximation within the Hartwigsen– Goedecker– Hutter (HGH) pseudopotentials.[29] We treat 9 valence electrons for K (3s23p64s1) and Na (2s22p63s1), 6 for O (2s22p4), and 13 for Nb (4s24p64d45s1). The Brillouin zone (BZ) is sampled by a 5 × 5 × 5 Monkhorst– Pack mesh of k points for cubic, rhombohedra, orthorhombic, and tetragonal K1− xNaxNbO3. To achieve high accuracy in the total energy calculation and ionic relaxation, we use a very high plane-wave energy cutoff of 3591 eV, and the convergence for the maximum tolerance on the wave function squared residual is set to 10− 18. Firstly, the appropriate cell parameters of K1− xNaxNbO3 are obtained after optimizing the lattice constants and internal atom coordinates based on the maximal force tolerance (0.0257  eV/Å ). Secondly, the total energy, piezoelectric constant, Born effective charge, and elastic property calculations for K1− xNaxNbO3 are performed by the density functional perturbation theory (DFPT).[30]

3. Results and discussion
3.1. Structure properties

Above the Curie temperature (around 624  K), KNbO3 (KNO) is a typical perovskite structure with the cubic symmetry. Cooling KNO from the Curie temperature leads to a continuous phase transition in a cubic-tetragonal-orthorhombic-rhombohedral sequence. The calculated lattice parameter, bulk modulus, and band gap for these phases are presented in Table  1. The lattice constant a and the cell volume V of the cubic phase are calculated to be 3.968  Å and 62.46  Å 3, respectively, which are in good agreement with the experimental values (4.021  Å and 65.01  Å 3).[32] The calculated c/a ratio for the tetragonal KNO is 1.034, which agrees with the available experimental result (1.017).[36] According to the present DFT calculations, the unit-cell parameters of a = 5.629  Å , b = 3.959  Å , and c = 5.633  Å are predicted for the orthorhombic KNO. The rhombohedral KNO with the space group R3m only exists below room temperature.

Table 1. Lattice parameter, volume, bulk modulus, band gap, and piezoelectric constant of the cubic, tetragonal, orthorhombic, and rhombohedral KNbO3, compared with experimental and theoretical data. For each phase, the error of lattice parameter with respect to the experimental data is given.

Compared with the reported experimental data, [32, 36, 37] our calculated cell parameters listed in Table  1 have an accuracy with a relative error of ± 1.6%. The predicted bulk modulus for the cubic and orthorhombic KNO are 132  GPa and 185  GPa, respectively, which are consistent with the experimental data (138  GPa and 172  GPa).[34, 38] The calculated piezoelectric coefficients d33 are also in reasonable agreement with the experimental and theoretical data.[33, 39]

3.3. Piezoelectric properties

At the MPB, the piezoelectric materials would show the high piezoelectric properties. In order to further verify our prediction, we calculate the piezoelectric coefficients (d33 and d15). Figure  2 only shows d33 of the orthorhombic K1− xNaxNbO3 in the range of x = 0.0– 0.7, because the elastic constant calculations suggest that the orthorhombic K1− xNaxNbO3 becomes mechanically instable in the range of x > 0.7. The d33 of the orthorhombic structure shows a maximum value at x = 0.5. This result is in agreement with the experimental measurements.[912] The reason for the peak value at x = 0.2 is unclear yet. There may be another MPB or phase boundary. Similarly, the d33 also shows a maximum at x = 0.2 in both BiFe1− xCoxO and (NBT)1− x-(K1/2BiTiO3)x, which has been thoroughly verified by the experimental[4143] and theoretical[44] studies.

Fig.  2. d33 of the orthorhombic K1− xNaxNbO3 as a function of x.

The shear piezoelectric coefficient d15 is also very important for piezoelectrics. Therefore, we calculate d15 as a function of x for the orthorhombic K1− xNaxNbO3. As shown in Fig.  3, there are three peaks located at x = 0.2, 0.464, and 0.7. Although d15 at x = 0.7 is the highest among the three peaks, the corresponding d33 is the lowest. At x = 0.464, d15 is calculated to be 153.52 pC/N, which further identifies the MPB in K1− xNaxNbO3.

Fig.  3. d15 of the orthorhombic K1− xNaxNbO3 as a function of x.

3.4. Orientation dependence of piezoelectric coefficient

According to the suggestion of Davis et al., [45] the largest longitudinal piezoelectric response would be found away from the polar direction when the ratio d15/d33 is larger than 1.5, and the crystals can be categorized into “ rotator” ferroelectrics with d15/d33 > 1.5 and “ extender” with d15/d33 < 1.5. According to the predicted piezoelectric coefficients d33 and d15, the d15/d33 value is 7.29 for the orthorhombic K0.5Na0.5NbO3. Thus, the orthorhombic K0.5Na0.5NbO3 belongs to “ rotator” ferroelectrics and is highly anisotropic. The large anisotropy is also observed in other single-domain relaxor-PT-based crystals near their MPBs.[8] It is necessary for investigating the crystal orientation dependence of the piezoelectric properties. From the calculated single-domain data of the orthorhombic K1− xNaxNbO3 with x = 0.5 and 0.0, the orientation dependence of the intrinsic piezoelectric coefficient is calculated through the transformation of the single domain piezoelectric matrix, [8, 46, 47] and the results are presented in Fig.  4. One can see that the highest value of the piezoelectric coefficient of the orthorhombic K1− xNaxNbO3 with x = 0.5 is on the order of 58.4  pC/N, being along the θ max = 50° away from the spontaneous polarization [011] direction (close to the [001] direction), larger than that withx = 0.0. Thus, the piezoelectric coefficients exhibit strong orientation dependence.

Fig.  4. Orientation dependence of piezoelectric coefficient of the orthorhombic K1− xNaxNbO3 with x = 0.5 and x = 0.0.
3.5. Bulk modulus and shear modulus

For the orthorhombic K1− xNaxNbO3 with x = 0.5, there are nine independent elastic constants (c11, c22, c33, c44, c55, c66, c12, c13, and c23), which are obtained using DFPT. From the calculated elastic constants, both the bulk modulus (B) and shear modulus (G) can be derived using Voigt– Reuss– Hill (VRH) approximation, [48] and the results are presented in Fig.  5. The VRH bulk modulus is defined as the average B = (BV + BR)/2 between the Voigt upper and Reuss lower bounds

The VRH shear modulus G = (GV + GR)/2 is taken from the average between the Voigt upper GV and Reuss lower GR bounds

From Fig.  5, three minimum value of the bulk B and shear G modulus are located at x = 0.2, 0.464, and 0.7. The bulk B and shear G modulus at x = 0.7 are the smallest among the three minimum values. The corresponding d33 is also the smallest, as shown in Fig.  2. This implies that the small bulk and shear modulus would not necessarily give rise to a high piezoelectric response. Otherwise, the exceptionally high bulk and shear modulus are also not conducive to a high piezoelectric response. Therefore, the moderate bulk and shear modulus are better for the piezoelectric response.

Fig.  5. Bulk and shear modulus of the orthorhombic K1− xNaxNbO3 as a function of x.

3.6. Partial density of state

The first-principles studies have suggested that the orbital hybridization between the Pb 6s state and O 2p states plays a crucial role for large ferroelectricity in tetragonal PbTiO3.[49] The prediction is later verified experimentally, [50] i.e., the Pb– O bonds in tetragonal PbTiO3 show rather strong covalency. In order to investigate the interaction between Nb– O and K– O bonds, we plot the partial density of state (PDOS) of the Nb, K, and O orbitals for the orthorhombic K1− xNaxNbO3 with x = 0.5. As shown in Fig.  6, the K 2s and 2p, O 2p orbits are mainly located in the range from − 5.6 to 0  eV, indicating that the K– O bonds are mainly ionic. The calculated Born effective charge of K is 1.13, which is close to its oxidation state, further confirming the characteristic of ionic interaction between K and O ions. The Nb 4d orbits mixed with O 2p orbits are located in the valence band and conduction band, suggesting that the Nb– O bonds are both ionic and covalent. The calculated average Born effective charge of Nb is 8.79, larger than its oxidation state (+ 5), indicating that the Nb– O bonds are not purely ionic. Thus, the high piezoelectric for orthorhombic K1− xNaxNbO3 with x = 0.5 is related to the Nb– O covalency. This agrees with the recent theoretical study.[51]

Fig.  6. Partial density of state of the orthorhombic K1− xNaxNbO3.

3.7. Energy barrier on polarization rotation paths

To understand the physical mechanism of the high piezoelectric response for the orthorhombic K1− xNaxNbO3, we calculate the energy barrier on polarization rotation paths along the [001] direction through a simple approach: E(θ ) = E0(θ ) − E0(0), where E(θ ) is the energy barrier with θ being the angle of polarization rotation, E0(θ ) is the total energy of the orthorhombic K1− xNaxNbO3, and E0(0) denotes the total energy with θ = 0° . The adopted model for the energy barrier calculations is illustrated in Fig.  7, and the results are depicted in Fig.  8. As shown in Fig.  8, the energy barrier of the orthorhombic K1− xNaxNbO3 with x = 0.0 and x = 0.5 increases with the increase of the angle θ . However, the energy barrier of the orthorhombic K1− xNaxNbO3 with x = 0.5 increases more slowly than the one with x = 0.0, indicating that at MPB, the phase instability of the orthorhombic K1− xNaxNbO3 with x = 0.5 is enhanced and the polarization rotation becomes easier than the pure KNO, i.e., the flat free energy profile is the most important factor for inducing a high piezoelectric response.

Fig.  7. Schematic of polarization rotation of the orthorhombic K1− xNaxNbO3 resulting from [001] electric field.

Fig.  8. Energy barrier of the orthorhombic K1− xNaxNbO3 with x = 0.5 and 0.0 on polarization rotation paths along the [001] direction.

4. Conclusions

First-principles density functional theory within the virtual crystal approximation is applied to investigate the morphotropic phase boundary and piezoelectric properties of K1− xNaxNbO3. The structure parameter, bulk modulus, band gap, and piezoelectric coefficients have been calculated. The MPB in K1− xNaxNbO3 is predicted to be located at x = 0.521, which is in good agreement with the experimental data. The dependence of the longitudinal piezoelectric coefficient on the polarization orientation proves that the highest value of the piezoelectric coefficient of the orthorhombic K1− xNaxNbO3 with x = 0.5 is close to the [001] direction. According to the comparison of the energy barrier of polarization rotation from [011] to [001] direction, we prove that the polarization rotation for K0.5Na0.5NbO3 becomes easier at MPB, and a flat free-energy is the most important factor for high piezoelectricity. The present calculation also provides a way to predict the MPB in lead-free piezoelectric materials.

Reference
1 Sun E and Cao W 2014 Prog. Mater. Sci. 65 124 DOI:10.1016/j.pmatsci.2014.03.006 [Cited within:2] [JCR: 23.194]
2 Gao Y, Zhang J, Zong X, Wang C and Li J 2010 J. Appl. Phys. 107 074101 DOI:10.1063/1.3369435 [Cited within:1] [JCR: 0.71]
3 Guo Y, Kakimoto K I and Ohsato H 2004 Solid State Commun. 129 279 DOI:10.1016/j.ssc.2003.10.026 [Cited within:1] [JCR: 1.534]
4 Hollenstein E, Davis M, Damjanovic D and Setter N 2005 Appl. Phys. Lett. 87 182905 DOI:10.1063/1.2123387 [Cited within:1] [JCR: 3.794]
5 Shrout T R and Zhang S J 2007 J. Electroceram. 19 113 DOI:10.1007/s10832-007-9047-0 [Cited within:1] [JCR: 1.143]
6 Zhang J, Zong X, Wu L, Gao Y, Zheng P and Shao S 2009 Appl. Phys. Lett. 95 022909 DOI:10.1063/1.3182725 [Cited within:1] [JCR: 3.794]
7 Saito Y, Takao H, Tani T, Nonoyama T, Takatori K, Homma T, Nagaya T and Nakamura M 2004 Nature 432 84 DOI:10.1038/nature03028 [Cited within:1] [JCR: 38.597]
8 Zhang S and Li F 2012 J. Appl. Phys. 111 031301 DOI:10.1063/1.3679521 [Cited within:3] [JCR: 0.71]
9 Wang R, Xie R-J, Hanada K, Matsusaki K, Band o H, Sekiya T and Itoh M 2006 Ferroelectrics 336 39 DOI:10.1080/00150190600695321 [Cited within:2] [JCR: 0.415]
10 Tennery V J and Hang K W 1968 J. Appl. Phys. 39 4749 DOI:10.1063/1.1655833 [Cited within:1] [JCR: 0.71]
11 Zhang B P, Li J F, Wang K and Zhang H 2006 J. Am. Ceram. Soc. 89 1605 DOI:10.1111/jace.2006.89.issue-5 [Cited within:1] [JCR: 2.107]
12 Jaeger R E and Egerton L 1962 J. Am. Ceram. Soc. 45 209 DOI:10.1111/jace.1962.45.issue-5 [Cited within:2] [JCR: 2.107]
13 Dai Y J, Zhang X W and Chen K P 2009 Appl. Phys. Lett. 94 042905 DOI:10.1063/1.3076105 [Cited within:2] [JCR: 3.794]
14 Li Q, Huang D H, Cao Q L and Wang F H 2013 Chin. Phys. B 22 037101 DOI:10.1088/1674-1056/22/3/037101 [Cited within:1] [JCR: 1.148] [CJCR: 1.2429]
15 Luo F, Cheng Y, Cai L C and Chen X R 2013 J. Appl. Phys. 113 033517 DOI:10.1063/1.4776679 [Cited within:1] [JCR: 0.71]
16 Ao B-Y, Shi P, Guo Y and Gao T 2013 Chin. Phys. B 22 037103 DOI:10.1088/1674-1056/22/3/037103 [Cited within:1] [JCR: 1.148] [CJCR: 1.2429]
17 Guo Y, Ai J J, Gao T and Ao B Y 2013 Chin. Phys. B 22 057103 DOI:10.1088/1674-1056/22/5/057103 [Cited within:1] [JCR: 1.148] [CJCR: 1.2429]
18 Hu C, Wang F, Xia C, Zheng Z and Ren W 2011 Mod. Phys. Lett. B 25 333 DOI:10.1142/S0217984911025705 [Cited within:1] [JCR: 0.479]
19 Gonze X, Rignanese G M, Verstraete M, Beuken J M, Pouillon Y, Caracas R, Jollet F, Torrent M, Zerah G, Mikami M, Ghosez P, Veithen M, Raty J Y, Olevano V, Bruneval F, Reining L, Godby R, Onida G, Hamann D R and Allan D C Z 2005 Kristallogr. 220 558 DOI:10.1524/zkri.220.5.558.65066 [Cited within:1] [JCR: 1.241]
20 Gonze X, Amadon B, Anglade P M, Beuken J M, Bottin F, Boulanger P, Bruneval F, Caliste D, Caracas R and Cote M 2009 Comput. Phys. Commun. 180 2582 DOI:10.1016/j.cpc.2009.07.007 [Cited within:1] [JCR: 3.078]
21 Gonze X and Lee C 1997 Phys. Rev. B 55 10355 DOI:10.1103/PhysRevB.55.10355 [Cited within:1]
22 Hamann D, Wu X, Rabe K M and Vand erbilt D 2005 Phys. Rev. B 71 035117 DOI:10.1103/PhysRevB.71.035117 [Cited within:1]
23 Bellaiche L and Vand erbilt D 2000 Phys. Rev. B 61 7877 DOI:10.1103/PhysRevB.61.7877 [Cited within:1]
24 Ramer N J and Rappe A M 2000 Phys. Rev. B 62 R743 DOI:10.1103/PhysRevB.62.R743 [Cited within:1]
25 Liu S Y, Shao Q S, Yu D S, Lv Y K, Li D J, Li Y and Cao M S 2013 Chin. Phys. B 22 17702 DOI:10.1088/1674-1056/22/1/017702 [Cited within:1]
26 Deng Y, Wang R Z, Xu L C, Fang H, Yang X, Yan H and Chu P K 2011 Appl. Phys. A 104 1085 DOI:10.1007/s00339-011-6375-3 [Cited within:1] [JCR: 1.545]
27 Haumont R, Dkhil B, Kiat J, Al-Barakaty A, Dammak H and Bellaiche L 2003 Phys. Rev. B 68 014114 DOI:10.1103/PhysRevB.68.014114 [Cited within:1]
28 Geneste G, Kiat J M and Malibert C 2008 Phys. Rev. B 77 052106 DOI:10.1103/PhysRevB.77.052106 [Cited within:1]
29 Hartwigsen C, Goedecker S and Hutter J 1998 Phys. Rev. B 58 3641 DOI:10.1103/PhysRevB.58.3641 [Cited within:1]
30 Wu X, Vand erbilt D and Hamann D 2005 Phys. Rev. B 72 035105 DOI:10.1103/PhysRevB.72.035105 [Cited within:1]
31 Wu Z, Cohen R and Singh D 2004 Phys. Rev. B 70 104112 DOI:10.1103/PhysRevB.70.104112 [Cited within:1]
32 Shirane G, Newnham R and Pepinsky R 1954 Phys. Rev. 96 581 DOI:10.1103/PhysRev.96.581 [Cited within:2] [JCR: 6.583]
33 Wan L, Nishimatsu T and Beckman S 2012 J. Appl. Phys. 111 104107 DOI:10.1063/1.4712052 [Cited within:1] [JCR: 0.71]
34 Wiesendanger E 1973 Ferroelectrics 6 263 DOI:10.1080/00150197408243977 [Cited within:1] [JCR: 0.415]
35 Shigemi A and Wada T 2005 Jpn. J. Appl. Phys. 44 8048 DOI:10.1143/JJAP.44.8048 [Cited within:1] [JCR: 1.067]
36 Hewat A 1973 J. Phys. C: Solid State Phys. 6 2559 DOI:10.1088/0022-3719/6/16/010 [Cited within:2]
37 Katz L and Megaw H 1967 Acta Cryst. 22 639 DOI:10.1107/S0365110X6700129X [Cited within:1]
38 Kalinichev A G, Bass J D, Zha C S, Han P D and Payne D A 1993 J. Appl. Phys. 74 6603 DOI:10.1063/1.355099 [Cited within:1] [JCR: 0.71]
39 Günter P 1977 Jpn. J. Appl. Phys. 16 1727 DOI:10.1143/JJAP.16.1727 [Cited within:1] [JCR: 1.067]
40 Jaffe B 2012 Piezoelectric CeramicsVol.  3 (Elsevier) [Cited within:1]
41 Nakamura Y, Kawai M, Azuma M and Shimakawa Y 2010 Jpn. J. Appl. Phys. 49 051501 DOI:10.1143/JJAP.49.051501 [Cited within:1] [JCR: 1.067]
42 Hiruma Y, Nagata H and Takenaka T 2008 J. Appl. Phys. 104 124106 DOI:10.1063/1.3043588 [Cited within:1] [JCR: 0.71]
43 Royles A, Bell A, Jephcoat A, Kleppe A, Milne S and Comyn T 2010 Appl. Phys. Lett. 97 132909 DOI:10.1063/1.3490235 [Cited within:1] [JCR: 3.794]
44 Kaoru M Makoto K Masaki A Hiroshi F 2010 Jpn. J. Appl. Phys. 49 09ME07 DOI:10.1143/JJAP.49.09ME07 [Cited within:1] [JCR: 1.067]
45 Davis M, Budimir M, Damjanovic D and Setter N 2007 J. Appl. Phys. 101 054112 DOI:10.1063/1.2653925 [Cited within:1] [JCR: 0.71]
46 Li F, Zhang S, Xu Z, Wei X, Luo J and Shrout T R 2010 J. Am. Ceram. Soc. 93 2731 DOI:10.1111/jace.2010.93.issue-9 [Cited within:1] [JCR: 2.107]
47 Liu X, Zhang S, Luo J, Shrout T R and Cao W 2010 Appl. Phys. Lett. 96 012907 DOI:10.1063/1.3275803 [Cited within:1]
48 Hill R 1963 J. Mech. Phys. Solid. 11 357 DOI:10.1016/0022-5096(63)90036-X [Cited within:1] [JCR: 3.406]
49 Cohen R E 1992 Nature 358 136 DOI:10.1038/358136a0 [Cited within:1] [JCR: 38.597]
50 Kuroiwa Y, Aoyagi S, Sawada A, Harada J, Nishibori E, Takata M and Sakata M 2001 Phys. Rev. Lett. 87 217601 DOI:10.1103/PhysRevLett.87.217601 [Cited within:1] [JCR: 7.943]
51 Shi J, Grinberg I, Wang X and Rappe A M 2014 Phys. Rev. B 89 094105 DOI:10.1103/PhysRevB.89.094105 [Cited within:1]