The intrinsic piezoelectricity of monoclinic Pb(Zr_1−xTi_x)O₃

Abstract

The phase diagram of Pb(Zr_1−xTi_x)O₃ (PZT) is quite clear; however, the existence and intrinsic piezoelectricity of low-symmetry structures near the MPB (morphotropic phase boundary) remain subjects of considerable uncertainty. The impact of the MPB on the piezoelectric properties of PZT was investigated via first principles study. Different PZT structures were established by virtual crystal approximation (VCA). The density functional theory (DFT) with plane-wave basis set method was utilized to calculate the energy, elastic, and piezoelectric properties. The presence of the monoclinic phase in MPB was confirmed near x ∼ 0.53. This phase exhibited distinct piezoelectric properties. Although the coefficients d₃₃ and d₁₅ were lower, the maximum value of d₃₁ reached −198.33 pm V⁻¹. By exhaustively exploring all possible structures at the MPB, we demonstrate that the intrinsic piezoelectricity of monoclinic structures is not the main contributor to the MPB effect. Further study shows that this is related to the hybridization between the O2p orbital and the d orbital of the central atom. Finally, analysis of the energy barrier along the polarization rotation paths implies a flattening of free energy in the monoclinic phase, suggesting a high intrinsic piezoelectric coefficient. The intrinsic piezoelectric properties of the monoclinic phase may bring new ideas to the study of the MPB effect.

---

1. Introduction

Lead-based piezoelectric materials have garnered significant attention owing to their exceptional ferroelectric properties.^1–7 PZT, as a perovskite-type ferroelectric material, exhibits different phases near Curie temperature (T_c).⁸ Above T_c, PZT crystals exist in a cubic phase. Below T_c, PZT consists of a continuous solid solution of PbTiO₃ and PbZrO₃. At Ti concentration greater than 0.53, the microstructure primarily consists of the tetragonal phase. Conversely, at Ti concentration of less than 0.53, the microstructure consists of the rhombohedral phase. Currently, the most widely used PZT material is located at the MPB.^1,2,6,9,10

There is no consensus on the MPB effect, although there are two main conjectures regarding this effect. The first theory is “field-induced phase transition”. Shrout et al. found that applying a high electric field in the [001] direction induces a transformation from the rhombohedral to the tetragonal phase, accompanied with a large strain in relaxor single crystals.¹¹ Another theory is polarization rotation, first proposed by Fu et al.¹² They compared the energy of different polarization rotation paths in BaTiO₃. Noheda et al. proved the polarization rotation path along a–g–d–e in Pb(ZnNb)O₃–PbTiO₃ single crystals.¹³ The flattening of the free energy surface was considered a major factor contributing to high piezoelectricity. At present, polarization rotation theory is recognized by a large number of researchers.^14,15 However, the intricacies of the MPB exceed conventional understanding. Even pure PbTiO₃, which is usually considered unrelated to the MPB, demonstrated an MPB containing two monoclinic structures under a certain pressure in theory.¹⁶

In recent years, beyond the study of domain wall structures,^17–21 researchers have focused on the contribution of intrinsic piezoelectricity of the low symmetry phase, especially the monoclinic phase.^22–30 Bai et al. proved that the shear strain of a monoclinic structure has an important influence on its ferroelectric properties through thermodynamics.³¹ Damjanovic emphasized the crucial role played by intrinsic piezoelectricity in the MPB.³² The monoclinic phase, serving as a structural entity along the phase transition path, constitutes a significant subject for investigating lead-based piezoelectric crystals through the polarization rotation theory. In 1999, Noheda et al. confirmed that the monoclinic PZT brings unique ferroelectric mechanisms.³³ Cohen et al. confirmed the existence of the R–M_A–T transition utilizing DFT.¹² In 2016, the diffraction peak of a certain Cm structure in PZT was found via SXRD.³⁴ Neutron diffraction also confirmed the existence of low-temperature Cc structure and high-temperature Cm structure in PZT.^35,36 However, the academic community still lacks a complete understanding of the stability and intrinsic piezoelectric properties of various monoclinic phases. Liu et al. believed that compared with other structures, the M_A structure is the key factor in generating the giant piezoelectric response of a perovskite crystal.^37,38 On the phase boundaries of the two monoclinic structures, the larger displacement of Pb in the M_A structure in the external field or strain might be more advantageous in piezoelectric response.³⁹ Jin et al. proposed the ferroelectric adaptive theory, which suggests that the monoclinic phase is composed of rhombohedral phase or tetragonal phase nano-twins, which appear as monoclinic phase as a whole, and such diffraction spots have been observed in Pb(MnNb)O₃–PbTiO₃ ceramics and Pb(ZnNb)O₃–PbTiO₃ single crystals. The emergence of these studies raised new doubts about the intrinsic piezoelectric properties of the monoclinic phase.^40,41 However, evidence was found that the pure Cm structure existed on the side of the tetragonal phase, not the adaptive phase in the PZT, by convergent beam electron diffraction in 2011.⁴

