Polarization twist in perovskite ferrielectrics

Because the functions of polar materials are governed primarily by their polarization response to external stimuli, the majority of studies have focused on controlling polar lattice distortions. In some perovskite oxides, polar distortions coexist with nonpolar tilts and rotations of oxygen octahedra. The interplay between nonpolar and polar instabilities appears to play a crucial role, raising the question of how to design materials by exploiting their coupling. Here, we introduce the concept of ‘polarization twist’, which offers enhanced control over piezoelectric responses in polar materials. Our experimental and theoretical studies provide direct evidence that a ferrielectric perovskite exhibits a large piezoelectric response because of extended polar distortion, accompanied by nonpolar octahedral rotations, as if twisted polarization relaxes under electric fields. The concept underlying the polarization twist opens new possibilities for developing alternative materials in bulk and thin-film forms.

material because the BNT-BT system contains various structural features [21][22][23][24] . In situ X-ray diffraction (XRD) analysis using high-energy synchrotron radiation demonstrates that electric fields induce an extended polar displacement associated with nonpolar octahedral rotations in ferrielectric crystals, as if twisted polarization relaxes and stretches. The twisted polarization creates a large piezoresponse and eventually becomes a ferroelectric polarization. Our simulations based on density functional theory (DFT) and phenomenological theory show that this concept stems from a structural coupling between nonpolar octahedral rotation and polar distortion.

Results
Strain and polarization properties of ferrielectric single crystals. Figure 1 shows the strain (S) and polarization (P) properties under an electric field (E) applied along the [001] direction. The crystals present a butterfly-type S-E curve ( Fig. 1a) with large jumps at E = ± 20 kV/cm and a double-hysteresis-like P-E loop (Fig. 1b). Unlike antiferroelectrics with a remanent polarization (P r ) of zero 25 , the crystals have a certain P r . This apparent P r associated with the pinched hysteresis suggests that the crystals have the ferrielectric phase in space group P4bm at E = 0, as revealed for BNT-7%BT powders by neutron diffraction structural analysis 26 .
In Fig. 1c,d, we display the S-E and P-E curves measured under unipolar electric fields. Prior to the measurements, an E of 100 kV/cm (E poling ) was applied for the poling treatment, and the unipolar E was applied in the same direction as E poling . The E-induced S with increasing E (Fig. 1c) features an extremely large hysteresis with an abrupt jump at E = 17 kV/cm, and the S value reaches 1.0% at E = 100 kV/cm. In the high E region of 60-100 kV/cm, the S-E curve exhibits a linear response. A decrease in E below E = 20 kV/cm diminishes S substantially, decreasing it to 0% at E = 0 kV/cm. We observed a similar hysteresis in the P-E curve (Fig. 1d). These E-induced properties associated with the abrupt jumps in S and P are repeatedly identified. Figure 2 shows the results of the in situ synchrotron radiation XRD (SR-XRD) measurements under electric fields. In the poled state at E = 0 kV/cm (Fig. 2a), all diffraction spots arise from a ferrielectric tetragonal structure in space group P4bm 22,26,27 . Among the fundamental hkl reflections from the pseudocubic cell, we detect 1/2{o o e} superlattice reflections, where o denotes an odd number and e an even number. These superlattice reflections are caused by the tilting of the TiO 6 octahedra of the P4bm phase. The crystals in the virgin state present the same diffraction pattern including the superlattice reflections (Supplementary Note 1). We verify that the crystals belong to the tetragonal P4bm phase in both the virgin and poled states. The P4bm structure is characterized by a small P s along the c axis, which is associated with the octahedral tilting around the c axis; the details are explained in Supplementary Notes 2 and 3. The superlattice reflections exhibit a substantial decrease in intensity under electric fields and disappear at E = 100 kV/cm (Fig. 2b). The structural analyses of the SR-XRD patterns reveal that the crystals under E = 100 kV/cm have a single phase of ferroelectric tetragonal P4mm without octahedral tilting. We note that the application of E induces the phase transition from the tetragonal P4bm (ferrielectric) to the tetragonal P4mm (ferroelectric) phase.

