Electrostriction coefficient of ferroelectric materials from ab initio computation

Electrostriction is an important material property that characterizes how strain changes with the development of polarization inside a material. We show that \textit{ab initio }techniques developed in recent years can be exploited to compute and understand electrostriction of ferroelectric materials. Here, electrostriction coefficients of ferroelectric BaTiO$_{3}$, PbTiO$_{3}$, as well as dielectric BaZrO$_{3}$, are obtained and analyzed. Possible causes of the difference between experimental and numerical results are discussed. We also identified that relative displacements between certain ions at a given polarization could be a good indicator of a material's electrostriction property.


I. INTRODUCTION
Electrostriction effect describes how a material deforms as a polarization develops inside it. Being one of the nonlinear phenomena particularly important in ferroelectric materials, it is the physical basis for their application in actuators, sensors, 1,2 phase shifters, and sonar projectors. 3 The effectiveness of electrostriction is often gauged by the electrostrictive coefficients Q i j k l (i = x, y, z represents Cartesian coordinate) that characterizes how the strain of a material changes quadratically with its polarization, i.e., x i j = Q i j k l P k P l , 4,5 where x i j represents strain tensors (in the following we will use Voigt notation for strain tensor, i.e., x i represents strain with i = 1 . . . 6). Electrostriction is a four-rank polar tensor, which can be observed in essentially all crystals regardless of their symmetries.
Development of materials with large electrostriction is an important task in material science research. For instance, it was suggested that the large piezoelectric coefficient of the lead-free ceramic system BZT-BCT 6,7 is originated from its large electrostriction coefficients. It is obvious that accurately obtaining electrostriction coefficient is a key step in such research. Experimentally, there are several ways to obtain electrostriction coefficient Q: 5 (i) Direct approaches, such as measuring the change in strain as a function of polarization, or (ii) Indirect approaches, such as measuring the change of permittivity (via change of capacitance) as a function of mechanical stress. Computationally, the electrostriction coefficient can be obtained by applying an electric field to the target material, and compute the change of polarization and strain at each step as the electric field increases.

065122-2
Jiang et al. AIP Advances 6, 065122 (2016) established. 8,9 There are two main obstacles: (i) The application of an electric field in a periodic structure needs to be treated carefully; 10,11 (ii) Electrostriction is often obtained by assuming the target material is in high symmetry phases, either because it is above the phase transition temperature, T C (above which many perovskite ferroelectrics are in the paraelectric phase), or because polycrystals or ceramics are used, in which the overall phase is isosymmetric. Such a situation causes the second difficulty in ab initio calculation because the mimicked high symmetry phase is unstable for ferroelectric materials (e.g., the paraelectric cubic phase for BaTiO 3 is not the ground state). Such an unstable system can hardly be used in ab initio computation as it quickly jumps to other more stable phases (e.g., the tetragonal phase in some perovskites) in structural relaxations, forbidding us to set the system at a given polarization and track how the strain changes.
Recently, these two obstacles have been overcome due to the progress in ab initio computation methods, which makes the calculation of electrostriction coefficients easier. The first breakthrough came when Souza et al. and Umari et al. 12,13 implemented the code for applying electric field in periodic structures. With this advancement, it is possible to compute in general how a material responds to applied electric field. As a matter of fact, such method of fixing the electric field E has been applied widely in the study of dielectric, piezoelectric, and ferroelectric behavior of materials. [14][15][16][17] To address the second problem of the unstable high-symmetry phase, Dieguez et al. 18 developed a method that is able to impose a constant polarization, P (fixed-P method). In this way, one can increase the polarization of the system gradually and track the change of strain, which paves the way for numerically obtaining electrostriction coefficient. However, further developments show that a better way is to impose the constraint of a constant electric displacement, D, (fixed-D method). 11 Among other reasons, this approach addresses a few subtle points that are not fully appreciated in previous studies. For instance, (i) The fixed-D method corresponds to imposing an open-circuit electrical boundary conditions, which is especially useful for studying ferroelectric capacitors, 19 superlattices 20-22 and metal-oxide interfaces; [23][24][25] (ii) This approach made it possible to mimic the unstable cubic phase of ferroelectric materials, and extract useful parameters (e.g., electrostriction coefficient). More importantly, with the implementation of such a method in the ABINIT software package, 26 it enables more researchers to obtain electrostriction coefficients via ab initio computation.
Here we use the fixed-D method to study electrostriction in prototypical ferroelectric materials BaTiO 3 (BTO) and PbTiO 3 (PTO). As a comparison, we also present results from BaZrO 3 (BZO), which has the cubic phase as its ground state. Our numerical results indeed shows that this approach can obtain electrostriction coefficients consistent with experimental results, therefore providing a valuable predictive power in discovering materials with better electrostrictive properties, and an insight to the factors that determine electrostriction coefficients.