Additionally, DFT is considered to be an effective method to study piezoelectric effects.^42,43 Bellaiche et al. used the DFT method to study the properties of PZT in a certain temperature range at MPB and proposed the possibility of a low-temperature monoclinic phase as a polarization rotation mesophase in 2000.⁴⁴ Liu et al. determined the position of MPB of PZT. Furthermore, the piezoelectric and elastic properties can be calculated.⁴⁵ Duerloo et al. predicted four novel high piezoelectric materials using this DFT method.⁴⁶ Ghosez et al. evaluated the differences in the crystal dynamics of the BaTiO₃, PbTiO₃, and PbZrO₃ lattices via the calculation of phonon scattering.⁴⁷

Finally, monoclinic structure as the transition state of ferroelectric phase transition is a crucial part of polarization rotation theory. The intrinsic piezoelectricity of the monoclinic phase is a significant aspect for researchers to understand the MPB effect. Therefore, in this study, we determined the existence and intrinsic piezoelectricity of monoclinic phases by first principles calculations in detail.

2. Computational methods

This study employed the open-source ABINIT software. For the purpose of investigating the intrinsic piezoelectric effect in different 2 × 2 × 2 VCA lattices, the exchange-correlation potential for electrons was described using GGA implemented in the PBE scheme.⁴⁸ The Monkhorst–Pack method was employed for sampling the K point in the Brillouin zone, with a grid density of 6 × 6 × 6. The BFGS algorithm was employed for geometry optimization during the calculations. The cutoff energy was set to 600 eV, and the maximum displacement was 0.0005 Å. To enhance accuracy, the convergence for the maximum tolerance on the plane wave function squared residual was 10⁻¹⁸. The maximum force tolerance was 0.01 eV Å⁻¹. The electric polarizations were calculated using the Berry phase method.⁴⁹

3. Results and discussion

Fig. 1 illustrates the total energy differences between different phases to P4mm. Table S1† shows the lattice parameters of these structures and Fig. S1† shows the total energy of this phase. The energy difference is defined via the formula E_D(x) = E_P4mm(x) − E_phase(x). Firstly, in the Zr-rich region (x < 0.53), the R3̄m phase exhibits the lowest total energy compared to other phases, with a maximum difference of around 0.04 eV. In the Ti-rich region (x > 0.53), the P4mm phase shows the lowest total energy, which indicates stability in this range. Secondly, high Ti content is not conducive to the existence of a rhombohedral phase; when x > 0.8, the total energy of R3̄m and R3̄c structure will fail to converge. Furthermore, at lower Ti concentration, the R3̄m phase exhibits lower total energy compared to the R3̄c phase. Similarly, this phenomenon is also reflected in the three monoclinic structures; in the range of x < 0.2 and x > 0.7, their structures also fail to converge. It is noteworthy that the Cm(M_B) phase has higher total energy at any Ti concentration, except for a local extremum near x = 0.51, which closely aligns with MPB (x ∼ 0.53). Additionally, the total energy and convergence are not the only verification methods for the existence. In view of this, the study attempts to calculate the free energy to study the existence of Cm(M_B) structures. For other monoclinic phases (M_C, M_A), they exhibit an excessive total energy at the MPB, with a significant energy barrier to even the nearest stable M_B structures (0.012 eV, 0.023 eV).

Fig. 1 The difference in total energy between P4mm and other phases. The spectral line closest to the top in the figure will have the lowest total energy to obtain a stable structure.

By calculating the phonon scattering situation at the Γ point, the relationship between energy and temperature can be determined, expressed as:

The zero-point energy is:

While the free energy F can be expressed as:

At 0 K, the free energy is equal to the zero-point energy.