Structural analyses via in situ synchrotron radiation XRD.
When the external field is decreased from 100 kV/cm to 2 kV/cm, the superlattice reflections reappear (Fig. 2c), and then the crystals return to the initial state at E = 0 kV/cm. Figure 2d shows the relative intensity of the 3/2 5/2 1 (superlattice) reflection as a function of E obtained by the in situ measurements. To compensate for the time-dependent change in X-ray strength, the reflection intensity in Fig. 2d is normalized to the background intensity around the 3/2 5/2 1 spot. The right axis indicates the rotation angle of the TiO 6 octahedra, which is approximately proportional to the square root of the 3/2 5/2 1 reflection intensity. Here, we adopt the rotation angle of 3.0° at E = 0 kV/cm determined for BNT-7%BT powders by neutron diffraction structural analysis 26 .
With increasing E, the superlattice intensity decreases markedly and reaches 14 ± 25% at E = 100 kV/cm. With decreasing E, a significant intensity cannot be detected down to 30 kV/cm. The superlattice intensity recovers below E = 20 kV/cm and returns to its initial value (91 ± 24%) at E = 0 kV/cm. We found the following experimentally: the application of an E exceeding 20 kV/cm stabilizes the P4mm phase, and the E-induced P4mm phase returns to the P4bm phase at E = 0 kV/cm. These results clarify that the E-induced P4bm-P4mm phase transition has a reversible switching path, as discussed in detail below.
We found that both the superlattice reflections and the fundamental hkl reflections vary under electric fields. In Fig. 2e, we show the evolution of the 1 3 8 reflection recorded in the same region of the imaging plate as a function of E. The application of E yields neither a split nor a satellite of the fundamental reflections and only leads to a shift of the peak positions. The line profile beside each spot in Fig. 2e displays the cross-sectional intensity in the vertical direction. The 1 3 8 reflection exhibits a sudden shift between E = 10 kV/cm and 20 kV/cm with increasing E, which is attributed to the transition from the P4bm phase to the P4mm phase. The P4mm phase persists when E > 20 kV/cm with increasing and decreasing electric fields. When E is below 25 kV/cm, the spot moves to the line of the P4bm phase and finally returns to its initial position at E = 0 kV/cm. These results show that neither a multidomain state nor a mixed-phase appears in our crystals within the resolution of the SR-XRD measurements. Figure 3a displays the lattice parameters of the pseudocubic unit cell as a function of increasing E. The crystal lattice expands along the c axis (|| E), accompanied by a shrinkage in the a-a plane (⊥ E) because of the converse piezoelectric effect. The P4bm phase has a small tetragonal distortion of c/a = 1.0007 at E = 0 kV/cm, which is consistent with the results obtained for BNT-7%BT powders (c/a = 1.0003) via neutron diffraction 26 . The crystals feature an abrupt elongation along the c axis at E = 17 kV/cm. The P4bm-P4mm phase transition leads to a discontinuous change not only in the lattice parameters but also in the tetragonal distortion. Our structural analyses demonstrate that the crystals exhibit an E-induced reversible phase transition between the P4bm and P4mm phases.  Fig. 1c) measured with increasing E. It is clear that the macroscopic response of S bulk originates from the microscopic deformation of the unit cell, i.e., the change in S unit-cell caused by the application of E. The quantitative agreement between S bulk and S unit-cell provides direct evidence that the large jump in S bulk at E = 17 kV/cm stems from the E-induced phase transition from the P4bm phase (c/a ≤ 1.002) to the P4mm phase (c/a ≥ 1.010).

