Suppressing The Ferroelectric Switching Barrier in Hybrid Improper Ferroelectrics

The A₃B₂O₇ compounds considered in this study are the $n=2$ members of the Ruddlesden-Popper series Ruddlesden1957 ; Ruddlesden1958 . They can be considered as layered perovskites with an extra AO layer inserted after every 2 perovskite bi-layers (i.e. 4 atomic layers) along the [001] direction (Fig. 1a). The extra AO layers cause a shift by ($a/2$, $a/2$, 0) on the $ab$ plane, and hence the structure becomes body centered tetragonal with space group $I4/mmm$ (#139). This shift also breaks the connectivity of the oxygen octahedra, and the AO double layer is held together by mostly ionic bonds between the A-site cations and O anions. The resulting dimensional reduction has important consequences on the electronic structure and lattice response (For example, Ref.’s Birol2011 ; Zhang2013Iridate ; Wang2013 ; Li2019Janotti ). Apart from the dimensional effects, the different periodicity of the Ruddlesden-Popper phases along the layering direction ($c$ axis, or the [001] direction) leads to a smaller Brillouin zone than ABO₃ perovskites. The equivalents of various structural instabilities that are at different points of the Brillouin zone in the ABO₃ perovskites can fold back onto the same point in A₃B₂O₇ Ruddlesden-Poppers, which leads to interesting couplings between them as discussed below. (This point can be qualitatively understood in analogy to a subduction problem, where a zone boundary mode of the parent group corresponds to a zone center mode of the subgroup. For example, when the unitcell of a cubic perovskite is doubled along the [001] axis as a result of cation order, the spacegroup becomes $P4/mmm$ and the zone boundary $X_{5}^{-}$ mode splits into $\Gamma_{5}^{-}\oplus X_{2}^{-}\oplus X_{3}^{-}$, where $\Gamma_{5}^{-}$ is polar. While there is no direct group-subgroup relationship between the Ruddlesden-Popper and perovskite structures, the $n=2$ Ruddlesden-Poppers have 2 perovskite blocks in their unitcells, and it is thus possible to recognize some phonon modes folded onto the $k_{z}=0$ plane.)

By far the most common structural distortions that decrease the symmetry of oxide perovskites is the oxygen octahedral rotations: About 90% of all oxide perovskites have this type of distortion in their crystal structures, which reduces the symmetry of the parent $Pm\bar{3}m$ phase Lufaso2001 . These distortions can be described in terms of symmetry-adapted-modes, which can be classified by irreducible representations (irreps) of the parent spacegroup $Pm\bar{3}m$ miller1967tables . The phonon modes that correspond to these distortions are the M point mode M${}_{2}^{+}$, which is an in-phase rotation of octahedra around one axis, and the R point mode R${}_{5}^{-}$, which is an out-of-phase rotation of octahedra around one axis. The former is denoted by a ‘$+$’ superscript in the Glazer notation, such as a⁰a⁰c⁺, and the latter is denoted by a ‘$-$’ superscript, such as a^-a^-a^-. The most common rotation pattern that more than half of all oxide perovskites have is a^-a^-c⁺, which leads to the space group $Pnma$ (#62) Woodward1997 . Another distortion that is often significant in the $Pnma$ structure is the X${}_{5}^{-}$ out-of-phase A-site displacement. Unlike the M${}_{2}^{+}$ and R${}_{5}^{-}$, the X${}_{5}^{-}$ often does not show up as an unstable phonon mode in the high symmetry ($Pm\bar{3}m$) phase. Rather, it is an improper order parameter, which attains a nonzero magnitude only because of a trilinear coupling in the Landau free energy

$$\mathcal{F}_{\textrm{trilinear}}=\gamma M_{2}^{+}R_{5}^{-}X_{5}^{-}.$$

(1)

The presence of $\mathcal{F}_{\textrm{trilinear}}$ in the free energy expansion, which is imposed by group theory, guarantees a nonzero X${}_{5}^{-}$ distortion whenever the octahedral rotations M${}_{2}^{+}$ and R${}_{5}^{-}$ are present, no matter the sign of the coupling $\gamma$.

Instabilities in the A₃B₂O₇ Ruddlesden-Poppers that are similar to the M${}_{2}^{+}$ and R${}_{5}^{-}$ normal modes in the ABO₃ perovskites give rise to a wider range of different combinations and resultant symmetries. (For simplicity, we follow the convention to refer octahedral rotations around the out-of-plane ($c$) axis as ‘rotations’ (OOR), and the rotations around the in-plane axes as ‘tilts’ (OOT).) One reason for this is that there is a new degree of freedom, since the body-centered primitive cell now contains two oxygen octahedra. Also, the double AO layers break the connectivity of oxygen octahedra, and hence the relative phase of neighboring octahedra on either side of the double layer is not fixed. As an example, we consider the modes relevant to the $A2_{1}am$ phase observed in Ca₃Ti₂O₇ and many other HIF Ruddlesden-Popper compounds in Fig. 2. In ABO₃ perovskites, there are two possible rotation patterns around, for example, the $c$ axis: in-phase (M${}_{2}^{+}$, a⁰a⁰c⁺) or out-of-phase (R${}_{5}^{-}$, a⁰a⁰c^-). In the A₃B₂O₇, on the other hand, there are four possibilities: The X${}_{2}^{+}$ mode corresponds to an in-phase rotation of the two octahedra in one perovskite slab that consists of 5 atomic layers, and is the primitive unit cell. However, X${}_{2}^{+}$ is a two dimensional irrep, and depending on its direction a particular pair of octahedra on either side of a double AO layer can have either in-phase or out-of-phase rotations, as shown on the left two panels of Fig. 2c. Similarly, the rotations that are out-of-phase within one perovskite slab transform as the two dimensional irrep X${}_{1}^{-}$, as shown in the right panels of Fig. 2c.

The most relevant octahedral rotation modes in A₃B₂O₇ all have the same wavevector: they correspond to X point normal modes. This leads to a richer set of possibilities for the modes induced by trilinear couplings compared to ABO₃ perovskites. In the trilinear coupling terms in ABO₃ perovskites, an M and an R mode has to couple with an X mode due to the translational symmetry. In A₃B₂O₇ compounds, on the other hand, the trilinear couplings that contain two separate X modes can contain either an M mode or a $\Gamma$ mode as the third mode. (M point is denoted as the Z point in the convention of Ref. Bradley2010 .) The reason is that there are two separate X points on the Brillouin zone that are related to each other via a four-fold rotation, and depending on which pair of X wavevectors are chosen, their sum can either give the $\Gamma$ or the M wavevector. In Table 1, we list the possible trilinear couplings between two X modes and a third mode in the A₃B₂O₇ structure, and in Fig. 3, we display the polarization patterns of some of these structures.

Hybrid improper ferroelectricity in the A₃B₂O₇ compounds emerges due to the trilinear coupling between X${}_{2}^{+}$ and X${}_{3}^{-}$ modes, which induces a polar displacement $\Gamma_{5}^{-}$. In the HIF structure with space group $A2_{1}am$ (#36), each AO layer has a polarization, which are in alternating directions within each perovskite slab, and hence cancel each other - but only partially. As a result, every perovskite slab between the double AO layers have a net dipole moment. These moments order in parallel and give rise to a macroscopic polarization (Fig. 3a). A different combination of the same X modes can couple to the M${}_{5}^{-}$ mode, leading to anti-parallel slab dipoles, and hence to an anti-polar phase shown in Fig. 3b. (We refer to phases with nonzero dipole moments of each perovskite slab as either polar or antipolar.) Other combinations of the X modes couple with different M modes, such as M${}_{5}^{+}$ or M${}_{2}^{+}$, and give rise to nonpolar phases, where the dipole moments of each atomic layer cancel each other within each perovskite slab between to double AO layers (Fig. 3c-d). (We refer to phases where dipole moments of each slab are zero as ‘nonpolar’.) Many of these phases are observed to emerge in various A₃B₂O₇ oxides under biaxial strain or equivalent doping, and are also shown to be important as intermediate states in the coherent switching of polarization Lu2016 ; Li2018 ; Nowadnick2016 ; Munro2018 ; Yoshida2018 ; Li2017 . This is in addition to single-tilt systems observed, for example, at finite temperature Pitcher2015 . In the next subsection, we draw the strain–tolerance factor phase diagram of these compounds to identify regions where these antipolar and nonpolar multi-tilt phases emerge.

Most –more than half– of oxide perovskites have a tolerance factor of $\tau<1$, and attain the space group $Pnma$ at low temperatures Lufaso2001 . The corresponding octahedral rotation pattern a^-a^-c⁺ is also common in A₃B₂O₇ Ruddlesden-Poppers, and gives rise to the polar space group $A2_{1}am$ observed in HIFs. In addition to the polar phase, strain phase diagrams of these compounds often abound with transitions to nonpolar phases introduced in the preceding subsection. As an example, in Fig. 4a, we present the energy of three lowest energy structures for Sr₃Sn₂O₇ as a function of biaxial strain Lu2017 . The zero temperature DFT calculations reproduce the experimentally observed room temperature phase $A2_{1}am$ in the unstrained compound. Both tensile and compressive strain decrease the energy difference between this phase and the next lowest energy state, and there are phase transitions to nonpolar phases for strain $\gtrsim 2.5\%$ on either direction. Similar strain driven transitions have been predicted for Sr₃Zr₂O₇ and Ca₃Ti₂O₇ HIF compounds previously, and the pattern of octahedral rotations often change under strain in the ABO₃ compounds as well. A common trend in A²⁺B⁴⁺O₃ perovskites is that tensile biaxial strain suppresses OOR around the out-of-plane axis, whereas compressive strain enhances it. Sr₃Sn₂O₇ follows a similar trend: The transition under tensile strain is to the $P4_{2}/mnm$ phase, which has only X${}_{3}^{-}$ tilts, whereas the transition under compressive strain is to the $Aeaa$ phase, which has only the X${}_{1}^{-}$ rotations around the $c$ axis. The transition to these nonpolar phases is not a result of a continuous suppression of polarization by strain: the magnitudes of polarization in the $A2_{1}am$ phase on both phase boundaries are sizable, and is even enhanced under tensile strain, as shown in Fig. 4(b).

In order to elucidate the behavior of different HIF compounds under strain, in Fig. 5 we map out the strain – tolerance factor phase diagram by considering 11 different A₃B₂O₇ compounds. (We do not include 2 compounds with larger tolerance factors, since they do not display any OOR or OOT. Most of these compounds have been studied from first principles before, but to the best of our knowledge, this is the first time that this information is compiled to display all compounds together. We consider a strain range of $\pm 4\%$, which covers the experimentally feasible range. For most of the compounds with $\tau<1$ that we consider, the lowest energy unstrained structure is $A2_{1}am$, which corresponds to the HIF phase. For $0.92\lesssim\tau\lesssim 1$, nonpolar structures emerge under both tensile and compressive strain. We observe three different nonpolar structures: $Pnab$ and $P4_{2}/mnm$ under tensile strain, and $Pnab$ and $Aeaa$ under compressive strain. They correspond to the following changes in the octahedral rotations and tilts:

• •

  Compressive strain induced OOT suppression (leads to $Aeaa$): This is observed in Sr₃Sn₂O₇ and Ca₃Ge₂O₇. The OOT mode amplitude drops to zero and OOR mode phase changes under compressive strain as shown in Fig. 6(c),(d).

• •

  Tensile strain induced OOR suppression (leads to $P4_{2}/mnm$): This is observed in Sr₃Zr₂O₇ and Sr₃Sn₂O₇. Similar to the first situation, but the OOR mode drops to zero under tensile strain instead of OOT mode, as shown in Fig. 6(b-c).

• •

  Tensile/compressive strain induced OOR phase change (leads to $Pnab$): This is observed in Ca₃Ti₂O₇ under both tensile and compressive strain, in Sr₃Zr₂O₇ under compressive strain, or in Cd₃Ti₂O₇ under small tensile as well as compressive strains. (Fig. 6(a-b)). Amplitudes of both the OOR and OOT mode retain non-zero, but the in-phase OOR mode changes into out-of-phase manner. This structure is shown in figure  3(b). The A-site cations around two interfaces move in the opposite direction, which cancels the polarization in bulk.

Some of these transitions are explained by local measures such as the global instability index (GII), which is known to predict the octahedral rotation patterns and angles in ABO₃ perovskites successfully Lufaso2001 ; Salinas-Sanchez1992 . For example, the transition to $P4_{2}/mnm$ in Sr₃Sn₂O₇ is coincident with the strain value above which the GII of this phase is the smallest Supplement . However, the GII by itself does not explain why the polar $A2_{1}am$ structure is preferred over the $Aeaa$ one, for these two phases have very similar GII values under compressive strain. It is possible that the interplay of GII with the long-range Coulomb interaction (which is an important factor in stabilizing the polarization in proper ferroelectrics such as BaTiO₃ghosez1996coulomb ) is responsible of the transition to the $Aeaa$ phase.

The transition to a single-tilt system can be explained phenomenologically by the cross term between OOR and OOT – a large OOR might suppress OOT and vice versa. All compounds in the $A2_{1}am$ follow the same aforementioned trend as many ABO₃ perovskites that compressive strain enhances OOR, whereas tensile strain enhances OOT (Fig. 6a-d). (For example, see Ref’s Lu2016 ; Yang2012 ; Zayak2006 .) This trend is likely the result of the strain reducing particular B–O bond lengths, which can be increased by the OOT or OOR distortions. The lowest order cross term between the OOR and OOT in the free energy is $F\sim\beta R^{2}T^{2}$ (where we denote the amplitudes of rotations and tilts by $R$ and $T$ respectively). For fixed value of $R$, this term renormalizes the coefficient of the $T^{2}$ term $\sim\alpha T^{2}$ as $\sim(\alpha+\beta R^{2})$, and hence for large OOR $R^{2}>-\alpha/\beta$, the tilting instability is suppressed, and it becomes energetically favorable to have no tilts, as is the case in compressively strained Sr₃Sn₂O₇ in the $Aeaa$ phase.

A phenomenological explanation of the strain induced transition to nonpolar $Pnab$ structure requires not only the biquadratic terms between the OOR and the OOT modes, but also various trilinear terms that couple these modes to other antiferrodistortive displacements Lu2016 . It is particularly interesting that in Ca₃Ti₂O₇, this transition is re-entrant in the sense that it happens under both tensile and compressive strains. The GII does not have an obvious trend that explains this transition Supplement , and the electrostatic interaction between the O ions on different layers is possibly important Supplement . We leave the microscopic explanation of this transition to a future study.

Enhanced susceptibilities near second order phase transitions can be exploited to design materials with large responses, for example, magnetic permeability or dielectric constants. While no such enhancement of linear susceptibility is mandated near first order transitions, it is nevertheless possible to obtain large response near a first order phase boundary if the external field can induce the transition. Examples of demonstrations of this approach include Terfenol, Pb(Zr,Ti)O₃, and BiFeO₃ Zeches2009 ; Newnham1998 ; Birol2012 . The phase boundaries of structural transitions depend on strain very sensitively, and as a result, this approach is a promising means to enhance the response of materials via strain.

The question we focus on in this subsection is whether the ferroelectric polarization switching barrier is affected when strain is used to tune the materials to the vicinity of the polar-nonpolar phase transitions. In order to answer this question, we use the minimum energy barrier for coherent polarization switching as a proxy to the coercive electric field. While in an actual experiment defects, domain structure, as well as size and shape effects significantly alter the coercive field, trends of coherent switching barrier can be used as a first principles proxy to the trends of the coercive field Beckman2009 as explicitly shown in HfO₂ Clima2014 . (Finite element methods which take into account the domain structure provide much lower switching barriers Dittrich2002 .) In practice, the coherent switching field calculated from the first principles energy barrier by assuming that the dipole moment in every unit cell of an infinite crystal switches at the same time is a gross overestimate. As a result, we don’t report the electric field required for switching, but instead report only the energy barriers.

Since in the hybrid improper ferroelectric A₃B₂O₇ compounds the polarization emerges as an improper order parameter through a trilinear coupling with rotation and tilting modes, switching one of these two modes is necessary to switch the polarization. It was recognized as early on as in the first HIF paper that this makes different switching pathways possible, and that the corresponding energy barriers can be tuned by strain Benedek2011 . Later, the work of Nowadnick and Fennie Nowadnick2016 analyzed the possible roles of different switching mechanisms, and Munro et al. used the idea of distortion symmetry groups to identify other switching pathways Vanleeuwen2015 ; Munro2018 . Since then, the energetics of switching in various HIF compounds have been studied, for example in Ref. Liu2019 . However, to the best of our knowledge, a comparison of different compounds and their strain dependence have not been performed yet.

In Fig. 7(a-d), we show four possible polarization switching pathways. We follow the convention of the distortion symmetry groups to name these pathways Vanleeuwen2015 . This process involves identifying not only the symmetry operations shared by all images on the pathway, but also those operations that reverse the distortion, which is the polarization in this case. The latter are referred to as distortion reversal symmetries, and are denoted by a ‘*’ superscript. For example, $Pn^{*}ab$ means that each image along the switching path has two glide planes with translations along a and b axes; and the glide plane $n^{*}$ reverses the distortion. Three of the switching paths we consider ($Pn^{*}ab$, $Pb^{*}nm$ and $Pn^{*}am$) have a similar name as their intermediate phase (up to the asterisks), because the spatial symmetry elements of the intermediate phase either remain unchanged or become reversal symmetry operation for other images. But this is not the case for $Pn^{*}2_{1}^{*}m$. All of the four are so-called 2-step switching pathways, where there exists a local minimum of energy on the switching path, as seen from Fig. 8(a), and they are the lowest ones among such paths for the 3 compounds we considered. They each have distinct intermediate states, but the same initial and final states. Since the Ruddlesden-Popper structure consists of weakly bound perovskite blocks separated by an interface between two rock-salt AO layers, it is possible to consider supercells extended along the [001] direction, and polarization being switched in one perovskite block at a time. This, in principle, gives rise to an infinite number of different switching pathways, the barrier energy per formula unit can be arbitrarily small (since only one block out of arbitrarily many switches at each step.) This has been observed in Ref’s Munro2018 ; Liu2019 , where typically the 4-step switching paths have lower (but comparable) barriers than the 2-step ones, which in turn have lower barriers than the single-step paths. (The path with a very large number of steps can be considered to be a simple model of domain wall in motion along the [001] direction.) However, this does not necessarily imply that the pathway with the highest number of steps determines the coercive field, because what is more important for the switching under an electric field is the slope of the energy vs. polarization curve Beckman2009 . For simplicity, as well as computational manageability, we focus only on 2 step switching pathways.

Each of the four pathways can be reproduced within the same doubled conventional cell as the polar structure. The $Pn^{*}ab$ and $Pn^{*}am$ pathways (Fig. 7(a-b)) involve changes in the direction of the OOR mode, and both of them have nonpolar intermediate structures, with space groups $Pnab$ and $Pnam$ respectively. The out-of-phase displacements of the A-site cations are along the polar axis in both of these intermediate structures. The $Pb^{*}nm$ pathway involves switching the direction of the OOT mode ($X_{3}^{-}$), whereas in the $Pn^{*}2_{1}^{*}m$ both OOR and OOT change directions, as shown in Fig. 7(c-d). Mode decompositions of these switching pathways are given in the supplementary information Supplement .

In Fig. 8(a), we plot the energy as a function of the reaction coordinate for these four switching pathways in unstrained Sr₃Sn₂O₇. The energy barriers are comparable and the lowest one is for the $Pn^{*}am$ pathway. Results presented in Fig. 8(b) show how the energetics of this path behaves under tensile strain: Tensile strain monotonically decreases the $Pn^{*}am$ switching barrier, thus lowering the expected coercive field required for switching. This is not a surprising result, since the OOR’s weaken under tensile strain, as shown in Fig. 6(c) and the $Pn^{*}am$ pathway involves a change in the OOR character. What is interesting, and important for applications, is that this reduction in the switching barrier is not accompanied with a lower polarization under tensile strain (Fig. 4). Thus, strain can be used as a means to lower the coercive field of hybrid improper ferroelectrics.

The strong strain dependence of the switching barrier is not specific only to Sr₃Sn₂O₇, or the $Pn^{*}am$ pathway. In Fig. 9(a-b), we show the barrier for different switching paths of Sr₃Sn₂O₇ and Ca₃Ti₂O₇ as a function of strain throughout the strain range that the HIF phase is stable. While the error bars in the energy barriers from the NEB calculations cause the curves to be rather rugged, two trends are evident: (i) under tensile strain, the barriers for pathways that involve changing the direction of the OOR mode ($Pn^{*}am$ and $Pn^{*}ab$) are lowered, and (ii) under compressive strain, the barrier for the pathway that only involve changing the direction of the OOT mode ($Pb^{*}nm$) is lowered. These are consistent with the tendencies towards OOR and OOT distortions becoming weaker under tensile and compressive strain as discussed earlier. Near 0% strain, the lowest barrier pathway switches from $Pb^{*}nm$ or $Pn^{*}2_{1}^{*}m$ to either $Pn^{*}ab$ or $Pn^{*}am$, and either strain direction leads to a lower coherent switching energy barrier. The lowest barriers are obtained near the phase boundaries between the polar and nonpolar phases, and the maximum suppression is about 50% in both compounds.

Ca₃Sn₂O₇ has a lower tolerance factor than Sr₃Sn₂O₇ and Ca₃Ti₂O₇, and it does not display a strain induced phase transition in the strain range we considered. It does not show a strain induced change in the switching pathway, or a significant decrease in the switching barrier either (Fig. 9(c)). This is likely because this compound is very far from the phase boundaries, and with its small tolerance factor, it has such large OOR and OOT that the strain induced changes in the instabilities are inconsequential.

Density functional theory calculations are performed using the projector augmented wave approach blochl1994projector as implemented in the Vienna Ab-initio Simulation Package (VASP) VASP1 ; VASP2 , and using the PBEsol generalized gradient approximation PBEsol . All calculations are done in a 48-atom (4 formula unit) supercell, which can be viewed as a $\sqrt{2}\times\sqrt{2}\times{2}$  multiple of the primitive cell of the reference $I4/mmm$ structure. A $\Gamma$-centered $6\times 6\times 2$ grid of k-points is used for the Brillouin zone integrals.

We consider all A₃B₂O₇ compounds with A = Ca, Sr, Ba and B = Ti, Zr, Sn, Ge, as well as Cd₃Ti₂O₇ Wang2017 ; Yoshida2018 ; Yoshida2018-2 ; lee2013exploiting ; Li2018 ; Kennedy2011 ; henriques2007ab ; moriwake2011first . These compounds are all band insulators with sizable gaps, so using the PBEsol generalized gradient approximation is expected to reproduce the crystal structures with reasonable accuracy. Biaxial strain boundary conditions are simulated by fixing the in-plane lattice constants, and allowing the out of plane component, as well as internal atomic positions, to relax with an force threshold of 2 $meV$/Å. The zero strain is defined for each compound by the $a$ lattice constant obtained by completely relaxing the structure in the reference high symmetry structure $I4/mmm$.

The Goldschmidt tolerance factor Goldschmidt1926 , which is used as a simple measure to predict tendency towards octahedral rotations, and is originally defined in terms of the ionic radii $r$ using

$$\tau=\frac{r_{A}+r_{O}}{\sqrt{2}\left(r_{B}+r_{O}\right)}$$

(2)

is instead calculated using the bond lengths for 12 coordinated A-site ($d_{AO}$) and 6 coordinated B-site ($d_{BO}$) ions from the bond valence model as

$$\tau=\frac{d_{AO}}{\sqrt{2}\cdot d_{BO}}.$$

(3)

(This approach is following Ref. Lufaso2001 .)

In order to calculate the minimum energy barrier for polarization switching, climbing-image nudged elastic band (CI-NEB) method was used to further relax linearly interpolated switching paths to the minimum energy path NEB . The spring constant was set to $5$ eV/Ų, and a convergence criterion of $1$ meV per supercell was used. Distortion symmetry groups Vanleeuwen2015 ; Padmanabhan2020 are used to enumerate and name the possible initial pathways following Ref. Munro2018 with the help of the DiSPy package Munro2019 . All the switching pathways reported in the text retain their symmetry for all values of the reaction coordinate under NEB calculation.

As various points in this paper, symmetry and group theory related arguments are built using the Isotropy Software Packageisotropy2007 and the Bilbao Crystallographic Serverbilbao1 ; bilbao2 ; bilbao3 . VESTA software was used for visualization of crystal structures.VESTA2008