Based on the calculations of the free energy in Fig. 2, the P4mm phase has the lowest free energy for x > 0.512, while the R3̄m phase has the lowest free energy for x < 0.475.⁵⁰ The range in which the free energy of both takes its minimum is very close to the range commonly considered for the MPB of the PZT. Within a narrow range of 0.475 < x < 0.512, the Cm(M_B) phase exhibits the lowest free energy, indicating its stable existence in this chemical composition. These two methods (total energy and free energy) yield very close results. In order to further explore the origin of this stability of Cm(M_B), the paper first selects the structure of x = 0.5 as the research object. Subsequently, phonon scattering spectra are calculated via the linear response method.

Fig. 2 The free energy near the MPB. The detailed scope is shown in the accompanying figure.

In Fig. 3, firstly, it is found that the degeneracy of the phonon structure of the Cm(M_B) phase is destroyed compared with the literature. The first reason for this result is that the Cm(M_B) phase belongs to the lower crystal group, and the symmetry of the structure is reduced. Secondly, it is due to the doping effect of Zr atoms in the crystal. The softest frequency is 2.72 cm⁻¹, and the second softest frequency is 2.9 cm⁻¹. They are all located at Γ point, which belongs to a vibration of the optical transverse mode of the central atom of the lattice. Finally, since the lowest frequency vibration mode of Cm(M_B) phase is above the virtual frequency region, it can be considered that Cm(M_B) phase can exist stably in theory.

Fig. 3 The phonon spectrum of the Cm(M_B) phase Pb(Zr_0.5Ti_0.5)O₃. The abscissa is the high symmetry point located in the first Brillouin zone (Γ, R, F, L, H, K, M, and A). Their local coordinates are (0, 0, 0), (0, ζ, 0), (ζ, ζ, 0), (ζ, 0, 0), (ζ, ζ, ζ), (0, 0, ζ), (0, ζ, ζ), (ζ, 0, ζ), respectively.

This stability of the monoclinic structure is reflected in its ferroelectric property. For the Cm(M_B) phase located at x = 0.5, the component of the polarization intensity vector is P_x = 0.133, P_y = 0.125, and P_z = 0.0457, respectively. This indicates that Cm(M_B) phase may have different displacement polarization characteristics and polarization rotation characteristics from other structures. This phenomenon should first be reflected in the chemical bond changes in the structure. For PZT, the covalent bond between A position atoms and O with the covalent bond between the central atoms and O2p orbit are equally important contributions to the stability of structure.^51–53 The Born effective charge in the range of 0.45 < x < 0.55 and the displacement of the atom with respect to the symmetric position are shown in Table S2.† The Z* stands for the average Born effective charge, and the , . These two very close values indicate that the two atoms have very close chemical environments.⁴⁷ Additionally, both and are greater than its valency of +4; this suggests that the chemical bond between Ti, Zr, and O is mainly covalent (Z_O1x = −2.494, Z_O1y = −4.429, Z_O1z = −2.209, Z_O3x = −2.381, Z_O3y = −2.341, Z_O3z = −4.684, Z_O4x = −2.75, Z_O4y = −4.846, Z_O4z = −1.966, Z_O6x = −2.106, Z_O6y = −2.041, Z_O6z = −4.058). Meanwhile, the , which may indicate the weakening of the covalent bond between Pb6s and O2p orbit compared with the P4mm structure in this chemical compound. Finally, the most stable structure around x = 0.5 is the R3̄m phase,⁵⁴ which has a similar band gap to the Cm(M_B) structure shown in Fig. S5.† This suggests that the Cm(M_B) structure may have an electronic structure with a reduced DOS of conduction bond and an increased DOS of valence band.

