EElectrostriction coefficient of ferroelectric materials from ab initio computation
Z. Jiang,
1, 2
R. Zhang, F. Li, L. Jin, N. Zhang, D. Wang, ∗ and C.-L. Jia
1, 4 Electronic Materials Research Laboratory–Key Laboratory of the Ministry of Education andInternational Center for Dielectric Research, Xi’an Jiaotong University, Xi’an 710049, China Physics Department and Institute for Nanoscience and Engineering,University of Arkansas, Fayetteville, Arkansas 72701, USA Department of Physics, State Key Laboratory of Photoelectric Technology andFunctional Materials (Cultivation Base), Northwest University, Xi’an 710069, China Peter Grnberg Institute and Ernst Ruska Center for Microscopy andSpectroscopy with Electrons, Research Center Jlich, D-52425 Jlich, Germany
Electrostriction is an important material property that characterizes how strain changes with the developmentof polarization inside a material. We show that ab initio techniques developed in recent years can be exploitedto compute and understand electrostriction of ferroelectric materials. Here, electrostriction coefficients of fer-roelectric BaTiO , PbTiO , as well as dielectric BaZrO , are obtained and analyzed. Possible causes of thedifference between experimental and numerical results are discussed. We also identified that relative displace-ments between certain ions at a given polarization could be a good indicator of a material’s electrostrictionproperty. I. INTRODUCTION
Electrostriction effect describes how a material deforms asa polarization develops inside it. Being one of the nonlin-ear phenomena particularly important in ferroelectric mate-rials, it is the physical basis for their application in actuators,sensors, phase shifters, and sonar projectors. The effective-ness of electrostriction is often gauged by the electrostrictivecoefficients Q ijkl ( i = x, y, z represents Cartesian coordi-nate) that characterizes how the strain of a material changesquadratically with its polarization, i.e., x ij = Q ijkl P k P l , where x ij represents strain tensors (in the following we willuse Voigt notation for strain tensor, i.e., x i represents strainwith i = 1 . . . ). Electrostriction is a four-rank polar tensor,which can be observed in essentially all crystals regardless oftheir symmetries.Development of materials with large electrostriction is animportant task in material science research. For instance, itwas suggested that the large piezoelectric coefficient of thelead-free ceramic system BZT-BCT is originated from itslarge electrostriction coefficients. It is obvious that accuratelyobtaining electrostriction coefficient is a key step in such re-search. Experimentally, there are several ways to obtain elec-trostriction coefficient Q : (i) Direct approaches, such as mea-suring the change in strain as a function of polarization, or (ii)Indirect approaches, such as measuring the change of permit-tivity (via change of capacitance) as a function of mechanicalstress. Computationally, the electrostriction coefficient can beobtained by applying an electric field to the target material,and compute the change of polarization and strain at each stepas the electric field increases.However, for ferroelectric materials, ab initio computationof electrostriction coefficients is not as straightforward as onemay expect even after the modern theory of polarization wasestablished. There are two main obstacles: (i) The appli-cation of an electric field in a periodic structure needs to betreated carefully; (ii) Electrostriction is often obtained byassuming the target material is in high symmetry phases, ei-ther because it is above the phase transition temperature, T C (above which many perovskite ferroelectrics are in the para-electric phase), or because polycrystals or ceramics are used,in which the overall phase is isosymmetric. Such a situationcauses the second difficulty in ab initio calculation becausethe mimicked high symmetry phase is unstable for ferroelec-tric materials (e.g., the paraelectric cubic phase for BaTiO isnot the ground state). Such an unstable system can hardly beused in ab initio computation as it quickly jumps to other morestable phases (e.g., the tetragonal phase in some perovskites)in structural relaxations, forbidding us to set the system at agiven polarization and track how the strain changes.Recently, these two obstacles have been overcome dueto the progress in ab initio computation methods, whichmakes the calculation of electrostriction coefficients easier.The first breakthrough came when Souza et al. and Umari et al . implemented the code for applying electric fieldin periodic structures. With this advancement, it is possi-ble to compute in general how a material responds to ap-plied electric field. As a matter of fact, such method offixing the electric field E has been applied widely in thestudy of dielectric, piezoelectric, and ferroelectric behavior ofmaterials. To address the second problem of the unstablehigh-symmetry phase, Dieguez et al . developed a methodthat is able to impose a constant polarization, P (fixed- P method). In this way, one can increase the polarization ofthe system gradually and track the change of strain, whichpaves the way for numerically obtaining electrostriction co-efficient. However, further developments show that a betterway is to impose the constraint of a constant electric displace-ment, D , (fixed- D method). Among other reasons, this ap-proach addresses a few subtle points that are not fully ap-preciated in previous studies. For instance, (i) The fixed- D method corresponds to imposing an open-circuit electricalboundary conditions, which is especially useful for studyingferroelectric capacitors, superlattices and metal-oxideinterfaces; (ii) This approach made it possible to mimicthe unstable cubic phase of ferroelectric materials, and extractuseful parameters (e.g., electrostriction coefficient). More im-portantly, with the implementation of such a method in the a r X i v : . [ c ond - m a t . m t r l - s c i ] J un ABINIT software package, it enables more researchers toobtain electrostriction coefficients via ab initio computation.Here we use the fixed- D method to study electrostric-tion in prototypical ferroelectric materials BaTiO (BTO) andPbTiO (PTO). As a comparison, we also present results fromBaZrO (BZO), which has the cubic phase as its ground state.Our numerical results indeed shows that this approach can ob-tain electrostriction coefficients consistent with experimentalresults, therefore providing a valuable predictive power in dis-covering materials with better electrostrictive properties, andan insight to the factors that determine electrostriction coeffi-cients. II. METHOD
The ab initio calculations were performed using the open-source ABINIT software package. Our calculations wereperformed within density-functional theory in the local-density approximation (LDA) and the projector-augmented-wave (PAW) method is used along with pseudopotentialsimplemented in the GBRV package. Ba
5s 5p 5d 6s , Ti
3s 3p 4s 3d , Pb
5d 6s 6p , Zr
4s 4p 4d 5s and O
2s 2p or-bitals were treated as valence orbitals. A plane wave basiswith kinetic energy cutoff of 40 Hartree was used to ensurethe convergence in all the calculations. k -point sampling of × × Monkhorst-Pack grid was used for all the 5-atomunit cell. The atomic coordinates of the five-atom unit cellwere relaxed until all atomic-force components were smallerthan − Hatree/Bohr, and the cell size-and-shape was varieduntil all stress components were below − Hatree/Bohr . III. RESULTS AND DISCUSSION
The fixed- D calculation implemented by Stengel et al . is exploited in this work, where we gradually increase thereduced electric displacement field d along the [001] direc-tion. This method enables the relaxation of both atom po-sitions and cell size-and-shape of a given crystal at a given d , the flux of the electric displacement field, defined as d = a × a D/ (4 π ) , where a and a are the in-plane latticevectors, and D is the electric displacement field. At each stepof the calculation process, the internal energy of the system, U , the polarization of the system, as well as the relaxed cell(determined by a , a , and a ) are computed. After the com-putation, we have the relation between strain and polarization,which can be used to compute electrostriction coefficients.Note that for the calculation of strain, we use the initial re-laxed cubic phase a as the reference. The system is initiallyset to be cubic ( P m ¯3 m ), but with the application of the fixed-displacement-field D along the [001] direction, it becomestetragonal ( P mm ). The electric field, polarization, internalenergy, ionic positions and lattice parameters are all variedat a given value of d . To obtain electrostriction coefficientsof each system (BTO, PTO and BZO), a series of computa-tions are performed at given electric displacement, d , alongthe z direction. Internal energy, polarization and strain are then obtained in each process, and related to the constrainedvariable, d . We therefore obtain the out-of-plane strain, x ,and the in-plane strain x as functions of d , which are thenfitted to calculate electrostriction coeffcients. Moreover, inthis process we also track ion displacements versus d , whichprovides some insight to the electrostriction effects. We notethat to define response properties it is important to specify theboundary condition or constraint used in their definition. Here our calculation of electrostriction coeffcients is under theconstraint of zero stress.
A. Ferroelectric materials BaTiO and PbTiO -8081624 E n e r g y ( m e V ) (a) -70-60-50-40-30-20-100 (e) P o l a r i z a ti o n ( C / m ) d ( e ) (b) d ( e ) (f) x (c) (g) -0.006-0.004-0.00200 0.1 0.2 0.3 0.4 0.5 x Polarization (C/m ) (d) -0.016-0.012-0.008-0.00400 0.2 0.4 0.6 0.8 1 1.2 Polarization (C/m ) (h)BTO PTO Figure 1: Internal energy (a), polarization (b), out-of-plane strain (c),and in-plane strain (d) of BTO, as a function of d , the reduced electricdisplacement [(a), (b)] and polarization [(c), (d)]. Internal energy(e), polarization (f), out-of-plane strain (g), and in-plane strain (h) ofPTO, as a function of d [(e), (f)] and polarization [(g), (h)]. We first show the results of internal energy U , strain inthe [001] ( x ) and [100] direction ( x ) versus d , which liesin the [001] direction. For BTO, we start from the relaxedcubic structure (the lattice constant a c = 3 . ˚A), and grad-ually increase d in steps of . e until it reaches 0.48 e ( − e is the electron charge). For PTO, we increase d in steps of . e until it reaches . e . The structure was fully relaxedwith respect to both ionic positions and lattice parameters ateach d value for both BTO and PTO. Strains are calculatedusing x = ( a c − a c ) /a c , where a c is relaxed lattice con-stant at a given d . Numerical results for BTO are shown inFig. 1. In Fig. 1(a), we plot the internal energy as a func-tion of d . Clearly, there is an energy minimum at d (cid:39) . e ,indicating that the cubic BTO is unstable and such instabilitywill introduce a ground state with nonzero d . In Fig. 1(b),we plot the polarization versus d , and find that the polariza-tion changes linearly with d . The reason is that in the regionof interest, it is shown that the ε E ( ε is the vacuum per-mittivity) is much smaller than P . Therefore P dominatesthe relation D = ε E + P , ensuring this linear relation. Figures 1(c) and 1(d) show that the strain versus polarizationcurves are quadratic, which implies that electrostriction, notpiezoelectricity, dominates the mechanical response to polar-ization. To quantify the electrostriction effect, we employ theexpression x ij = Q ijkl P k P l to fit the curves in Figs. 1(c) and1(d). More precisely, we use the equation x = a + a P for fitting. It is found that Q = 0 . m /C , which is quiteclose to the experimental value of Q = 0 . m /C , and Q = − . m /C , which is comparable to the experi-mental value of Q = − . m /C . The piezoelectricvoltage coefficients g and g can be obtained from the rela-tion g ijk = 2 Q ijkl P s . For BTO, we obtained g = 0 . Vm/N and g = − . V m/N, which are comparable to theexperimental values of g = 0 . V m/N and g = − . V m/N. This comparison supports the idea that electrostric-tion may be used to estimate the electromechanical properties(e.g., piezoelectric coefficients) of ferroelectric materials withspontaneous polarization.Numerical results for PTO is also shown in Fig. 1, wherewe plot various properties of PTO under different d . Fig-ure 1(e) shows that PTO also has an energy minimum, whichagrees well with previous studies in terms of the double-well potential and the position of the energy minimum. Fig-ure 1(e) indicates that the depth of energy well of PTO( ∼ meV) is much deeper than that of BTO ( ∼ meV)and the energy minimum corresponds to a larger polariza-tion ( P = 0 . C/m versus P = 0 . C/m in BTO.). InFigs. 1(g) and 1(h), we show the out-of-plane and in-planestrain versus polarization, and found Q = 0 . m /C and Q = − . m /C by fitting the two quadratic curves.Comparable experimental values of Q = 0 . m /C and Q = − . m /C have been reported. We note that thefitting range and whether including higher order terms can po-tentially affect the final results. For instance, if we include the P term to fit Figs. 1(g) and 1(h), we find Q = 0 . m /C and Q = − . m /C . B. Dielectric material BaZrO Having examined ferroelectric materials, we now move toBZO, which is a dielectric material with perovskite structurethat has larger lattice constant, smaller thermal expansion, andsmaller dielectric permittivity than BaTiO . Experimen-tal and theoretical results indicate that BZO remains paraelec-tric at low temperatures, which means it is not piezoelectric,but can have electrostriction effect. Figure 2 shows the in-ternal energy U and polarization versus d , and strains versuspolarization. Figure 2(a) shows that, unlike BTO and PTO,the internal energy U always increases with d . This is con- E n e r g y ( m e V ) (a) P o l a r i z a ti o n ( C / m ) d ( e ) (b) x (c) -0.005-0.004-0.003-0.002-0.00100 0.1 0.2 0.3 0.4 x Polarization (C/m ) (d) Figure 2: Internal energy (a), polarization (b), out-of-plane strain(c), and in-plane strain (d) of BZO, as a function of d [(a), (b)] andpolarization [(c), (d)]. sistent with the the fact that cubic BZO has no ferroelectricinstability. The electrostriction coefficients we obtained are Q = 0 . m /C and Q = − . m /C . Since di-rect experimental value on Q and Q are not available, wealso use Q h = Q + 2 Q to calculate the hydrostatic elec-trostrictive coefficients and find Q h = 0 . m /C . As acomparison, the experimental value of Q h is found to be . m /C in Ref. where the technique of measuring capacitancechange under hydrostatic pressure is used. The difference be-tween experimental Q h and the value we find here may berelated to the facts that (i) We did not use the constant vol-ume constraint, and/or (ii) There are anti-phase motions ofoxygen octahedra found in BZO. It is interesting to notethat the numerical Q of BZO obtained here is larger thanPTO (but smaller than BTO). However, it is worth noting thatsince BZO is not ferroelectric, it may be difficult to convert thelarge electrostriction effect into piezoelectricity, unlike BTOor PTO in which spontaneous ferroelectricity automaticallybreak the centrosymmetry without external electric field.Since BZO does not have a ferroelectric instability, it isalso possible to obtain its electrostriction coefficient by di-rectly applying an electric field. Here, we apply an electricfield, E , along the [001] direction. We applied electric fieldup to E = 1 . × V/m. At every step of the appliedelectric field, the structure of BZO, including its lattice con-stant, atomic positions, and cell shape, are optimized, similarto what is done in the fixed- D approach. Once the optimizedstructure of BZO at a given E is obtained, its polarization iscalculated using the Berry phase approach. In this way, theelectrostriction coefficients are found to be Q = 0 . m /C and Q = − . m /C , which agree well withthe results obtained via the fixed- D method. C. Ion displacement
We also obtained the displacement of each ions versus po-larization in Fig. 3(a), where we show the reduced coordi-nates of BTO versus the polarization. A prominent feature ofthis figure is that all the atoms shift linearly with polarization. -0.03-0.02-0.0100.010.020.030.040 0.1 0.2 0.3 0.4 0.5 R e du c e d c oo r d i n a t e s (a) -0.06-0.04-0.0200.020.040.060.080 0.2 0.4 0.6 0.8 1 (b) -0.03-0.02-0.0100.010.020.030.040.050 0.1 0.2 0.3 0.4 R e du c e d c oo r d i n a t e s Polarization [C/m ] (c) -0.03-0.02-0.0100.010.020.030.040 0.1 0.2 0.3 0.4 Polarization [C/m ] (d) ZrBaOO ZrBaOOTiPbOOTiBaOO
Figure 3: The out-of-plane displacements of each ion in BTO, PTOand BZO from zero-field positions change with polarization obtainedwith the fixed- D method are shown in panels (a), (b), and (c), respec-tively. For BZO, the method of directly applying an electric field isalso shown [Panel (d)]. As expected, Ba and Ti ions (being positively charged) movealong the direction of the reduced electric displacement field d , while O ⊥ and O (cid:107) ions (note the Ti-O (cid:107) and Ti-O ⊥ bondsare parallel and perpendicular to the z direction, respectively),which are negatively charged, move in the opposite direction.Together, they induce a polarization along the z direction. Wefurther note that the displacements of Ti, Ba, O ⊥ , O (cid:107) obtainedhere with the fixed- D method are compatible with that ob-tained from force-constant matrix, where the eigenvector isgiven by (cid:0) ξ Ba , ξ Ti , ξ O ⊥ , ξ O (cid:107) (cid:1) = (0 . , . , − . , − . . Here we find from the slope of the lines shown in Fig. 3(a)that the eigenvector is (0 . , . , − . , − . . The close-ness of these two results and the linear correlation shown inFig. 3(a) further validate the accuracy of the effective Hamil-tonian approach in the study of BTO. It is worth notingthat the eigenvector shown in Ref. is obtained by diagonaliz-ing the force-constant matrix. With a series of similar fixed- D calculations, we believe the static field-induced ion displace-ments agree better with the force-constant matrix eigenvectorsthan the dynamical matrix eigenvector. Figure 3(b) shows thedisplacement of each ions in PTO as a function of the polar-ization. Similar to BTO, the displacements of the ions arelinearly dependent on P , and consistent with the eigenvectorobtained in frozen phonon calculation and experiments. However, one can see that in PTO, the A-site ion Pb movesin a larger distance than the B-site ion, Ti, which is oppositeto what happens in BTO. Since the displacement of Pb ionis about twice as large as that of Ti, this result also indicatesthat in PTO Pb has larger contribution to polarization than Ti.Since Pb is outside oxygen octahedrons, it likely has moreroom to move without distorting the unit cell along z , unlikeTi ions. In other words, to achieve the same polarization, Ti inPTO has a smaller displacement than BTO, and thus a smaller deformation, which results in a smaller Q . Such a differencemay explain why Q of PTO is smaller than that of BTO. Wethus suggest that B-site driven ferroelectrics may in generalhave better electrostriction than A-site driven ones. Figure3(b) also shows an anomalous change around P = 0 . C / m .The inverse capacitance of PTO was also found to have sucha change, which is discussed in Ref. . Figure 4: The B -O (cid:107) bond length (in reduced coordinates) at P =0 . C / m and Q of BTO, KNO (KNbO ), STO (SrTiO ), BZOand PTO. The ion displacements in BZO versus d also has a linearbehavior as shown in Fig. 3(c) (fixed-D) and Fig. 3(d) (fixed- E ), where Ba and Zr ions move along d or E , while O ionsmove in the opposite direction. It is worth noting that, inBZO, the A-site atom, Ba, shifts more than the B-site atom,Zr, which is opposite to what happens in BTO, but similarto PTO. Comparison of Fig. 3(c) and Fig. 3(b) to Fig. 3(a)also reveals another important difference between BTO andBZO/PTO, i.e., (cid:12)(cid:12) ξ O (cid:107) (cid:12)(cid:12) > | ξ O ⊥ | in BTO while the opposite istrue in BZO/PTO. The larger ξ O ⊥ in BZO/PTO is likely dueto the larger displacement of the A-site atom (Ba in BZO andPb in PTO) than that of the B-site atom ( Zr in BZO and Ti inPTO) as the displacement A-site ion may cause neighboringoxygen octahedrons to tilt. The importance of such differenceon electrostriction becomes apparent when the bond lengthof the B-site atom and O (cid:107) along with Q is summarized inFig. 4 for five different perovskites (Table I summarizes thenumerical results of Q , Q , g , g , C , C , C andspontaneous polarization P , the results of elastic constantsand spontaneous polarization are in general agree with exist-ing computational and experimental results. ), whichindicates that Q has positive correlation with the B -O (cid:107) bondlength. In other words, longer B -O (cid:107) means larger theoreticalelectrostriction coefficients when other things (e.g. polariza-tion) being equal. This relation can be understood by realizingthe fact that the B -O (cid:107) bond length characterizes and, to someextent, determines the unit cell deformation along z , whichin turn decides the magnitude of Q . In practice, such anindicator can potentially facilitate the investigation of elec-trostriction by paying more attention to the B -O (cid:107) bond length. BaTiO KNbO SrTiO BaZrO PbTiO Q (m /C ) 0.137 0.106 0.088 0.094 0.065 Q (m /C ) -0.026 -0.028 -0.020 -0.025 -0.012 g (V m/N) 0.071 0.066 - - 0.108 g ( V m/N) -0.014 -0.017 - - -0.020 C (GPa) 322.7 422.6 355.2 332.4 324.5 C (GPa) 110.3 73.6 101.6 79.9 124.0 C (GPa) 126.6 94.1 113.9 86.6 101.1 P s (C/m ) 0.26 0.31 0.00 0.00 0.83B-O (cid:107) Q , Q , piezoelectric volt-age coefficients g and g , cubic phase elastic constants C , C , C , spontaneous polarization P s , and the B-O (cid:107) bond length (in re-duced coordinates) at P = 0 . C/m for BTO, KNO, STO, BZO andPTO. D. Discussion
With the aid of ab initio computation in general and thefixed- D method in particular, we now have an additional per-spective to understand electrostriction. As we know, elec-trostriction represents a relation between displacements ofatoms (that determines the polarization) and structural param-eters (that determines strain). It is often taken for granted thatthe polarization is caused by an external electric field. How-ever, it has been shown that the relation between P and E israther complex for ferroelectric materials. For instance,to constrain the system to a given P , E need to be negative(i.e., opposite in the direction of P ) to stabilize the system ata given configuration. The reason is not hard to understand– it is the internal electric field arising from electric dipolesthat drives the ferroelectric instability, which in turn resultsin the spontaneous polarization. For ferroelectric materi-als, it is also found from ab initio computation that for P val-ues of our interest, the required electric field is much smallerthan P , i.e., ε E (cid:28) P , which explains why the relation x ij = Q ijkl P k P l is adopted more than x ij = M ijkl E k E l tocharacterize electrostriction in ferroelectric materials.As we have shown, the fixed- D approach is effective toobtain electrostriction coefficients via ab initio computations.This computational approach is obviously useful for discover-ing materials or structures with desired electrostrictive prop-erties. In order to employ this method properly, it is importantto understand potential causes of difference between experi-mental work and ab initio calculation. First, in many exper-iments on electrostriction, polycrystals or ceramics are used,which has two implications: (i) The effect of domains, domainwalls, and/or grain boundaries in real materials cannot be re-flected in ab initio calculation with a rather small number ofunit cells; and (ii) The sample may not be in its high symmetry phase due to mechanical stress or low temperature (lower than T C ) – for such a situation, a direct comparison is not ideal.More sophisticated ab initio computation of other phases (e.g.the tetragonal phase) and some statistics on the results maybe mandatory. Second, in comparison with experimental val-ues, one needs to be aware of some electric boundary condi-tions between different ab initio computations. For instance,the fixed- E computation corresponds to a closed circuit elec-tric boundary condition with a biased voltage and assumesthat E is continuous through the system, while the fixed- D method corresponds to an open circuit electric boundary con-dition and continuity of D . This subtlety warrants care-ful consideration before adopting one or the other approach.For complex non-uniform structures, we suggest to use thefixed- D approach, which is particularly important for layeredstructures, such as superlattices or metal/insulator interfaces,as pointed in Ref. . This practice should be considered evenfor non-ferroelectric materials where applying electric fieldis not a problem. Third, more complex structures (e.g. theanti-phase oxygen octahedron tiltings in BZO) could affectthe accuracy of ab initio results. However, we note that in-troducing such complex crystal structures will substantiallyslow down the calculation with the fixed- D approach. Finally, ab initio computations assume zero temperature, therefore theoptimized lattice constant may not corresponds to the experi-mental values at finite temperatures. Experimentally verifiedvalues (e.g., lattice constant) may be needed to reduce differ-ences in numerical values. IV. CONCLUSION
In summary, we have shown that the fixed- D approach canbe applied to obtain electrostriction coefficients of ferroelec-tric materials. By investigating three materials (BTO, PTO,and BZO), we find the fixed- D approach provide a powerfulmethod to compute and understand electrostrictive effects ofvarious materials, and is potentially useful for future predic-tions of materials or designs. With the aid of this approach,we have also shown that, given a perovskite material, the bondlength of the B-site atom and the O (cid:107) ion is an important indi-cator of the magnitude of its electrostriction, which could be auseful criterion in screening a large number of perovskites tofind those with good electrostrictive effects. In addition, pos-sible difference between electrostriction coefficients obtainedtheoretically and experimentally are discussed in detail. Wethus hope that the approach we demonstrat here will be widelyused in discovering and exploring materials with better elec-trostrictive effects. Acknowledgments
This work is financially supported by the National NaturalScience Foundation of China (NSFC), Grant No. 51390472,11574246, and National Basic Research Program of China,Grant No. 2015CB654903. F.L. acknowledges NSFC GrantNo. 51572214. Z.J. also acknowledges support from ChinaScholarship Council (CSC No. 201506280055) and the 111Project (B14040). We thank Shanghai Supercomputer Center for providing resources for ab initio computations. ∗ Electronic address: [email protected] K. Uchino, Piezoelectric Actuators and Ultrasonic Motors(Kluwer Academic Publishers, Boston, 1996). K. Uchino, S. Nomura, L. E. Cross, R. E. Newnham, and S. J.Jang, J. Mater. Sci. , 569 (1981). S.-E. Park and W. Hackenberger, Curr. Opin. Solid State Mater.Sci. , 11 (2002). R. E. Newnham, V. Sundar, R. Yimnirun, J. Su, and Q. M. Zhang,J. Phys. Chem. B , 10141 (1997). R. E. Newnham, Properties of Materials: Anisotropy, Symmetry,Structure (Oxford, New York, 2005). W. Liu and X. Ren, Phys. Rev. Lett. , 257602 (2009). F. Li, L. Jin, and R. Guo, Appl. Phys. Lett. , 232903 (2014). R. D. King-Smith and D. Vanderbilt, Phys. Rev. B , 1651(1993). R. Resta, Rev. Mod. Phys. , 899 (1994). N. Sai, K. M. Rabe, and D. Vanderbilt, Phys. Rev. B , 104108(2002). M. Stengel, N. A. Spaldin, and D. Vanderbilt, Nat. Phys. , 304(2009). I. Souza, J. iguez, and D. Vanderbilt, Phys. Rev. Lett. , 117602(2002). P. Umari and A. Pasquarello, Phys. Rev. Lett. , 157602 (2002). H. Fu and L. Bellaiche, Phys. Rev. Lett. , 057601 (2003). M. Veithen, X. Gonze, and P. Ghosez, Phys. Rev. Lett. , 187401(2004). A. Antons, J. B. Neaton, K. M. Rabe, and D. Vanderbilt, Phys.Rev. B , 024102 (2005). M. Stengel and N. A. Spaldin, Phys. Rev. B , 205121 (2007). O. Dieguez and D. Vanderbilt, Phys. Rev. Lett. , 056401 (2006). M. Stengel, D. Vanderbilt, and N. A. Spaldin, Phys. Rev. B ,224110 (2009). X. Wu, M. Stengel, K. M. Rabe, and D. Vanderbilt, Phys. Rev.Lett. , 087601 (2008). X. Wu, K. M. Rabe, and D. Vanderbilt, Phys. Rev. B , 020104(2011). M. Stengel, C. J. Fennie, and P. Ghosez, Phys. Rev. B , 094112(2012). M. Stengel and D. Vanderbilt, Phys. Rev. B , 241103 (2009). M. Stengel, Phys. Rev. Lett. , 136803 (2011). C. Cancellieri, D. Fontaine, S. Gariglio, N. Reyren, A. D. Cav-iglia, A. Fte, S. J. Leake, S. A. Pauli, P. R. Willmott, M. Stengel, P.Ghosez, and J. M. Triscone, Phys. Rev. Lett. , 056102 (2011). J. Hong and D. Vanderbilt, Phys. Rev. B , 115107 (2011). X. Gonze, J. M. Beuken, R. Caracas, F. Detraux, M. Fuchs, G.M. Rignanese, L. Sindic, M. Verstraete, G. Zerah, F. Jollet, M.Torrent, A. Roy, M. Mikami, P. Ghosez, J. Y. Raty, and D. C.Allan, Comp. Mat. Sci. , 478 (2002). J. P. Perdew and Y. Wang, Phys. Rev. B , 13244 (1992). P. E. Blchl, Phys. Rev. B , 17953 (1994). K. F. Garrity, J. W. Bennett, K. M. Rabe, and D. Vanderbilt, Com-put. Mater. Sci. , 446 (2014). H. J. Monkhorst and J. D. Pack, Phys. Rev. B , 5188 (1976). X. Wu, D. Vanderbilt, and D. R. Hamann, Phys. Rev. B ,035105 (2005). J. J. Wang, F. Y. Meng, X. Q. Ma, M. X. Xu, and L. Q. Chen, J.Appl. Phys. , 034107 (2010). F. Li, L. Jin, Z. Xu, and S. Zhang, Appl. Phys. Rev. , 011103(2014). We note that if constant volume is assumed in the process of mea-suring electrostriction coefficients, it can be shown that, for thetetragonal symmetry, Q (cid:39) − Q (1 − Q ) . We did not usethis constraint in structural relaxations since such an assumptionis not suitable for all situations. For instance, for thin films grownon a substrate, the volume can indeed change with d applied in theout-of-plane direction. D. Berlincourt and H. Jaffe, Phys. Rev. , 143 (1958). P. S. Dobal, A. Dixit, R. S. Katiyar, Z. Yu, R. Guo, and A. S.Bhalla, J. Appl. Phys. , 8085 (2001). A. M. Azad, S. Subramaniam, and T. W. Dung, J. Alloys Compd. , 118 (2002). A. R. Akbarzadeh, I. Kornev, C. Malibert, L. Bellaiche, and J. M.Kiat, Phys. Rev. B , 205104 (2005). K. Uchino, L. E. Cross, R. E. Newnham, and S. Nomura, PhaseTransitions , 333 (1980). A. I. Lebedev and I. A. Sluchinskaya, J. Adv. Dielect. , 1550019(2015). R. D. King-Smith and D. Vanderbilt, Phys. Rev. B , 5828(1994). W. Zhong, D. Vanderbilt, and K. M. Rabe, Phys. Rev. B , 6301(1995). N. Zhang, H. Yokota, A. M. Glazer, and P. A. Thomas, Acta Cryst.B , 386 (2011). T. Nishimatsu, M. Iwamoto, Y. Kawazoe and U. V. Waghmare,Phys. Rev. B , 134106 (2010). D. Berlincourt and H. Jaffe, Phys. Rev. , 143 (1958). Y. Liu, G. Xu, C. Song, Z. Ren, G. Han, Y. Zheng, Mater. Sci.Eng., A , 269 (2008). R. Terki, H. Feraoun, G. Bertrand, and H. Aourag, Phys. Stat. Sol.(b) , 1054 (2005). Z. Li, M. Grimsditch, C. M. Foster, S. K. Chan, J. Phys. Chem.Solids , 1433 (1996). S. Yamanaka, H. Fujikane, T. Ha maguchi, H. Muta, T. Oyama, T.Matsuda , S. Kobayashi, and K. Kurosaki, J. Alloys Compd. ,109 (2003). A. G. Kalinichev, J. D. Bass, C. S. Zha, P. D. Han and D. A. Payne,J. Appl. Phys. , 6603 (1993). H. H. Wieder, Phys. Rev. , 1161 (1955). M. Fechner, S. Ostanin, I. Mertig, Phys. Rev. B , 094112(2008). W. Kleemann, F. J. Schfer, M. D. Fontana, Phys. Rev. B , 1148(1984). V. Gavrilyachenko, Sov. Phys. Solid State , 1203 (1970). C. Z. Wang, R. Yu, H. Krakauer, Ferroelectrics194