Discussion
We investigated the P4bm (ferrielectric)-P4mm (ferroelectric) phase transition induced by electric fields using density functional theory (DFT) calculations. The structural optimizations revealed that BNT with a rock-salt-like A-site ordering forms two types of tetragonal structures with comparable energies. The detailed methodology is described in the Calculation Methods section. Figure 4 depicts the crystal structures of the two tetragonal phases. One exhibits a weak-polar structure associated with octahedral tilting around the c axis ( Fig. 4a,b), which is denoted by the 'T' phase' . Because the T' phase has an octahedral tilt involving a small P s along the c axis as schematised in Fig. 4c (that is, it features the P4bm structure), we can regard the T' phase as the P4bm phase. The other has a strong-polar structure exhibiting a large P s along the c axis without octahedral tilting ( Fig. 4d-f), which is called the 'T phase' . Because the T phase has the essential structural elements of the P4mm phase and the symmetry of its TiO 6 octahedra can be approximated in space group P4mm, the T phase corresponds to the P4mm phase.
Here, we discuss the total free energy (G DFT ) obtained by the DFT calculations. The formula for G DFT is expressed as DFT t otal h cell where U total denotes the total energy of the perovskite unit cell, TS is the entropy term, and V cell is the cell volume. We employed zero-temperature DFT calculations; thus, TS = 0. We define two structural parameters, g and m, representing the crystallographic features of the T' and T phases: g is the degree of the polar displacement leading to P s along the c axis, and m is the degree of the octahedral rotation around the c axis. The fractional coordinates of the atoms in the unit cell can be described by the parameters (g, m) and are defined such that the T' and T phases (Fig. 4) are located at (g, m) = (0, 1) and (1, 0), respectively. The detailed definitions of these parameters are provided in the Calculation Methods section. Figure 5 shows the G DFT profiles as functions of g and m. For each profile, the parameter g (the polar displacement) is fixed, and the parameter m (the octahedral rotation) is varied. At g = 0, the G DFT profile shows a double-minimum potential with respect to m. The potential minima correspond to the T' phase exhibiting the octahedral rotation, the degree of which is defined as m = 1. Increasing g from 0 to 1 lowers the potential height between the minima, and the m values at the minima approach zero. At g = 1.0 (that is, for the T phase with the large polar displacement), the G DFT profile presents a single-minimum potential. We note that the g and m parameters exhibit a trade-off relationship. Large values of both g and m cannot coexist; thus, a large g is achieved by sacrificing m and vice versa.
To gain further insight into the relation between the T' and T phases, we examine the free-energy potential based on the phenomenological Landau-Ginzburg-Devonshire (LGD) theory 28 . We start from a free-energy function of G LGD to express a ferrielectric-ferroelectric phase transition by LGD 0  where two order parameters, P and R, are included: P is the net polarization along the polar axis, R denotes the degree of the octahedral rotation around the polar axis, and the rests are the independent parameters. By choosing appropriate parameters (see Supplementary Note 4), we can obtain the G LGD potential involving two local minima at (P, R) = (P a , R a ) and (P b , 0) in the P ≥ 0 and R ≥ 0 region, where P a and P b indicate small and large polarizations, respectively, and R a represents the degree of the octahedral rotation. Therefore, we can consider the former minimum as the T' phase (a small P s with octahedral tilting) and the latter as the T phase (a large P s without octahedral tilting). Next, we introduce the following relations to correlate the order parameters (P and R) of the G LGD function with the structural parameters (g and m):    Figure 6a shows the G LGD landscape obtained by fitting the G LGD function to the G DFT potential. The difference between these values can be smaller than 0.05 eV, especially by virtue of the coupling terms (γ 22 P 2 R 2 + γ 42 P 4 R 2 ) in Eq. 2 (see Supplementary Note 4). In addition, a reasonable relation of the polarizations, 0 < P a < P b , is realized for the T' and T phases. Figure 6b represents the cross-sectional profiles of the G LGD surface at various g values. The overall features of the G DFT profile (Fig. 5) are well reproduced by the G LGD theory. The agreement between the two calculations (Figs 5 and 6a) confirms that the LGD theory with the two order parameters (P and R) can describe the T'-T phase transition in a simple analytical manner.
The calculations based on the LGD theory enable us to identify the structural factor dominating the T'-T transition. Figure 6c,d displays the energy curves with respect to m, corresponding to the octahedral rotation (G R ) and the polarization-rotation coupling (G c ) in the G LGD function: G R = α 2 R 2 /2 + α 4 R 4 /4 and G c = γ 22 P 2 R 2 + γ 42 P 4 R 2 . The sum of the two terms, G R + G c , determines the overall shape of the G LGD landscape. Because G R is independent of g, G c is decisive in the G LGD function. The γ 22 and γ 42 coefficients are positive; therefore, the G c term yields a positive quadratic variation with m. At g = 0 (Fig. 6c), G c plays no role in G LGD because of the small P (= P a ), and hence, the G LGD function at g = 0 directly reflects the double-minimum G R character. At g = 1.0 (Fig. 6d), the large G c quadratic negates the double minimum in G R and changes the G LGD function to the single-minimum potential. The LGD theory demonstrates that the g-dependent G c character governs the free-energy feature in the T'-T phase transition. The G c term (namely, the polarization-rotation coupling), dominates the competing energy relation between the T' and T phases.
We then analyse the detailed structural variation under electric fields using DFT calculations. The influence of E on the electronic energy U total in Eq. 1 can be included by the perturbation expansion after the discretization (PEAD) approach 29 using the macroscopic polarization, as defined in the 'modern theory of polarization' 30,31 . Figure 7 shows the G DFT maps in the (g, m) subspace under an electric field (E [001] ) applied along the [001] direction (the c axis). The red dashed line depicts the energy valley exhibiting the switching path between the T' and T phases. At E = 0 MV/cm (Fig. 7a), the G DFT map has a deep minimum at (g, m) = (0, 1) and a shallow minimum at (g, m) = (1, 0). The former corresponds to the T' phase, which is characterized by a small polarization with apparent octahedral tilting. The latter is equivalent to the T phase, exhibiting a large polarization without octahedral tilting. Figure 7d presents the G DFT variations along the energy valley as a function of g. When an E [001] is applied, the T' phase moves towards the T phase (denoted by I). Moreover, E [001] lowers the energy barrier from the T' to T phase (Δ G T'→T ). The Δ G T'→T value at E [001] = 3 MV/cm (Fig. 7b) is 9.0 meV, which is considerably less than 56 meV at E = 0 MV/cm (Fig. 7a). Given that E [001] reaches 4 MV/cm (Fig. 7c), Δ G T'→T is as small as 1.3 meV, which is one order of magnitude smaller than the thermal energy at room temperature (k B T RT ): 26 meV. The T' phase exceeds this small Δ G T'→T by the thermal energy and eventually transits to the T phase under such high fields (II in Fig. 7d).
With decreasing E [001] , the T phase remains stabilized in the high-E [001] region with a slight decrease in its polar displacement (III in Fig. 7d), and the energy barrier from the T to T' phase (Δ G T→T' ) is lowered. The Δ G T→T' value decreases to 4.5 meV at E = 1 MV/cm, which is much smaller than k B T RT . This change results in a reverse transition from the T to T' phase (IV in Fig. 7d). As E [001] decreases further, the T' phase returns to its initial state (g, m) = (0, 1) at E = 0 (V in Fig. 7d). Our DFT calculations taking electric fields into account show that the switching path for the reversible phase transition lies between the T' and T phases, as identified by the experimental results (Fig. 2). Figure 8a,b displays the polarization (P DFT ) and strain of the unit cell (S DFT ) as a function of E [001] obtained by our G DFT calculations. The S DFT values are estimated from the unit-cell structures at the local minima moving along the energy valley in the (g, m) space. We assume here that the transition between the T' and T phases occurs once the barrier of Δ G T'→T or Δ G T→T' is below 5 meV. The P DFT variation (Fig. 8a) shows that the crystal lattice undergoes a reversible phase transition with an apparent hysteresis. As E [001] increases, the T' phase with a small P DFT changes to the T phase with a large P DFT at a threshold E [001] . As E [001] decreases, the P DFT value remains large and then drops sharply because of the reverse transition from the T to T' phase. The S DFT curve (Fig. 8b) also exhibits hysteresis with abrupt jumps with both increasing and decreasing E [001] . The jumps in S DFT are attributed to the phase transitions between the T' and T phases, as can be also seen in the P DFT variation. The polarization and strain curves simulated by the DFT calculations are in qualitative agreement with the experimental results (Fig. 1).
Here, we discuss the ferrielectric features under unipolar electric fields and the associated piezoelectric response. Our experiments demonstrate that an application of E enhances P (Fig. 1d) and reduces the rotation angle (R) of the TiO 6 octahedra (Fig. 2d). Our DFT calculations also show that an E [001] application increases the polarization (Fig. 8c) and the tetragonality (c/a, Fig. 8d) accompanied by a decrease in R. It is worth noting that the free-energy surface of the T' phase is incredibly flat compared with that of the T phase (see Fig. 7 and Supplementary Note 5). With increasing E, the T' phase moves substantially in the (g, m) subspace along the energy valley with a gentle slope. An external field induces an extended polar displacement with a suppressed rotational distortion of the oxygen octahedra, where the P s vector behaves as if the twisted polarization relaxes and stretches. The polarization twist yields a linear increase in S DFT with E [001] (Fig. 8b) because of the converse piezoelectric effect in the low-field region. We observed this behaviour also in the S unit-cell data from the SR-XRD analyses (Fig. 3b). The piezoelectric strain constant of the ferrielectric P4bm phase estimated from the slope of S unit-cell is as high as 1,000 pm/V, which is, to our knowledge, the highest for lead-free ferroelectrics. The higher strain constant of the T' phase is also validated by the DFT calculations (Fig. 8b). The superior piezoelectric character in the ferrielectric phase stems from the flat free-energy profile unique to the polarization twist compared with the polarization extension in the ferroelectric phase. We expect that ferrielectrics involving octahedral rotation potentially possesses enhanced piezoelectric properties compared with ferroelectrics. Thus, we identify the polarization twist as a new type of piezoresponse, in addition to polarization rotation 11,12 and polarization extension 32 . Exploiting the interplay between nonpolar and polar instabilities in ferrielectrics is expected to provide a new degree of freedom for developing high-performance perovskite materials for use in bulk and thin-film forms.