Fig. 4 shows the structure of covalent bonds by calculating the PDOS method. It can be observed firstly that both phases exhibit large band gaps, with E_gCm(MB) = 2.987 eV, and E_gPm3m = 2.463 eV, respectively. Fig. S7† shows the electronic structure of the M_B-Cm phase. Additionally, the PDOS spectra of the two phases exhibit similar peak shapes and identical positions. The peak position of the Cm(M_B) phase is closer to lower energy (Cm(M_B): −5.5 eV to 0 eV, Pm3m: −5 eV to 0 eV) indeed due to the reduction of the total DOS visible in the conduction band. Fig. S8† provides some electronic structures of PZT in different phases. Unlike phase transitions in general materials, the ferroelectric phase transition in the PZT cells is entirely determined by the distortion of oxygen octahedra in the lattice. As can be seen in Fig. S8,† the band structures and DOS between different phases are very similar. Therefore, despite the inherent lack of rigor in this approach, it is possible to determine band shifts by assessing the band gaps of these structures. This difference is mainly the contribution of Ti and Zr atoms in the valence band. This can be attributed to atomic displacements caused by the phase transition process. In PDOS, the range from −5 eV to 0 eV primarily consists of hybridization peaks between the oxygen and the central atoms. The main contributors to the p-orbit DOS are O2p orbit, while the main contributors to the d-orbit DOS are Ti3d and Zr4d orbits. Above the Fermi level, the O2p with Ti3d and Zr4d orbits exhibit similar peak shapes and peak positions. Meanwhile, even though in the conduction band, Cm(M_B) has lower DOS than Pm3m phase, Cm(M_B) crystal still occupies higher energy positions (Cm(M_B): 3 eV, Pm3m: 2 eV); this effect plays a significant role in the stability of the Cm(M_B) phase. This indicates a weak hybridization between the O2p and d orbits of central atoms, suggesting a weakening bonding energy in the Ti–O and Zr–O bonds. At −8 eV to −6 eV, the DOS of Pb6s orbit is drastically reduced compared to Pm3m. Table S3† shows the –COHP of the different PZT phases. Due to the limitations of the VCA, these calculations are based on 2 × 4 × 4 supercells, where the atomic distribution belongs to the T_d point group, representing a crystal with a Zr/Ti ratio of 1. In the Pm3m phase, the –COHP values for Zr–O and Ti–O are 0.0431 and 0.0427, respectively. Similarly, for commonly used ferroelectric structures such as P4mm and R3m, the –COHP values also show a decrease. In the Cm(M_B) structure, the –COHP is the lowest, with values of 0.0406 and 0.0395 for Ti–O₁ and Ti–O₂, respectively, and 0.0412 and 0.0409 for Zr–O₁ and Zr–O₂, respectively. This indicates a reduction in the overlap of bonding orbits between the central atoms and the oxygen atoms in the Cm(M_B) structure, which is consistent with the conclusion of weakened covalent bonds obtained from the PDOS analysis.

Fig. 4 The density of states (DOS) for the Pm3m phase (above) and Cm(M_B) (bottom) phases located at x = 0.5.

The above conclusions based on PDOS are consistent with the conclusions drawn by the study of Born effective charge in the previous paragraphs of the paper. In the electron density image presented in Fig. S2† the monoclinic structure was selected as the study object. In the 2 × 2 × 2 VCA supercell constructed in this paper, all parameters are fully relaxed. By examining the electron density images of the central atom plane and the O3, O4 atom plane in the (001) direction of the crystal, the charge density transfer between the central atom and the O3, O4 atoms are found. This effect may relate to the contraction or expansion of the oxygen octahedron structure in this direction. In addition, the antiferroelectric distortion belonging to the M₂⁺ mode can be observed in the fully relaxed supercell.^55–58 This is preliminarily believed to be related to the monoclinic structure reducing its own energy and improving the stability.

Based on the calculated ferroelectric and elastic properties of the crystal, the piezoelectric constants can be obtained from the rate of change of polarization intensity with respect to strain,^49,59,60 where the strain is set to not exceed 1%. Firstly, when the external electric field is zero, the piezoelectric stress constants can be expressed as:

In practice, this relationship is usually expressed in two parts:In these equations, Ω represents the unit cell volume, R and G denote real-space and reciprocal-space vectors, respectively, η stands for macroscopic strain, and the summation range α extends over three dimensions, u is the intracell strain or called the intrinsic parameter, and the range of k includes all the atoms. Under these conditions, polarization is contributed by the following two components:

p_i = p^s_i + e_ijη_j (6)

where p^s represents the material's spontaneous polarization. Thus, for a monoclinic structure material, the change in polarization along the three directions can be expressed as follows:

For the monoclinic structure in m space group, the elastic constant matrix is:

The relationship between them is as follows:

e_ij = d_ikc_kj (9)

This relationship can be reduced to:

Taking the polarization intensity along the [001] direction as an example, it can be expressed as:

Δp₃ = e₃₁η₁ + e₃₂η₂ + e₃₃η₃ + e₃₅η₅ (11)