II. METHOD
The ab initio calculations were performed using the open-source ABINIT software package. 27 Our calculations were performed within density-functional theory in the local-density approximation (LDA) 28 and the projector-augmented-wave (PAW) method 29 is used along with pseudopotentials implemented in the GBRV package. 30 Ba 5s 5p 5d 6s, Ti 3s 3p 4s 3d, Pb 5d 6s 6p, Zr 4s 4p 4d 5s and O 2s 2p orbitals were treated as valence orbitals. A plane wave basis with kinetic energy cutoff of 40 Hartree was used to ensure the convergence in all the calculations. k-point sampling of 6 × 6 × 6 Monkhorst-Pack grid 31 was used for all the 5-atom unit cell. The atomic coordinates of the five-atom unit cell were relaxed until all atomic-force components were smaller than 10 −5 Hatree/Bohr, and the cell size-and-shape was varied until all stress components were below 10 −7 Hatree/Bohr 3 .

III. RESULTS AND DISCUSSION
The fixed-D calculation implemented by Stengel et al. 11,26 is exploited in this work, where we gradually increase the reduced electric displacement field d along the [001] direction. This method enables the relaxation of both atom positions and cell size-and-shape of a given crystal at a given d, the flux of the electric displacement field, defined as d = a 1 × a 2 D/(4π), 11 where a 1 and a 2 are the in-plane lattice vectors, and D is the electric displacement field. At each step of the calculation process, the internal energy of the system, U, the polarization of the system, as well as the relaxed cell (determined by a 1 , a 2 , and a 3 ) are computed. After the computation, 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 relaxed cubic phase a 0 as the reference. The system is initially set to be cubic (Pm3m), but with the application of the fixed-displacement-field D along the [001] direction, it becomes tetragonal (P4mm). The electric field, polarization, internal energy, ionic positions and lattice parameters are all varied at a given value of d. To obtain electrostriction coefficients of each system (BTO, PTO and BZO), a series of computations are performed at given electric displacement, d, along the z direction. Internal energy, polarization and strain are then obtained in each process, and related to the constrained variable, d. We therefore obtain the out-of-plane strain, x 3 , and the in-plane strain x 1 as functions of d, which are then fitted to calculate electrostriction coeffcients. Moreover, in this process we also track ion displacements versus d, which provides some insight to the electrostriction effects. We note that to define response properties it is important to specify the boundary condition or constraint used in their definition. 32,33 Here our calculation of electrostriction coeffcients is under the constraint of zero stress.