Methods
Crystal growth and electric measurements. Single BNT-7%BT crystals were grown by a top-seeded solution growth method at a high oxygen pressure (Po 2 = 0.9 MPa) 33,34 . First, BNT-BT powders prepared via a solid-state reaction were mixed with a flux composed of Bi 2 O 3 , NaF and BaCO 3 and placed in a platinum crucible. The mixture was soaked at 1,100 °C for 4 h, slowly cooled to 1,070 °C at a rate of − 2 °C/h, and then cooled to room temperature. The details of the crystal growth are described in refs 34-36. We confirmed the chemical composition of BNT-7%BT via X-ray fluorescence and inductively coupled plasma-atomic emission spectrometry.
The crystals were cut into plates with thicknesses of 0.2 mm. The axis normal to the crystal plates was along the < 100> direction. Gold electrodes were sputtered on both cut surfaces of the platelet crystals. We investigated the electrical properties along the [001] direction. Strain properties were measured using a laser Doppler displacement meter.

SR-XRD measurements.
We performed the SR-XRD analyses with a transmission geometry using a large cylindrical two-dimensional imaging plate camera at BL02B1 in the SPring-8 synchrotron radiation facility 37,38 . We adopted a high SR energy of 35 keV [wavelength: 0.035313 (15) nm] such that the X-ray could penetrate through the crystals and a high-angle diffraction pattern could be observed. The X-ray beam incident on the crystals was 150 μ m in diameter, and the measurement temperature was 25 °C. To investigate the intrinsic response of the unit cell with respect to E, we conducted in situ measurements of the diffraction patterns under the application of a static E parallel to the [001] direction. We obtained four X-ray oscillation photographs with an oscillation angle of 10° at different angles between the (100) face and the incident X-ray beam at each electric field. We refined the lattice parameters via the least-squares method using approximately 150 peaks of the hkl reflections (d ≥ 0.04 nm). We analysed the diffraction data by imposing the following structural constraints on the lattice parameters: a = b and α = β = γ = 90°.
Calculation methods. First-principles calculations based on density functional theory (DFT) were performed to investigate the E-induced P4bm (T')-P4mm (T) phase transition. While BNT and its related materials do not exhibit a long-range ordering of the A-site cations 39,40 , we need to apply a periodic arrangement of Bi and Na to the BNT model structure for the DFT calculations. We chose pure BNT with a rock-salt-like A-site ordering 40,41 as a model substance. First, we constructed the supercell consisting of × × 2 2 2 perovskite, which is defined as the T super cell. The fractional coordinates of the constituent atoms are listed in Table 1. The T super cell contains two types of atomic displacements: polar displacement along the c axis (z i ) yielding P s and nonpolar displacement (x i ) leading to a tilting of the TiO 6 octahedra around the c axis (i is the index of atoms). The T and T' structures can be described by the T super cell as follows: the T phase has a relatively large z i with x i = 0, whereas the T' phase features a small z i with a certain degree of x i . The T phase constructed in the T super cell has I4mm symmetry, and the T' phase in the T super cell has P4 2 nm symmetry. These symmetries are lower than those of the actual P4mm and P4bm phases for the T phase and T' phase, respectively, due to the A-site ordering adopted in the T super cell: i.e., I4mm is a subgroup of P4mm while P4 2 nm is a subgroup of P4bm.
The experimental results presented in Fig. 3 show that the free energy (G) of the T' phase is small but comparable to that of the T phase at E = 0. We adopted a hydrostatic pressure (p h ) of 6.0 GPa in our DFT calculations, at which the G feature that the T' phase is stabilized compared to the T phase (metastable) is satisfied. The methodology is summarized in Supplementary Note 6. The structural optimization of the T super cell at p h = 6.0 GPa leads to the T' phase, in which the fractional coordinates are denoted by z i (T′ E=0 ) and x i (T′ E=0 ). The optimization of the T super cell under the constraint of x i = 0 results in the T phase, and its coordinate is defined as z i (T E = 0 ). We confirmed that the optimized structure of the T phase possesses a larger tetragonal distortion (c/a) than the T' phase (c/a = 1.047 for the T' phase and 1.053 for the T phase).
Here, we introduce the two-dimensional subspace (g, m) representing the crystal structures with various degrees of P s along the c axis and of the TiO 6 tilting around the c axis. The parameter g indicates the degree of the polar displacement, which is expressed as The parameter m expresses the degree of the TiO 6 tilting as Oa 0 In the (g, m) subspace at E = 0, the T phase is present at (g, m) = (1, 0), and the T' phase is located at (g, m) = (0, 1). If the parameters g and m are given, the fractional coordinates of the T super cell can be specified according to Eqs 5 and 6. The structure optimizations of the lattice parameters (a and c) were performed at p h = 6.0 GPa with the fixed fractional coordinates at each of the mesh points in the (g, m) subspace. After the optimizations, these T super cells are used to calculate the free energy G DFT in the (g, m) space according to Eq. 1.
The DFT calculations were performed according to the generalized gradient approximation 42 using a plane wave basis set, as implemented in the Vienna ab initio simulation package (VASP) 43 . We used the projector-augmented wave potentials 44 with the valence-electron configurations of 5d 10 6s 2 6p 3 for Bi, 2p 6 3s 1 for Na, 3p 6 3d 2 4s 2 for Ti, and 2s 2 2p 4 for O. The Perdew-Burke-Ernzerhof functional modified for solids (PBEsol) 45 was employed for the exchange-correlation potential. A plane-wave cut-off energy of 520 eV was adopted, and the total energy was converged to less than 10 −5 eV in all calculations. The Monkhorst-Pack k-mesh of 3 × 3 × 2 was adopted for geometrical optimization calculations of the T super cell. The crystal structures obtained by the DFT calculations (Fig. 4) were drawn using the three-dimensional visualization software VESTA 46 .