By applying strain along the c-axis, the piezoelectric constant can be obtained as:

The piezoelectric strain constant d_ij is related to this formula. The bold e is the meta charge. Z means Born effective charge.

d₃₃ ≈ e₃₃s₃₃ (13)

where s_ij represents the calculated elastic compliance constant.

Fig. 5 presents the relationship between piezoelectric strain constants and Ti concentration. The maximum value of the d₃₃ = 345.71 pm V⁻¹ is observed at x = 0.5. The maximum value of the d₁₅ = 466.84 pm V⁻¹ is found at x = 0.52. Similarly, the maximum value of the d₃₁ = −198.33 pm V⁻¹ is at x = 0.5. These positions are in close proximity to MPB. Compared to experimental values, the Cm(M_B) phase exhibits lower values.^1,2 However, although the highest performance of the structure is a shear piezoelectric coefficient, the relatively high d₃₁ in the Cm(M_B) phase suggests that the structure possesses distinctive transverse piezoelectricity. This is because the above three extreme values are not located in the same chemical composition. Due to the differing symmetries of various ferroelectric structures, utilizing a single metric such as the c/a ratio to measure lattice distortion is not universally applicable. In Fig. S9,† the volume and lattice parameters decrease with increasing x, while the displacement of the central atom relative to the crystal center remains essentially unchanged. This indicates that as x increases, the lattice distortion caused by the central atom intensifies, which is consistent with the behavior of PZT near MPB.

Fig. 5 The piezoelectric strain constants and elastic properties in the Cm(M_B) phase.

A comparison of the piezoelectric properties of several monoclinic structures at the MPB is presented in Table 1. Among these monoclinic structures, the M_B structure exhibits the highest piezoelectric performance; however, this performance still shows a significant gap compared to the rhombohedral and tetragonal structures. For instance, the most successful commercial PZT exhibits a d₃₃ of 750 pm V⁻¹ at the MPB, while the common P4mm and R3̄m phases of PZT have d₃₃ values of 260 pm V⁻¹ and 225 pm V⁻¹, respectively. By exhaustively exploring all possible structures at the MPB based on the VCA method, we demonstrate that the intrinsic piezoelectricity of monoclinic structures is not the main contributor to the MPB effect.

Table 1 Comparison of the piezoelectric properties of each structure at the MPB

pm V^−1	P4mm	R3̄m	R3̄c	Cm(MA)	Cm(MB)	Pm(MC)	Cc
d33	454.76	473.46	442.33	281.17	345.71	212.36	171.34
d15	583.67	615.14	501.82	301.12	466.84	251.37	204.71

---

Fig. S3† illustrates the total energy, the zero point energy, and the frequency of the softest phonon vibration mode at Γ point of the structure around x = 0.5. It is shown that x = 0.5 and x = 0.52 are both located at the bottom of the energy curve, while the softest frequency of x = 0.5 is lower and reaches the range of virtual frequencies. Fig. S4† shows the above properties for x = 0.52 versus c/a. These data suggest that the structure with x = 0.52 is more stable indeed, while the high piezoelectric coefficient d₃₃ and d₃₁ at x = 0.5 may be related to the smaller virtual frequency here, due to the softening of the transverse mode vibration of the central atoms. Fig. 5 also depicts the bulk modulus B and shear modulus G. These parameters are derived from the elastic stiffness constants c_ij and the elastic compliance constants s_ij. The computed results show the maximum and minimum values at x = 0.49, x = 0.53, and x = 0.51, respectively, for B and G. These extrema positions do not align with the position in which it has the direction of maximum polarization, indicating that such elastic properties may play a significant role in influencing the piezoelectric performance. Additionally, these extreme values stand vary close (s_31,x=0.5 = 31.28, s_31,x=0.52 = 29.53, s_33,x=0.5 = 34.86, s_33,x=0.52 = 34.24, s_15,x=0.5 = 22.02, s_15,x=0.52 = 17.99, unit: 10⁻¹² m² n⁻¹). The difference in piezoelectric properties is almost due to differences in electronic structure (e_31,x=0.5 = −6.34, e_31,x=0.52 = −6.23, e_33,x=0.5 = 9.92, e_33,x=0.52 = 9.88, e_15,x=0.5 = 19.41, e_15,x=0.52 = 25.94, unit: C m²). Based on this conclusion, the polarization rotation theory is adopted in this paper. By applying strain to the model, the position of the central atoms is changed, and the rotation of the polarization direction of the structure to the d position is simulated.