A. Ferroelectric materials BaTiO 3 and PbTiO 3
We first show the results of internal energy U, strain in the  Fig. 1(a), we plot the internal energy as a function of d. Clearly, there is an energy minimum at d ≃ 0.24 e, indicating that the cubic BTO is unstable and such instability will introduce a ground state with nonzero d. In Fig. 1(b), we plot the polarization versus d, and find that the polarization changes linearly with d. The reason is that in the region of interest, it is shown that the ε 0 E (ε 0 is the vacuum permittivity) is much smaller than P. Therefore P dominates the relation D = ε 0 E + P, ensuring this linear relation. 26 Figures 1(c) and 1(d) show that the strain versus polarization curves are quadratic, which implies that electrostriction, not piezoelectricity, dominates the mechanical response to polarization. To quantify the electrostriction effect, we employ the expression x i j = Q i j k l P k P l to fit the curves in Figs. 1(c) and 1(d). 5 More precisely, we use the equation x = a 0 + a 2 P 2 for fitting. It is found that Q 11 = 0.137 m 4 /C 2 , which is quite close to the experimental value of Q 11 = 0.11 m 4 /C 2 , and Q 12 = −0.026 m 4 /C 2 , which is comparable to the experimental value of Q 12 = −0.045 m 4 /C 2 . 34,35 The piezoelectric voltage coefficients g 33 and g 31 can be obtained from the relation g i j k = 2Q i j k l P s . 4 For BTO, we obtained g 33 = 0.071 V m/N and g 31 = −0.014 V m/N, which are comparable to the experimental values of g 33 = 0.058 V m/N and g 31 = −0.023 V m/N. 36 This comparison supports the idea that electrostriction may be used to estimate the electromechanical properties (e.g., piezoelectric coefficients) of ferroelectric materials with spontaneous polarization.
Numerical results for PTO is also shown in Fig. 1, where we plot various properties of PTO under different d. Figure 1(e) shows that PTO also has an energy minimum, which agrees well with previous studies 11,26 in terms of the double-well potential and the position of the energy minimum. Figure 1(e) indicates that the depth of energy well of PTO (∼ 60 meV) is much deeper than that of BTO (∼ 5 meV) and the energy minimum corresponds to a larger polarization (P = 0.86 C/m 2 versus P = 0.25 C/m 2 in BTO.). In Figs. 1(g) and 1(h), we show the out-of-plane and in-plane strain versus polarization, and found Q 11 = 0.065 m 4 /C 2 and Q 12 = −0.012 m 4 /C 2 by fitting the two quadratic curves. Comparable experimental values of Q 11 = 0.09 m 4 /C 2 and Q 12 = −0.03 m 4 /C 2 have been reported. 34 We note that the fitting range and whether including higher order terms can

