Antiferroelectricity in a family of pyroxene-like oxides with rich polymorphism

Antiferroelectrics have potential applications in energy conversion and storage, but are scarce, particularly among oxides that otherwise display rich ferroic behaviours. Are we overlooking potential antiferroelectrics, simply because we have not discovered their corresponding ferroelectric phase yet? Here we report a first-principles study suggesting this is the case of a family ABO$_3$ pyroxene-like materials, characterized by chains of corner-sharing BO$_4$ tetrahedra, a well-known member being KVO$_3$. The irregular tetrahedra have an electric dipole associated to them. In the most stable polymorph, the dipoles display an antipolar pattern with zero net moment. However, upon application of an electric field, half of the tetrahedra rotate, flipping the corresponding dipoles and reaching a ferroelectric state. We discuss the unique possibilities for tuning and optimization these antiferroelectrics offer. We argue that the structural features enabling this antiferroelectric behaviour are also present in other all-important mineral families.

There is currently considerable interest in finding new antiferroelectric (AFE) materials, owing to their technological importance and relative scarcity [1][2][3] . Applications of AFEs rely on their unique response to an applied electric bias, featuring a double hysteresis loop that is the result of a field-induced phase transition to a polar state. This double-loop makes them particularly efficient for energy applications 4,5 , as e.g. in pulsed-power capacitors 6,7 . The field-induced transformation usually results in a large mechanical response, suitable for transducers and actuators 8,9 . Finally, AFEs display an inverse electro-caloric effect [10][11][12] .
To identify new AFEs, one could proceed as follows: Chosen a materials family and a pertinent highsymmetry structure (e.g., the ideal cubic phase of perovskite oxides ABO 3 ), one could use first-principles simulations to find compositions that present similarly-strong polar and antipolar phonon instabilities of this parent phase. Any such compound is likely to present metastable polar and antipolar polymorphs of similar energy. Then, if the antipolar state were more stable, and given that an electric field will always favour the polar one, we would have a good candidate to display AFE behaviour. This is exactly the exercise that led to the discoveries here reported. We ran a high-throughput first-principles study of the dynamical stability of numerous compounds in the cubic perovskite phase, and found that alkali vanadates (NaVO 3 , KVO 3 , RbVO 3 and CsVO 3 ) show dominant and related polar and antipolar soft modes (see Supp. Note 1 and Suppl. Fig. S1). Then, when trying to identify the polar and antipolar polymorphs associated to these instabilities, we found they result in very large distortions of the cubic phase, to the point that the perovskite lattice is nearly destroyed. The resulting energy minima include many low-energy antipolar polymorphs, as well as their polar counterparts. Indeed, according to our calculations, they constitute a very promising family of AFE compounds.
Results.-In order to identify stable structures, we run first-principles molecular dynamics simulations (see Methods) of several alkali vanadates, which reveal the existence of numerous low-energy metastable polymorphs with some common features: a strong tetragonal distortion, with c/a ratios larger than 1.3; a sublattice of A cations that retain a perovskite-like configuration; and a 4-fold coordination of the vanadium atoms yielding corner-sharing VO 4 tetrahedra.
The connection between the perovskite phase and the ground state of these vanadates is sketched in Fig. 1, which shows the result of a nudged elastic band (NEB) calculation (see Methods for details). Starting from the perovskite structure (Fig. 1a), the system undergoes a transition to a phase with a strong tetragonal distortion,  in which a first V-O bond is broken in each AVO 3 unit (Fig. 1b, and Suppl. Video #1). The obtained structures are reminiscent of supertetragonal phases known for perovskites like BiFeO 3 and PbVO 3 , which show a large c/a ratio and consequently a layered structure formed by BO 5 pyramids 13,14 . In a subsequent step, the VO 5 square pyramids rotate along the V-O bond of the apical oxygen (Fig. 1c), while the V cations move towards one of the oxygen atoms at the base of the pyramid (Fig. 1d), ending in the rupture of a second V-O bond. Each of the resulting VO 4 groups is corner-linked to two adjacent tetrahedra, forming chains.
We obtain several polymorphs featuring the mentioned vertex-sharing oxygen tetrahedra, inspired by metastable phases observed in our molecular dynamics simulations. The tetrahedra can form one-dimensional (1D) zigzag (ZZ) structures along the [110] pseudo-cubic direction (Fig. 2a). The chains can alternatively follow a 1D battlement-like (BM) pattern along the [010] axis (Fig. 2b), or even form closed loops yielding zerodimensional ring (RG) structures (Fig. 2c). In the ZZ structure the apical oxygens point in opposite directions for adjacent chains, closely resembling what is found in inosilicate pyroxenes [15][16][17] , while the RG structures have close analogs in the cyclosilicate family [18][19][20] . The ZZ structure is the most stable one, the lowest-energy BM and RG structures lying 117 meV and 125 meV per formula unit (f.u.) above. The ZZ ground state we find for KVO 3 (as well as for RbVO 3 and CsVO 3 ) is orthorhombic with space group P bcm and coincides with the experimentally observed structure 16,[21][22][23][24][25][26] .
Importantly, the V-O bonds forming the backbone of the tetrahedra chains are longer than the remaining two hovering V-O bonds 16,27 -see bond lengths for the ZZ polymorph of KVO 3 in Fig. 2d. This difference of bond lengths induces a local electric dipole which lies along the direction defined by the V cation and the centre of the tetrahedron edge formed by the lingering oxygens, as shown in Fig. 2d.
We now focus on the ZZ phase of KVO 3 since it is the ground state of the best-established compound among those studied. For a given chain, the ZZ structure gives rise to an antipolar pattern of the in-plane component of the dipoles, while the out-of-plane component remains constant. Since the apical oxygens point in opposite directions for contiguous chains, the out-of-plane component changes sign from chain to chain, yielding a striped antipolar pattern (Fig. 2e). An obvious question arises: is KVO 3 antiferroelectric? A positive answer can only be given if a related ferroelectric (FE) phase accessible via an electric field is found.
We crucially realised that a rotation of a tetrahedron about the backbone edge is tantamount to switching the out-of-plane component of its electric dipole. We thus constructed a ZZ structure with FE inter-chain ordering (ZZ-F) that turned out to be stable and 135 meV/f.u. higher in energy than the antipolar ground state (ZZ-A). The polar order of the ZZ-F polymorph is shown schematically in Fig. 3g. We calculate its spontaneous polarisation 28 to be 0.093 C/m 2 (0.054 C/m 2 ) for the out-of-plane (in-plane) component, of the same order of magnitude as perovskite FEs like BaTiO 3 (0.43 C/m 2 ) 29 . Since the V-O distances in the ZZ-A phase differ by less than 0.5 % from those of the ZZ-F phase, we estimate the out-of-plane sublattice polarisation in the ZZ-A phase to be approximately 0.046 C/m 2 , i.e., half of the ZZ-F polarisation value. Using the tools described in Methods, we find that the optimized ZZ-F state has P m symmetry. An intermediate ferrielectric state (ZZ-i), with half of the tetrahedra in one of the chains flipped (Fig. 3f), was also found to be metastable for KVO 3 .
What remains in order to confirm the AFE character of KVO 3 is to find a connecting path for the field-induced transition. To this end we carry out NEB calculations between the ZZ-A and the ZZ-F polymorphs. The re- sults are shown in Fig. 3a (thick black line). A continuous energy path is found with an energy barrier of 147 meV/f.u. The switching of the dipole chain occurs in a step-wise fashion, with half of the tetrahedra rotating at a first stage -and ending up in the ZZ-i phase -, and the remaining ones rotating at a second stage (see Suppl. Video #2). Note that we extend the NEB path up to a FE state in which the polarisation lies fully outof-plane; for KVO 3 such a state, with space group P ma2, is a saddle point of the energy. The two components of the spontaneous polarisation that acquire non-zero values along the switching path are shown in blue solid lines in Fig. 3a. The out-of-plane polarisation increases abruptly in the first stage, and slower in the second one. Along the path a non-zero in-plane polarisation develops, perpendicular to the tetrahedron chains, showing an oscillating behaviour due to the rigid dipoles crossing the plane of the chains in each of the two switching steps. Finally, we estimate the behaviour of KVO 3 under application of an out-of-plane electric field. The response to the applied field can be approximated by constructing an electric enthalpy with the form where U 0 , P and v are, respectively, the energy, polarisation and volume at zero field at each step of the path, as obtained from first principles; E is the applied electric field. The ZZ-A state is chosen as the zero of energy for convenience. By introducing this approximated enthalpy we avoid running costly first-principles calculations explicitly considering an applied field. In Fig. 3a the energy profile of KVO 3 along the switching path is shown for different values of the field. The hysteresis loop of the polarisation under an outof-plane electric field can be numerically reconstructed from the enthalpies and polarisations in Fig. 3a. More specifically, we consider the case of T = 300 K, assuming a thermal energy of 26 meV/f.u., and obtain the results in Fig 3b. (See Supp. Note 2 and Supp. Fig. S2 for details on the calculation of these loops.) The AFE→FE switching, and the FE→AFE back-switching, occur when the applied field lowers the corresponding energy barrier below the thermal activation energy. At 300 K, this occurs for 40 MV/cm (switching) and 20 MV/cm (backswitching), and we observe a sizeable hysteresis. The efficiency of the material as an AFE capacitor would be about 50 %.
The computed switching fields are very large compared to those that can be applied experimentally in similar oxides before inducing leakage (i.e., about 1 MV/cm). To understand the implications of this result, let us first note that first-principles estimates like ours are known to exaggerate FE coercive fields by factors of up to two orders of magnitude 30,31 , probably because they miss effects (e.g., easier nucleation of the field-induced phase at defects) that play an important role at controlling the transformation kinetics. However, in the case of a AFE↔FE transformation, the coercive bias must be as large as to equalize the energies of the polar and antipolar states. Our simulations do suggest that a very strong field (of about 28 MV/cm) is needed to achieve this in KVO 3 ; hence, notwithstanding possible inaccuracies in our estimate, it seems unlikely KVO 3 can be experimentally switched. Having said this, as we discuss below, we have reasons to believe that there are promising ways to optimize the switching characteristics of KVO 3 and related compounds, e.g. by means of appropriate chemical substitutions, considerably reducing the fields required to achieve AFE behaviour.
Discussion.-Let us start by commenting on three possibilities for AFE optimization and tuning that these materials offer. First, our simulations indicate that the behaviour of RbVO 3 and CsVO 3 is very similar to that of KVO 3 . The ZZ-A phase is also the ground state for these compounds, and the ZZ-F phase is a metastable polymorph with an out-of-plane spontaneous polarisation of 0.078 C/m 2 and 0.052 C/m 2 , respectively. The decrease in the polarisation, as compared to the result for KVO 3 , can be ascribed to the larger unit cell volume. The ZZ-F state of RbVO 3 shows a non-negligible in-plane component (0.037 C/m 2 ), while for CsVO 3 the polarisation is fully out-of-plane. The calculated energy barrier between the ZZ-A and ZZ-F states decreases with increasing size of the A cation (being of 128 meV/f.u. and 112 meV/f.u. for RbVO 3 and CsVO 3 , respectively). Further, the energy difference between the ZZ-F and the ZZ-A states follows the same trend, being of 117 meV/f.u.and 97 meV/f.u. for RbVO 3 and CsVO 3 , respectively (see Suppl. Note 3 and Suppl. Fig. S3).
Second, regarding the possibility of having different B cations, further calculations indicate that ZZ-A is a metastable structure for all the alkali tantalates and niobates, in particular for the well-studied FE KNbO 3 . Therefore, a morphotropic phase boundary (MPB) between the FE perovskite phase and the ZZ-A phase must exist for the solid mixture KV 1−x Nb x O 3 . Our preliminary studies indicate that such an MPB occurs between x=0.375 and x=0.5. At such compositions, the energies of the AFE (ZZ-A) and FE (KNbO 3 -like) states become very close (cross). Hence, these solid solutions naturally provide us with AFE and FE states that are very close in energy -thus solving the main difficulty mentioned above to obtain AFE↔FE switching at moderate fields -, and seem ideal candidates to yield AFE materials with optimized properties.
Third, as already mentioned, we find polymorphs with different tetrahedral arrangements, forming BM and RG patterns. Interestingly, for the three alkali vanadates considered, the most stable dipolar order is FE in the BM case. Moreover, the (shear) strains are distinct for the ZZ, BM and RG geometries (see the different structures of the A cation sublattice in Fig. 2). Hence, different chain arrangements, and in turn the polar order, may be accessible by growing these materials on appropriate substrates that impose suitable epitaxial conditions. Further, as can be appreciated from Figs. 2e and 2g, the ZZ-F state is more square in-plane than the ZZ-A, which suggests that a suitable substrate may allow us to tune the corresponding energy difference.
It is also worth noting some chemical and structural aspects of the compounds studied in this work. The origin of the peculiar polymorphs found seems to be the small size of the V 5+ cation (nominally, 0.54Å for V 5+ in a 6-fold coordination, and 0.355Å in a 4-fold coordination) 32 relative to the O 2− ionic radius (1.35Å) 32 . The size difference is such that all vanadates lie below the octahedral limit 33 , which is known to lead to lower B-cation coordination in perovskites [33][34][35] . Further, we find that the V-O bond lengths change by less than 1 % across the switching process between the ZZ-A and ZZ-F phases. In fact, the deviations in V-O bond lengths among all the ZZ, BM, and RG polymorphs of KVO 3 are below 1 %, and even the differences among the three studied vanadates remain below 1 %. The V-O bonds, and consequently the local dipoles, thus prove to be very rigid 36 , suggesting that in these materials phase transitions involving dipoles will probably be of the orderdisorder type 37 . More importantly, this bond stiffness also ensures that the dipoles will not vanish; therefore, these compounds can be viewed as model AFEs whose behaviour is analogous to that of antiferromagnets 38 .
Also of note is the switching in these compoundsby quasi-rigid rotations of VO 4 molecular-like groups -, which is rather unique as compared to similar transformations in inorganic FE and AFE materials. Indeed, the identified switching path suggests that, in these materials, such transformations will typically occur through many steps. The present calculation shows a switching in only two steps, which seems a direct consequence of the finite size of the simulation box employed; however, larger simulation cells would probably reveal a many-step process. This is strongly reminiscent of memristors, in which the electric resistance can be tuned quasi-continuously by taking advantage of a controllable multi-step transformation. Therefore, our results suggest that pyroxene-like AFEs could find application in memristor devices 39 .
Finally, the findings here reported hint at a promising strategy to discover further AFE materials. Pyroxenes, pyroxenoids 40 3 49 , with structures akin to that of KVO 3also contain chains of irregular oxygen tetrahedra. Since such tetrahedra host an electric dipole, these compounds can be viewed to present an antipolar ground state, and are candidates to display AFE behaviour. The situation is similar to that of BiVO 4 , an extensively investigated material (for its catalytic properties 50 ) that is formed by irregular VO 4 tetrahedra whose corresponding electric dipoles order in an antipolar pattern, and which has recently been proposed as a possible AFE 51 . Hence, we hope the present work will stimulate experimental and theoretical activities to explore this intriguing possibility, namely, that some of the best-known and most-abundant minerals on Earth may be antiferroelectrics in disguise! Acknowledgements. Work funded by the Luxembourg National Research Fund through the project IN-TER/ANR/16/11562984/EXPAND/Kreisel. We would like to thank Enric Canadell (ICMAB-CSIC) for his valuable comments and for pointing to Ref. 27 Author Contributions. Work carried out by H.A. and supervised by J.Í.
Methods. In this work we employ first-principles calculations within the density functional theory (DFT) framework as implemented in the Vienna Ab-initio Simulation Package (vasp) 52,53 to obtain the crystal structure and relative energies of the different polymorphs and compounds. The plane-wave cut-off for the basis set is set to 500 eV in all cases. We choose the Perdew-Burke-Ernzerhof functional modified for solids (PBEsol) 54 as the exchange-correlation. For perovskite primitive cell calculations (5-atom cell) a 6 × 6 × 6 Monkhorst-Pack 55 grid is employed for sampling the Brillouin zone, which yields well-converged results for the three alkali vanadates. The primitive cell of the ZZ-A ground state of KVO 3 , RbVO 3 and CsVO 3 is a √ 2 × 2 √ 2 × 1 supercell with respect to the ideal perovskite; for such simulation cells we employ a 4 × 3 × 6 k-point sampling. In all struc-tural relaxations the crystals are allowed to fully relax until atomic forces become smaller than 0.01 eV/Å. The molecular dynamics calculations are carried out within the canonical ensemble using a Langevin thermostat, as implemented in vasp. A 2 × 2 × 2 supercell is employed for this purpose (with respect to the 5-atom cell perovskite), and the k-sampling is reduced to 2 × 2 × 2 to speed up the calculations. A wide range of temperatures is studied, from 10 K to 1000 K.
The space groups of the studied crystal structures are determined employing the spglib library through its im-plementation in the phonopy package 56 .
We used the vesta visualization package 57 to prepare some of the Figures and Supplementary Videos.
The calculations of the switching energy barrier and the connection between the perovskite structure and the ZZ-A state are obtained through the NEB method 58 along with the climbing-NEB (cNEB) 59 modification as implemented by the Henkelman group in the Virtual Transition State Theory (vtst) tools package for vasp 59 . A total of 19 (7) images are employed in the NEB calculations of the switching energy barrier (connection to the perovskite structure).