Fig. S6† illustrates the total energy, the zero point energy, and the frequency of the softest phonon vibration mode at Γ point of the structure located at x = 0.52. The energy of barriers along the polarization rotation path is defined via the formula E_B(θ) = E_{total energy}(θ) − E_{total energy}(0). θ is the angle of polarization rotation. It can be observed from Fig. 6 that the barrier height gradually increases with the increase in the rotation angle. Moreover, within the lower rotation angle, relatively low barrier heights are observed, corresponding to the flattening of free energy. This phenomenon is also found in the zero point energy of Fig. S6.† This serves as the origin of the intrinsic piezoelectricity in the Cm(M_B) phase.

Fig. 6 Portrayal of the energy barriers along the polarization rotation path of Cm(M_B) phase along the d direction at x = 0.52.

4. Conclusion

In this study, the intrinsic piezoelectric properties of monoclinic PZT were investigated via first-principles study. Firstly, crystal structures versus Ti concentrations were constructed based on VCA, and the position of the MPB was determined at x ∼ 0.53. The free energy indicates the relatively lower of Cm(M_B) phase at 0.475 < x < 0.512. The study of the phonon spectrum shows that the structure has the softest mode indicating stable existence. Secondly, studies on the electronic structure and Born effective charge show that this stability is inseparable from the bond energy enhancement between the O and central atoms. Subsequently, the elastic and piezoelectric properties of the crystal were calculated. The results reveal relatively lower values of the d₃₃ and d₁₅ piezoelectric strain constants in the Cm(M_B) phase, being 345.71 pm V⁻¹ (x = 0.5) and 466.84 pm V⁻¹ (x = 0.52), respectively. However, the d₃₁ component of the piezoelectric strain constant was found to be relatively high, at −198.33 pm V⁻¹ (x = 0.5). By exhaustively exploring all possible structures at the MPB, we demonstrate that the intrinsic piezoelectricity of monoclinic structures is not the main contributor to the MPB effect. Finally, the origin of this piezoelectricity was explained by calculating the energy barriers along the polarization rotation path. Due to the flattening of the free energy, the Cm(M_B) phase of PZT crystals exhibits relatively high intrinsic piezoelectricity. This may bring new ideas to the study of the MPB effect.

Data availability

All authors in this paper promise that all data do not exist in any published articles and all data are free from any restrictions involving privacy and secrecy.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

This research work was supported by calculation technology and computing resource from Shenzhen Key Laboratory of Advanced Thin Films and Applications, and Xi 'an Technological University Development Fund, No. 0103-207000026.

References