B. Dielectric material BaZrO 3
Having examined ferroelectric materials, we now move to BZO, which is a dielectric material with perovskite structure that has larger lattice constant, smaller thermal expansion, and smaller dielectric permittivity than BaTiO 3 . 37-39 Experimental and theoretical results indicate that BZO remains paraelectric at low temperatures, 39 which means it is not piezoelectric, but can have electrostriction effect. Figure 2 shows the internal energy U and polarization versus d, and strains versus polarization. Figure 2(a) shows that, unlike BTO and PTO, the internal energy U always increases with d. This is consistent with the the fact that cubic BZO has no ferroelectric instability. The electrostriction coefficients we obtained are Q 11 = 0.0944 m 4 /C 2 and Q 12 = −0.0254 m 4 /C 2 . Since direct experimental value on Q 11 and Q 12 are not available, we also use Q h = Q 11 + 2Q 12 to calculate the hydrostatic electrostrictive coefficients and find Q h = 0.0436 m 4 /C 2 . As a comparison, the experimental value of Q h is found to be 0.023 m 4 /C 2 in Ref. 40  and the value we find here may be related to the facts that (i) We did not use the constant volume constraint, 35 and/or (ii) There are anti-phase motions of oxygen octahedra found in BZO. 39,41 It is interesting to note that the numerical Q 11 of BZO obtained here is larger than PTO (but smaller than BTO). However, it is worth noting that since BZO is not ferroelectric, it may be difficult to convert the large electrostriction effect into piezoelectricity, unlike BTO or PTO in which spontaneous ferroelectricity automatically break the centrosymmetry without external electric field. Since BZO does not have a ferroelectric instability, it is also possible to obtain its electrostriction coefficient by directly applying an electric field. Here, we apply an electric field, E, along the [001] direction. We applied electric field up to E = 1.14 × 10 9 V/m. At every step of the applied electric field, the structure of BZO, including its lattice constant, atomic positions, and cell shape, are optimized, similar to what is done in the fixed-D approach. Once the optimized structure of BZO at a given E is obtained, its polarization is calculated using the Berry phase approach. 8 In this way, the electrostriction coefficients are found to be Q 11 = 0.0946 m 4 /C 2 and Q 12 = −0.0265 m 4 /C 2 , which agree well with the results obtained via the fixed-D results.
We also obtained the displacement of each ions versus polarization in Fig. 3(a), where we show the reduced coordinates of BTO versus the polarization. A prominent feature of this figure is that all the atoms shift linearly with polarization. As expected, Ba and Ti ions (being positively charged) move along the direction of the reduced electric displacement field d, while O ⊥ and O ∥ ions (note the Ti-O ∥ and Ti-O ⊥ bonds are 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. We further note that the displacements of Ti, Ba, O ⊥ , O ∥ obtained here with the fixed-D method are compatible with that obtained from force-constant matrix, where the eigenvector is given by 42 Here we find from the slope of the lines shown in Fig. 3(a) that the eigenvector is (0.20, 0.66, −0.32, −0.57). The closeness of these two results and the linear correlation shown in Fig. 3(a) further validate the accuracy of the effective Hamiltonian approach in the study of BTO. 42,43 It is worth noting that the eigenvector shown in Ref. 42 is obtained by diagonalizing the force-constant matrix. With a series of similar fixed-D calculations, we believe the static field-induced ion displacements agree better with the force-constant matrix eigenvectors than the dynamical matrix eigenvector. Figure 3(b) shows the displacement of each ions in PTO as a function of the polarization. Similar to BTO, the displacements of the ions are linearly dependent on P, and consistent with the eigenvector obtained in frozen phonon calculation 42 and experiments. 44 However, one can see that in PTO, the A-site ion Pb moves in a larger distance than the B-site ion, Ti, which is opposite to what happens in BTO. Since the displacement of Pb ion is about twice as large as that of Ti, this result also indicates that in PTO Pb has larger contribution to polarization than Ti. Since Pb is outside oxygen octahedrons, it likely has more room to move without distorting the unit cell along z , unlike Ti ions. In other words, to achieve the same polarization, Ti in PTO has a smaller displacement than BTO, and thus a smaller deformation, which results in a smaller Q 11 . Such a difference may explain why Q 11 of PTO is smaller than that of BTO. We thus suggest that B-site driven ferroelectrics may in general have better electrostriction than A-site driven ones. Figure 3(b) also shows an anomalous change around P = 0.8 C/m 2 . The inverse capacitance of PTO was also found to have such a change, which is discussed in Ref. 11. The ion displacements in BZO versus d also has a linear behavior 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 ions move in the opposite direction. It is worth noting that, in BZO, the A-site atom, Ba, shifts more than the B-site atom, Zr, which is opposite to what happens in BTO, but similar to PTO. Comparison of Fig. 3(c) and Fig. 3(b) to Fig. 3(a) also reveals another important difference between BTO and BZO/PTO, i.e., ξ O ∥ > ξ O ⊥ in BTO while the opposite is true in BZO/PTO. The larger ξ O ⊥ in BZO/PTO is likely due to the larger displacement of the A-site atom (Ba in BZO and Pb in PTO) than that of the B-site atom (Zr in BZO and Ti in PTO) as the displacement A-site ion may cause neighboring oxygen octahedrons to tilt. The importance of such difference on electrostriction becomes apparent when the bond length of the B-site atom and O ∥ along with Q 11 is summarized in Fig. 4 for five different perovskites (Table I summarizes the numerical results of Q 11 , Q 12 , g 33 , g 31 , C 11 , C 12 , C 44 and spontaneous polarization P, the results of elastic constants and spontaneous polarization are in general agree with existing computational and experimental results. 33,42,[45][46][47][48][49][50][51][52][53][54][55][56], which indicates that Q 11 has positive correlation with the B-O ∥ bond length. In other words, longer B-O ∥ means larger theoretical electrostriction coefficients when other things (e.g. polarization) being equal. This relation can be understood by realizing the fact that the B-O ∥ bond length characterizes and, to some extent, determines the unit cell deformation along z, which in turn decides

C. Discussion
With the aid of ab initio computation in general and the fixed-D method in particular, we now have an additional perspective to understand electrostriction. As we know, electrostriction represents a relation between displacements of atoms (that determines the polarization) and structural parameters (that determines strain). It is often taken for granted that the polarization is caused by an external electric field. However, it has been shown that the relation between P and E is rather complex 11,26 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 at a given configuration. The reason is not hard to understand -it is the internal electric field arising from electric dipoles that drives the ferroelectric instability, which in turn results in the spontaneous polarization. 8,43 For ferroelectric materials, it is also found from ab initio computation that for P values of our interest, the required electric field is much smaller than P, i.e., ε 0 E ≪ P, 26 which explains why the relation x i j = Q i j k l P k P l is adopted more than x i j = M i j k l E k E l to characterize electrostriction in ferroelectric materials.
As we have shown, the fixed-D approach is effective to obtain electrostriction coefficients via ab initio computations. This computational approach is obviously useful for discovering materials or structures with desired electrostrictive properties. In order to employ this method properly, it is important to understand potential causes of difference between experimental work and ab initio calculation. First, in many experiments on electrostriction, polycrystals or ceramics are used, which has two implications: (i) The effect of domains, domain walls, and/or grain boundaries in real materials cannot be reflected in ab initio calculation with a rather small number of unit 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 may be mandatory. Second, in comparison with experimental values, one needs to be aware of some electric boundary conditions between different ab initio computations. For instance, the fixed-E computation corresponds to a closed circuit electric boundary condition with a biased voltage and assumes that E is continuous through the system, while the fixed-D method corresponds to an open circuit electric boundary condition and continuity of D. [11][12][13] This subtlety warrants careful consideration before adopting one or the other approach. For complex non-uniform structures, we suggest to use the fixed-D approach, which is particularly important for layered structures, such as superlattices or metal/insulator interfaces, as pointed in Ref. 11. This practice should be considered even for non-ferroelectric materials where applying electric field is not a problem. Third, more complex structures (e.g. the anti-phase oxygen octahedron tiltings in BZO) could affect the accuracy of ab initio results. However, we note that introducing such complex crystal structures will substantially slow down the calculation with the fixed-D approach. Finally, ab initio computations assume zero temperature, therefore the optimized lattice constant may not corresponds to the experimental values at finite temperatures. Experimentally verified values (e.g., lattice constant) may be needed to reduce differences in numerical values.

IV. CONCLUSION
In summary, we have shown that the fixed-D approach can be applied to obtain electrostriction coefficients of ferroelectric materials. By investigating three materials (BTO, PTO, and BZO), we find the fixed-D approach provide a powerful method to compute and understand electrostrictive effects of various materials, and is potentially useful for future predictions of materials or designs. With the aid of this approach, we have also shown that, given a perovskite material, the bond length of the B-site atom and the O ∥ ion is an important indicator of the magnitude of its electrostriction, which could be a useful criterion in screening a large number of perovskites to find those with good electrostrictive effects. In addition, possible difference between electrostriction coefficients obtained theoretically and experimentally are discussed in detail. We thus hope that the approach we demonstrat here will be widely used in discovering and exploring materials with better electrostrictive effects.