1. J.-M. Cha, S.-J. Ha, J.-h. Kim, H.-A. Cha, J.-J. Choi, B.-D. Hahn, K.-H. Cho, S.-Y. Yoon and C.-W. Ahn, Ceram. Int., 2023, 49, 264–270.
2. B. Jaffe, R. S. Roth and S. Marzullo, J. Appl. Phys., 1954, 25, 809–810.
3. S. Kim, W.-J. Lee, Y.-H. Cho, M. Shim and S. Kim, Jpn. J. Appl. Phys., 2013, 52, 091101.
4. R. Schierholz and H. Fuess, Phys. Rev. B:Condens. Matter Mater. Phys., 2011, 84, 064122.
5. N. Zhang, H. Yokota, A. M. Glazer, Z. Ren, D. A. Keen, D. S. Keeble, P. A. Thomas and Z. G. Ye, Nat. Commun., 2014, 5, 5231.
6. M. J. Haun, E. Furman, S. J. Jang and L. E. Cross, Ferroelectrics, 1989, 99, 13–25.
7. X.-K. Wei, S. Prokhorenko, B.-X. Wang, Z. Liu, Y.-J. Xie, Y. Nahas, C.-L. Jia, R. E. Dunin-Borkowski, J. Mayer, L. Bellaiche and Z.-G. Ye, Nat. Commun., 2021, 12, 5322.
8. A. M. Glazer, P. A. Thomas, K. Z. Baba-Kishi, G. K. H. Pang and C. W. Tai, Phys. Rev. B:Condens. Matter Mater. Phys., 2004, 70, 184123.
9. A. M. George, J. Íñiguez and L. Bellaiche, Phys. Rev. Lett., 2003, 91, 045504.
10. R. Guo, E. Alberta, A. Thomas, B. A. Jones and L. E. Cross, Ferroelectrics, 2002, 270, 265–270.
11. S.-E. Park and T. R. Shrout, J. Appl. Phys., 1997, 82, 1804–1811.
12. H. Fu and R. E. Cohen, Nature, 2000, 403, 281–283.
13. B. Noheda, D. E. Cox, G. Shirane, S. E. Park, L. E. Cross and Z. Zhong, Phys. Rev. Lett., 2001, 86, 3891–3894.
14. K. Datta, R. B. Neder, J. Chen, J. C. Neuefeind and B. Mihailova, Sci. Rep., 2017, 7, 471.
15. X. Ke, S. Yang, Y. Wang, D. Wang, L. Zhao, J. Gao, Y. Wang and X. Ren, Phys. Rev. B, 2021, 103, 085132.
16. M. Ahart, M. Somayazulu, R. E. Cohen, P. Ganesh, P. Dera, H.-k. Mao, R. J. Hemley, Y. Ren, P. Liermann and Z. Wu, Nature, 2008, 451, 545–548.
17. Y. Cao, G. Sheng, J. X. Zhang, S. Choudhury, Y. L. Li, C. A. Randall and L. Q. Chen, Appl. Phys. Lett., 2010, 97, 252904.
18. T. Asada and Y. Koyama, Phys. Rev. B:Condens. Matter Mater. Phys., 2007, 75, 214111.
19. Z. Liu, B. Yang, W. Cao and T. Lookman, Phys. Rev. B, 2016, 94, 214117.
20. G. A. Rossetti, A. G. Khachaturyan, G. Akcay and Y. Ni, J. Appl. Phys., 2008, 103, 114113.
21. X.-K. Wei, C.-L. Jia, T. Sluka, B.-X. Wang, Z.-G. Ye and N. Setter, Nat. Commun., 2016, 7, 12385.
22. Y. Gu, F. Xue, S. Lei, T. T. A. Lummen, J. Wang, V. Gopalan and L.-Q. Chen, Phys. Rev. B:Condens. Matter Mater. Phys., 2014, 90, 024104.
23. T. Hyun Kim, J.-H. Ko, S. Kojima, A. A. Bokov, X. Long and Z.-G. Ye, Appl. Phys. Lett., 2012, 100, 082903.
24. A. A. Bokov, X. Long and Z.-G. Ye, Phys. Rev. B:Condens. Matter Mater. Phys., 2010, 81, 172103.
25. I. Lazar, D. Kajewski, A. Majchrowski, A. Soszyński, J. Koperski and K. Roleder, Ferroelectrics, 2016, 500, 67–75.
26. B. Noheda, D. E. Cox, G. Shirane, J. A. Gonzalo, L. E. Cross and S. E. Park, Appl. Phys. Lett., 1999, 74, 2059–2061.
27. B. Noheda, presented in part at the AIP Conference Proceedings, 2000.
28. A. G. A. Nisbet, F. Fabrizi, S. C. Vecchini, M. Stewart, M. G. Cain, T. Hase, P. Finkel, S. Grover, R. Grau-Crespo and S. P. Collins, Phys. Rev. Mater., 2021, 5, L120601.
29. K. Roleder, A. Majchrowski, I. Lazar, R. W. Whatmore, A. M. Glazer, D. Kajewski, J. Koperski and A. Soszyński, Phys. Rev. B, 2022, 105, 144104.
30. G. Catalan, A. Lubk, A. H. G. Vlooswijk, E. Snoeck, C. Magen, A. Janssens, G. Rispens, G. Rijnders, D. H. A. Blank and B. Noheda, Nat. Mater., 2011, 10, 963–967.
31. G. Bai, Y. Han, W. Li and C. Gao, Phys. Rev. B, 2022, 106, 184105.
32. D. Damjanovic, IEEE Trans. Ultrason. Ferroelectrics Freq. Control, 2009, 56, 1574–1585.
33. B. Noheda, J. A. Gonzalo, R. Guo, S.-E. Park, L. E. Cross, D. E. Cox and G. Shirane, AIP Conf. Proc., 2000, 535(1), 304–313.
34. L. Fan, J. Chen, Y. Ren, Z. Pan, L. Zhang and X. Xing, Phys. Rev. Lett., 2016, 116, 027601.
35. G. Fraysse, J. Haines, V. Bornand, J. Rouquette, M. Pintard, P. Papet and S. Hull, Phys. Rev. B:Condens. Matter Mater. Phys., 2008, 77, 064109.
36. A. K. Singh, D. Pandey and O. Zaharko, Phys. Rev. B:Condens. Matter Mater. Phys., 2006, 74, 024101.
37. H. Liu, J. Chen, L. Fan, Y. Ren, Z. Pan, K. V. Lalitha, J. Rödel and X. Xing, Phys. Rev. Lett., 2017, 119, 017601.
38. N. Zhang, H. Yokota, A. M. Glazer, D. A. Keen, S. Gorfman, P. A. Thomas, W. Ren and Z.-G. Ye, IUCrJ, 2018, 5, 73–81.
39. N. Zhang, H. Yokota, A. M. Glazer, Z. Ren, D. A. Keen, D. S. Keeble, P. A. Thomas and Z. G. Ye, Nat. Commun., 2014, 5, 5231.
40. Y. M. Jin, Y. U. Wang, A. G. Khachaturyan, J. F. Li and D. Viehland, J. Appl. Phys., 2003, 94, 3629–3640.
41. Y. M. Jin, Y. U. Wang, A. G. Khachaturyan, J. F. Li and D. Viehland, Phys. Rev. Lett., 2003, 91, 197601.
42. O. Hudak, Phase Transitions, 2008, 81, 1073–1079.
43. B. Völker, C. M. Landis and M. Kamlah, Smart Mater. Struct., 2012, 21, 035025.
44. L. Bellaiche, A. García and D. Vanderbilt, Phys. Rev. Lett., 2000, 84, 5427–5430.
45. Y. Y. Liu, L. Yang and J. Y. Li, J. Appl. Phys., 2013, 113, 183524.
46. K.-A. N. Duerloo, M. T. Ong and E. J. Reed, J. Phys. Chem. Lett., 2012, 3, 2871–2876.
47. P. Ghosez, E. Cockayne, U. V. Waghmare and K. M. Rabe, Phys. Rev. B:Condens. Matter Mater. Phys., 1999, 60, 836–843.
48. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 3865.
49. L. Dong, J. Lou and V. B. Shenoy, ACS Nano, 2017, 11, 8242–8248.
50. S. Baroni, S. de Gironcoli, A. Dal Corso and P. Giannozzi, Rev. Mod. Phys., 2001, 73(2), 515–562.
51. R. E. Cohen, Nature, 1992, 358, 136–138.
52. J. Shi, I. Grinberg, X. Wang and A. M. Rappe, Phys. Rev. B:Condens. Matter Mater. Phys., 2014, 89, 094105.
53. R. Tharamas, J. Jutimoosik, A. Bootchanont, P. Kidkhunthod, S. Rujirawat, Z.-G. Ye and R. Yimnirun, Integr. Ferroelectr., 2017, 177, 137–142.
54. S.-Y. Liu, Q.-S. Shao, D.-S. Yu, Y.-K. Lü, D.-J. Li, Y. Li and M.-S. Cao, Chin. Phys. B, 2013, 22, 017702.
55. B. K. Mani, S. Lisenkov and I. Ponomareva, Phys. Rev. B:Condens. Matter Mater. Phys., 2015, 91, 134112.
56. J. S. Baker and D. R. Bowler, Phys. Rev. B, 2019, 100, 140401(R).
57. N. Sicron, B. Ravel, Y. Yacoby, E. A. Stern, F. Dogan and J. J. Rehr, Phys. Rev. B:Condens. Matter Mater. Phys., 1994, 50, 13168–13180.
58. A. K. Tagantsev, K. Vaideeswaran, S. B. Vakhrushev, A. V. Filimonov, R. G. Burkovsky, A. Shaganov, D. Andronikova, A. I. Rudskoy, A. Q. R. Baron, H. Uchiyama, D. Chernyshov, A. Bosak, Z. Ujma, K. Roleder, A. Majchrowski, J. H. Ko and N. Setter, Nat. Commun., 2013, 4, 2229.
59. G. Sághi-Szabó, R. E. Cohen and H. Krakauer, Phys. Rev. Lett., 1998, 80, 4321.
60. Z. Wu and H. Krakauer, Phys. Rev. B:Condens. Matter Mater. Phys., 2003, 68, 014112.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d4ra06270a

---