Microscopic picture of paraelectric perovskites from structural prototypes

We show with first-principles molecular dynamics the persistence of intrinsic $\langle111\rangle$ Ti off-centerings for BaTiO$_3$ in its cubic paraelectric phase. Intriguingly, these are inconsistent with the Pm$\bar 3$m space group often used to atomistically model this phase using density functional theory or similar methods. Therefore we deploy a systematic symmetry analysis to construct representative structural models in the form of supercells that satisfy a desired point symmetry but are built from the combination of lower-symmetry primitive cells. We define as structural prototypes the smallest of these that are both energetically and dynamically stable. Remarkably, two 40-atom prototypes can be identified for paraelectric BaTiO$_3$; these are also common to many other $AB$O$_3$ perovskites. These prototypes can offer structural models of paraelectric phases that can be used for the computational engineering of functional materials. Last, we show that the emergence of B-cation off-centerings and the primitive-cell phonon instabilities are controlled by the equilibrium volume, in turn dictated by the filler A cation.

Compounds with the perovskite structure are a versatile class of functional materials exhibiting a wide range of properties, such as superconductivity [1], catalysis [2], photovoltaic energy harvesting [3] and ferroelectricity [4,5]. When ferroelectric, perovskites sustain a spontaneous polarization that can be switched with an electric field; as the temperature is raised, there is a transition above the Curie temperature to a paraelectric phase that has no net polarization. Early studies of BaTiO 3 (a prototypical ABO 3 ferroelectric perovskite) suggested for these transitions a microscopic "displacive" model, in which local displacements of the B-cation (titanium) align with the macroscopic polarization [6,7]. For BaTiO 3 this is along the 111 direction in the rhombohedral ground state; as the temperature increases there is a transition to an orthorhombic phase above 183K, with the polarization along 110 , then to a tetragonal phase above 278K, with the polarization along 100 , before reaching the paraelectric cubic phase above 393K, with no net polarization [7]. The results from diffuse X-ray scattering for all phases but the rhombohedral one [8,9] are somewhat inconsistent with such a displacive model. This has led to the application of the order-disorder model for the transitions [10,11] in which local polar displacements, driven by the pseudo Jahn-Teller effect [12], are in different ordered arrangements in the ferroelectric phases at low temperatures, and become disordered in the paraelectric phase. These two models can be reconciled if one considers the time-averaging inherent to most characterization techniques, which can effectively wash out the local displacements and present a higher-symmetry structure where the averaged displacements are aligned with the macroscopic polarization or cancel out [13].
To elucidate the microscopic picture of paraelectricity in these perovskites, we performed Car-Parrinello molecular dynamics (CPMD) simulations of cubic BaTiO 3 , finding clear microscopic evidence of Ti off-centerings, and associated dipoles, along the 111 directions which persist well above the Curie temperature, consistent with the order-disorder model. This is clearly apparent in Fig. 1a, where we present the results of CPMD simulations at 315K for a 4×4×4 cubic supercell; the histogram for the Ti displacements projected onto one of the equivalent [100] planes shows how the Ti atoms always occupy off-center 111 positions, rather than sitting at the center of their surrounding oxygen octahedron [35](for the Methodology see SI Sec. 1 [36]; for the associated data, see Ref. 37 on the Materials Cloud [38]). We observe these off-center displacements up to temperatures around 450K; furthermore, they can be suppressed with compressive hydrostatic strain, resulting in an isotropic distribu- tion (Fig. 1b) in agreement with experimental measurements [39,40]. This observation will be relevant to the later discussion of volume effects and the role of A-site cations.
Inspired by these results, we aim here to systematically explore the microscopic structure of the paraelectric phase of BaTiO 3 , to extend this exploration to other perovskites, and to lay the groundwork for a systematic analysis of phases that can possess "hidden order", including other ferroelectric or magnetic systems [41][42][43] or those displaying higher-order couplings [44]. With this goal in mind, we introduce first the concept of microscopic templates, defined as lower-symmetry supercells that preserve a desired point symmetry (e.g., cubic). We then define microscopic prototypes as the smallest of these templates that are both energetically and dynamically stable (i.e., lower in energy, per formula unit, than the highersymmetry primitive cell and with real, positive phonon dispersions), thus minimizing computational cost while identifying the highest-symmetry, stable structures possessing the requisite symmetry. This approach is distinct yet complementary to that of special quasirandom structures (SQSs) when used to describe a polymorphous network, in which a single, large SQS exhibiting many local motifs is used. In this context, slightly different from the original development of SQSs to characterize disordered alloys with first-principles calculation [45,46], SQSs have been recently used to study paramagnetic phases [47,48] and complex perovskite-based systems [34,[49][50][51][52]. In this work we develop instead a symmetry-based analysis and workflow, enumerating all possible supercells (up to a given size) with a desired point symmetry. In this way we identify not just local motifs, but more complex orderings which respect the desired global point symmetry. We describe it in the following and apply to structural microscopic prototypes, but these concepts can be equally applied to magnetic or electronic prototypes.
To identify structural prototypes we use groupsubgroup relations, as discussed in Ref. 53, to systematically enumerate all microscopic templates; here, we take the case of the cubic ABO 3 perovskite with space group Pm3m (international number 221), where the 1a, 1b and 3c Wyckoff positions are occupied by the A, B, and O atoms, respectively. For each cubic subgroup of Pm3m, we define a cubic microscopic template as a supercell that can host symmetry-allowed displacements of A, B, and O atoms relative to their positions in the highsymmetry parent structure (group Pm3m) with no net polarization (see SI Sec. 1 [36]for further details). Using 2×2×2 supercells of 40 atoms, we find 27 distinct cubic microscopic templates of group Pm3m, 10 of which host only oxygen displacements, while the remaining 17 allow the A and/or B cations to displace as well. Table I summarizes the subgroups in which the B cations can displace; see SI Sec. 2 & 3 [36]for the complete list as well as a list duplicates that correspond to microscopic templates with higher symmetry. The same analysis can be applied to supercells of any desired size, but we find that in BaTiO 3 these 2×2×2 supercells are already sufficient to identify structural prototypes.
We then determine which of these microscopic templates, if any, are energetically stable. Using variablecell first-principles relaxations performed with Quantum ESPRESSO [55,56]  instabilities are marked in red with a cross, a filled circle and an empty square, respectively. (b) Forty-atom undistorted supercell; X , M and R labels indicate the X, M and R points of the 2×2×2 supercell, respectively (the same labels are used also in panels c and d). All three instabilities marked in panel a fold at Γ in this supercell. (c&d) 40-atom supercell with the 4+4 and 2+6 displacement pattern, respectively. In panels c and d (top), red arrows indicate the direction of atomic displacements (only shown for B cations for clarity). We refer the reader to the SI Sec. 5 [36] for the displacement pattern associated with each unstable mode and the other modes that contribute to the 4+4 and 2+6 displacement patterns, which can also be visualized with the Interactive phonon visualizer tool on the Materials cloud [38]. symmetry remains cubic (See SI Sec. 1 [36] for further details). Remarkably, we find that two of the microscopic templates relax to supercells with non-trivial displacement patterns of the B cations; moreover, they display stable phonon dispersions across the entire Brillouin zone (see Fig. 2c,d). The remaining 25 templates either relax back to (the 2×2×2 supercell of) the fiveatom primitive cell, well-known to be dynamically unstable (see Fig. 2a,b) [61,62] or to one of these two non-trivial displacement patterns. hese energetically and dynamically stable structures are the structural prototypes. As they are locally stable structures of the 0K potential energy landscape they serve as minimal models possessing the signature of the paraelectric phase, namely a global cubic symmetry but with local Ti displacements. These displacements, driven by local chemistry, can then also acquire correlations (e.g. linear chains [8]) that can be studied with large-scale molecular dynamics simulations [27,30,63].
The two structural prototypes have symmetry I43m and Pa3, respectively (see Table I); their structure and B-atom (Ti) displacement patterns are shown in Fig. 2cd. We name these two prototypes 4+4 and 2+6 (for I43m and Pa3, respectively), since considering any Ba atom, in the 4+4 (2+6) structure there are 4 (2) surrounding Ti atoms that displace toward it, while the remaining 4 (6) displace outwards. We note that the 4+4 structure (I43m) has been previously discussed in the work of Zhang et al. [32]. The 4+4 and the 2+6 structural prototypes are lower in energy than the undistorted cubic structure by 11 and 15 meV/formula unit, respectively. Furthermore, there is an energy barrier of only 3 meV/formula unit between these two structural prototypes (as found by nudged-elastic-band calculations, see SI Sec. 4 [36]), suggesting that thermal fluctuations of the off-centerings do not require to go through the highsymmetry structure.
We contrast the phonon dispersions of the highsymmetry structure (Fig. 2a,b) with that of the 2 prototypes (Fig. 2c,d). The five-atom primitive cell displays instabilities at the zone-center Γ, belonging to the irreducible representation (irrep) Γ − 4 , and at the zoneboundary points X and M (irreps X + 5 and M − 2 , respectively). To gain further insight into the 4+4 and 2+6 patterns we analyze these with respect to the irreps of the five-atom-cell phonons using the ISODISTORT software [64,65]. We find that the displacements of both prototypes contain a mode with the symmetry of an unstable zone-boundary mode. Specifically, the 4+4 prototype can be constructed by adding the displacements having the symmetry of the M − 2 and M + 1 irreps, while the 2+6 prototype originates from the X + 5 and M + 5 irreps (see SI Fig. S2 [36] for the M + 1 and M + 5 modes). Most importantly, out of the 27 distinct cubic templates, the 4+4 (I43m) and 2+6 (Pa3) are the only ones with a displacement pattern that is constructed, in part, from a mode with the symmetry of an unstable mode of the parent structure, resulting in an appealing one-to-one correspondence between unstable zone-boundary phonon modes and prototypes with stable displacement patterns in 2×2×2 cubic supercells. Wenote that the displacement patterns must occur in combination with another mode in a cubic structure as they do not possess a cubic point symmetry. Furthermore, the Γ − 4 mode is the polar instability and can only occur in lower-symmetry polar phases of BaTiO 3 , which are therefore non-cubic.
We investigate in more detail the zone-boundary modes and the stability of the structural prototypes as a function of volume, prompted by the disappearance of the Ti off-centering under pressure in our CPMD simulations ( Fig. 1b) and in experiments [39,40]. We find that with increasing pressure the magnitude of the Ti displacements decreases for both prototypes, and disappears when the lattice parameter is reduced by ≈ 1.6%, as reported in Fig. 3a. We find that the Ti displacements as a function of volume can be fit by a double-well potential where the quadratic coefficient depends linearly on volume and changes sign at the onset of the displacements (see SI Sec. 6 [36]). This suggests that at least one phonon mode associated with this structural prototype becomes unstable at the same volume where the Ti displacement becomes energetically favorable. In Fig. 3b we plot as a function of the lattice parameter the phonon frequencies for the q = Γ, X and M modes that are unstable in the five-atom primitive cell (irreps Γ − 4 , X + 5 and M − 2 , respectively). We find that expanding the volume further softens these modes, while applying pressure stabilizes them, in agreement with previous calculations [66]. The fact that the Γ − 4 mode also stabilizes at a lower lattice parameter is indicative of a pressure at which the system could be ferroelectric below a critical temperature, but no Ti displacements would be observed in the cubic paraelectric phase. Notably, the M − 2 and X + 5 modes become unstable at the same lattice parameter where the 4+4 and 2+6 displacement patterns respectively emerge (gray arrows in Fig. 3b).
Thus, the 4+4 and 2+6 prototypes originate from the unstable M − 2 and X + 5 modes, which do not involve Acation displacements. To test the effect of the A cation we extend the study to PbTiO 3 , SrTiO 3 , and CaTiO 3 . We report in Fig. 4 the results for the 4+4 prototype, highlighting a universal trend where the B-site displacement as a function of lattice parameter is broadly independent of the chosen A cation. The stability of the prototype, and thus the nature of the paraelectric phase, is instead determined by the equilibrium lattice parameter -indicated in Fig. 4 by arrows -which is largely determined by the A cation. For Pb, Sr, and CaTiO 3 the lattice parameter is smaller than the critical value at which the displacement pattern becomes energetically favorable, ∼3.95Å. For all titanates studied, the displacement pattern onset occurs at the lattice parameter at which the M − 2 mode becomes unstable. A similar picture emerges for the 2+6 pattern (except for CaTiO 3 , due to its significantly smaller lattice parameter) -see SI Sec. 7 [36].
Testing a broader range of 49 perovskites from Ref. 67 shows that B-site off-centerings along 111 directions provides prototypes at the relaxed equilibrium lattice parameter not only for BaTiO 3 , but for most zirconates, niobates, and tantalates, CaHfO 3 and BiScO 3 , as reported in SI Sec. 8 [36]. However, the energetic stability of these prototypes as a function of lattice parameter is B-site specific. This is expanded on in SI Sec. 9 [36] where we further investigate the stability of the 4+4 and 2+6 displacement patterns as a function of lattice parameter in the titanates, niobates and zirconates, demonstrating the universality of the occurrence of B-cation displacements, and their strong, family-specific volume dependence.
The relationship we have observed between the unstable zone-boundary phonons of the primitive cubic structure and the displacement patterns as a function of lattice parameter indicates that, at a given volume, one could use the unstable phonon modes to predict which microscopic templates would result in structural prototypes. To verify the robustness of our conclusions against the choice of DFT functional, we tested for BaTiO 3 the dependence on the functional, finding that the Ti displacement amplitudes are independent of the functional choice (see SI Sec. 10 [36]), but that the functional determines the equilibrium volume. This highlights the need to choose a functional that accurately reproduces the experimental lattice parameter in order to correctly predict which prototypes occur at the equilibrium volume.
In summary, motivated by the observed persistence of 111 Ti off-centerings above the Curie temperature in our CPMD simulations of BaTiO 3 , we systematically identify microscopic structural prototypes of the paraelectric phase, i.e., the smallest supercells with cubic point symmetry that are simultaneously energetically and dynamically stable. These cubic prototypes, hosting stable local dipoles due to the 111 Ti displacements in a cubic paraelectric phase, are found through a symmetry analysis exploring all possible 40-atom microscopic templates, followed by density-functional theory and densityfunctional perturbation theory calculations to assess energetic and dynamical stability. Moreover, we highlight how off-centering amplitudes are strongly dependent on volume, and relate their patterns to the zone-boundary unstable phonons of the five-atom undistorted primitive cubic cell, suggesting a predictor for the identification of such prototypes. These cubic prototypes would be challenging to identify without the present symmetry-based approach, due to the combinatorial complexity of large supercells and the attractive basin of the rhombohedral five-atom ground state associated with the polar instability. We highlight that these prototypes can serve as minimal models of the paraelectric phase in first-principles calculations of response functions with the correct tensorial symmetry as they provide a faithful microscopic representation with key features: the persistence of local Ti displacements and the appropriate macroscopic point group.
We finally emphasize that this approach is general. Beyond its extension to study the prevalence of the B-site 111 off-centerings in ABO 3 perovskites, it can be used in any crystalline system to find candidate templates and efficiently search for prototypes that are local minima in the potential-energy surface, providing an indepth study of, for example, the electronic or magnetic properties of a polymorphic system. This approach lays the foundation to investigate dynamics, thermodynamics and chemical substitutions, as these prototypes could be used to capture subtle details of the energy landscape and to provide models to study the properties or transitions of disordered phases, such as alloys, paramagnetic phases, or defects in paraelectric phases. Our first-principles DFT-based calculations are performed using Quantum ESPRESSO [1,2] with the PBEsol [3] functional, RRKJ pseudopotentials [4] and wavefunction and charge-density energy cutoffs of 60 and 600 Ry, respectively.

I. Car-Parrinello Molecular Dynamics
CPMD calculations [5] are carried out using the cp.x code of Quantum ESPRESSO in the NVE ensemble (constant number of particles, volume and energy) on a 4 × 4 × 4 supercell (320 atoms) using a Γ-only sampling, where the initial atomic positions are chosen so that, after an initial equilibration (1.5 ps), the average temperature is the one reported in the figure caption. The lattice parameter is also increased with respect to the T = 0 K value so as to account for the experimentally-measured thermal expansion [6]. The wavefunction and charge-density energy cutoffs is set to 40 and 320 Ry, respectively. A time step of 0.25 fs and an electron mass of 400 m e are used.

II. Symmetry analysis to determine the microscopic templates
We search for all the subgroups with a desired point symmetry (here, cubic) of the high-symmetry space group with unique Wyckoff splittings in a recursive fashion. We use and compare both the ISOTROPY command-line tool [7] and the tools available on the Bilbao Crystallographic server [8,9], in particular the CELLSUB program to obtain the list of subgroups, the conjugacy classes and the transformation matrices, then used as input to the WYCKSPLIT program [10] to get the splitting of the relevant Wyckoff positions. These two set of rules yield the same results and allow us to generate the group-subgroup relationship used to build Tables 1 (main text), S1 and S2. To confirm the space group, removing the duplicate microscopic templates, we perform a final check with spglib [11].

III. Density functional theory calculations to search for structural prototypes
Using a random displacement on the order of 0.05Å for each Wyckoff-position free parameter, we construct cells for each of the 27 cubic microscopic templates associated with each unique subgroup for the cubic phase of BaTiO 3 and perform variable-cell relaxations constraining cubic point symmetry. At every ionic minimization, the forces were symmetrized, effectively enforcing symmetries detected at the onset of the calculation. This ensures that the system reaches the minimum-energy structure with the constraint of possessing the symmetry of the parent space group (e.g., a given cubic spacegroup), avoiding that the system relaxes to one of the lower-symmetry structures (e.g., the rhombohedral structure that is the ground state at T = 0 K for BaTiO 3 ). This enforces the cubic point symmetry possessed by all of the 27 microscopic templates. The calculations are managed using the AiiDA framework [12,13], a high-throughput platform that allows to automatically launch, retrieve, parse and organize the calculations, storing the results in a database and automatically managing sequences of calculations via its workflow engine. Given a microscopic template as a starting structure, we relax both the lattice parameters and the internal coordinates with a force tolerance of 5×10 −5 eV/Å and an energy tolerance of of 3×10 −11 eV using a 3 × 3 × 3 k-mesh in the 2 × 2 × 2 supercells (and equivalent meshes in different cell sizes, e.g., 6 × 6 × 6 in the 5-atom unit cell). To calculate the potential energy landscape, we perform a nudged-elastic-band (NEB) calculation [14] using the metastable 4+4 and 2+6 structures and the undistorted structure as the constrained points. To calculate the phonon dispersion we use density functional perturbation theory (DFPT) [15,16] as implemented in Quantum ESPRESSO with a 2 × 2 × 2 q-mesh imposing the acoustic sum rule. For the phonon dispersions, 50 points are used along each segment for the Fourier interpolation. We use ISODISTORT [17,18] from the ISOTROPY package to analyze the symmetry of the unstable modes and of the displacement patterns of the metastable lower-symmetry structures, i.e. the structural prototypes.   (221) possessing non-trivial splittings of the 1a, 1b and 3c Wyckoff positions are reported which each correspond to a microscopic template, including the 10 subgroups reported in the main text. We list the subgroup (international short symbol and number) followed by the subgroup index (in square brackets); the transformation matrix Tx; and the splittings of the three relevant Wyckoff positions. Tx is defined as the transformation from the primitive cell to the 2 × 2 × 2 cell (with matrix equal to twice the identity), and with an origin shift (x, x, x). In the first ten subgroups (above the dividing line), only oxygen atoms can displace from the high-symmetry structure. In the trivial subgroups, the resulting split Wyckoff positions have no degrees of freedom. The cell of the representation of the first four subgroups is the primitive cell and of the remaining six is the 2 × 2 × 2 cell. The duplicate subgroups, identified using spglib [11], when only the Wyckoff positions originating from the 1a, 1b and 3c positions in the parent group are occupied, actually have higher symmetry and fall back in one of the subgroups of Table S1, reported in the right-most columns. For these duplicate subgroups, there are 12 in which only the oxygens have degrees of freedom, and 3 in which the A and/or B sites have degrees of freedom. We report the subgroup using the international short symbol; the subgroup number followed by the subgroup index (in square brackets); the transformation matrix Tx; and the splittings of the three relevant Wyckoff positions. Tx is defined as the transformation from the primitive cell to the 2 × 2 × 2 cell (with matrix equal to twice the identity), and with an origin shift (x, x, x). For the four subgroups with primitive cell representations, the transformation matrix is trivial. The reported subgroup index and transformation matrix, in all cases, is given with respect to the parent group Pm3m (221) with the 1a, 1b, and 3c positions occupied. FIG. S1. Potential energy landscape calculated via two nudged-elastic-band calculations between the undistorted cubic BaTiO3 (used as the reference energy) and the relaxed 4+4 structure, and between the relaxed 4+4 and 2+6 structures. The calculations show that the barrier between the 4+4 and 2+6 structures is lower than the energy difference with respect to the undistorted cubic BaTiO3. We note that the non-centrosymmetric lower-temperature phases of BaTiO3, associated with the polar instability Γ − 4 are lower in energy that the 4+4 and 2+6 structural prototypes. Each mode is shown in its own primitive cell. We note that the Γ − 4 , X + 5 , and M + 5 modes are all two-fold degenerate. All data (inputs and outputs) from our self-consistent DFT calculations as well as the DFPT phonon calculations are available in a Materials Cloud Archive entry at the following DOI: https://doi.org/10.24435/materialscloud:pg-50, within the file materials 5atom phonon.zip file, in the directory materials 5atom phonon/BaTiO3/3.98; the displacement patterns associated with these phonon modes can be visualized in an interactive manner, using the "Interactive phonon visualizer" tool on the Materials Cloud, available at the following address: https://www.materialscloud.org/work/tools/interactivephonon. Using a simple double-well model we can fit the Ti displacement as a function of volume. Assuming a simple double-well potential where the central point is the undistorted cubic structure and the minima of the wells correspond to the 4+4 (or 2+6) metastable structure, we define where E is the total energy (with the zero set at the energy of the undistorted structure), x is the set of coordinates that takes the undistorted structure to the metastable state (we will use here the magnitude of Ti displacements), and the coefficients A and B depend on the volume of the system. Solving for the stationary points of Eq. (61), we can find the minima at (x(V ),Ẽ(V )), and then invert the solutions to obtain: From our DFT results of the relaxed displacement patterns as a function of lattice parameter (or, equivalently, volume) we can extract values forx(V ), and the corresponding total energyẼ(V ), ans thus obtain data points for A(V ) and B(V ) as a function of the volume from Eq. (62). These are reported in Fig. S4a,b (Fig. S5a,b) for the 4+4 (2+6) displacement pattern. Because of the dependence of the data points, we fit them with the following functional forms: A(V ) = A 0 · (V − V 0 ) (with A 0 and V 0 two constants) and B(V ) with a constant B 0 ; the resulting fits are reported in the plots (as well as a 0 = 3 √ V 0 , for convenience). In Fig. S4c (Fig. S5c) we then report the data forx as a function ofṼ for the 4+4 (2+6) pattern, where the solid line is the analytical expression for , obtained using the fitted values for A 0 , V 0 and B 0 in panels a and b, that reproduces very well the DFT data points. Similarly, in Fig. S4d (Fig. S5d) we plot the energy differenceẼ of the 4+4 (2+6) pattern with respect to the undistorted structure, and the corresponding analytical curvẽ    We explore the stability of the 4+4 and 2+6 displacement patterns for other perovskite oxides of the titanate family, namely PbTiO 3 , SrTiO 3 , and CaTiO 3 . For the sake of completeness, we include here also Fig. S6 with the results for BaTiO 3 , i.e., the same as Fig. 3 of the main text. These figures show that the 4+4 (2+6) displacement pattern onset coincides with the lattice parameter at which the M − 2 (X + 5 ) irrep of the 5-atom undistorted cubic cell becomes unstable, indicated by the gray arrows. For clarity, these figures do not include other zone-boundary phonon modes that are unstable within this range of lattice parameters; however, they are discussed in the material-specific sections below. The Γ − 4 irrep does not have any cubic subgroups and, therefore, does not correspond to a metastable cubic structure, so we omit it from the other plots.

I. Lead titanate
At the equilibrium lattice constant of cubic PbTiO 3 computed with PBEsol (3.93Å), we find unstable zone-boundary phonons at q = M and R as well as an unstable zone-center mode. These modes, omitted from Fig. S7 as they are not involved in the 4+4 or 2+6 displacement patterns, transform with the symmetry of irreps M + 2 , R − 5 , and Γ − 4 , respectively. The irreps at R and Γ have no cubic subgroups, and the mode that transforms as M + 2 is dominated by oxygen. Using instead a lattice constant of 3.97Å (the experimental cubic lattice parameter for the cubic phase, stable at temperatures higher than 763K [19]), we find additional unstable phonons at q = X and M that transform like X + 5 and M − 2 , respectively (like for BaTiO 3 ). These are the modes that are present in the 2+6 and 4+4 displacement patterns, respectively. As for BaTiO 3 , the Ti displacement onset for the 2+6 (4+4) displacement patterns correlates with the X + 5 (M − 2 ) mode becoming unstable, see arrows in Fig. S7. The stability of the 4+4 and 2+6 displacement patterns at the experimental lattice parameter of cubic PbTiO 3 suggests that the Ti atoms will displace along the 111 directions and is consistent with the non-negligible degree of order-disorder that has been observed experimentally [20][21][22][23][24]. As we continue to increase the lattice constant we find that additional modes become unstable, in particular the ones transforming like M − 5 (at ∼3.975Å) and like X − 5 (at ∼4.00Å).

II. Strontium titanate
At the equilibrium lattice constant of cubic SrTiO 3 computed with PBEsol (3.89Å), we find unstable zone-boundary phonons at q = M and R as well as an unstable zone-center mode. These modes, again omitted from Fig. S8 as they are not involved in the 4+4 or 2+6 displacement patterns, transform with the symmetry of irreps M + 2 , R − 5 , and Γ − 4 respectively. Again we show that the Ti displacement onset for the 2+6 (4+4) displacement patterns correlates with the X + 5 (M − 2 ) mode becoming unstable, see arrows in Fig. S8.

III. Calcium titanate
At the equilibrium lattice constant of cubic CaTiO 3 computed with PBEsol (3.85Å), we find unstable zone-boundary phonons at q = M and R as well as an unstable zone-center mode. These modes, again omitted from Fig. S9 as they are not involved in the 4+4 or 2+6 displacement patterns, transform with the symmetry of irreps M + 2 , M − 5 R − 5 , and Γ − 4 respectively. The displacement pattern associated with R − 5 has no cubic subgroups, while the unstable M modes do correspond with cubic subgroups. The displacement pattern associated with M + 2 only involves motion of the oxygen atoms; the displacement pattern associated with M − 5 involves all of the atoms; however, the displacements of the calcium and oxygen are more pronounced than the displacement of the titanium. This implies that, in the paraelectric phase, displacements from the high-symmetry structure would be dominated by the oxygen and calcium with only a small contribution from the titanium. This is not surprising as the tolerance factor of CaTiO 3 is well below 1 and has a well-known P bnm (spacegroup 62) orthorhombic phase that can be constructed through a linear combination of the M + 2 and R − 5 modes. It is worth noting that, according to the M − 5 mode, the titanium atoms are still restricted by symmetry to displace along the local 111 directions. Again we show that the Ti displacement onset for the 2+6 (4+4) displacement patterns correlates with the X + 5 (M − 2 ) mode becoming unstable, see arrows in Fig. S9. However, we see that the onset for the 2+6 pattern occurs at a smaller lattice parameter than the other titanates. Moreover, we see that as the lattice parameter increases, the titanium displacement reaches a maximum and then decreases, unlike all the other titanium displacement curves of the other titanates. In contrast to the behavior of the Ti displacement as a function of the lattice parameter, the X + 5 mode of the 5-atom primitive cell continues to soften monotonically, while at the same time, the other mode present in the structure (M + 5 , which only displaces oxygen atoms), remains stable. Using ISODISTORT to analyze the symmetry of the distortion present in the 2+6 displacement patterns under larger tensile strain, we find that the displacement associated with the X + 5 irrep (a six-dimensional irrep) becomes dominated by oxygen. In Fig. S10 we plot the supercell mode amplitude as a function of lattice parameter for the X + 5 irrep present in the 2+6 (40-atom) structures. In green are the supercell mode amplitudes associated with the oxygen displacements that transform according to the A 2u (solid circles) and E u (open squares) irreps of the point-group and in purple is the supercell mode amplitude associated with the titanium displacement. The total supercell X + 5 mode amplitude, in blue, is the norm of the three components, i.e., the square root of the sum of each component squared. To convert from supercell mode amplitude A s,Ti plotted in Fig. S10 (as defined within the ISODISTORT [17,18] program) to the titanium displacement ∆x Ti of Fig. S9 (the displacement inÅ from the cubic high-symmetry position along one coordinate), the following prescription is used: where a s is the supercell lattice parameter (i.e., the lattice parameter of the 2×2×2 supercell), n is a normalization factor, and m is the number of components of the displacement pattern associated with the irrep. Here, the number components for the titanium displacements is 24: the 8 Ti atoms displace along all three components, so we have m = 8 · 3 = 24. To investigate the universality of the paraelectric phases in other cubic perovskites, we explore the stability of the 4+4 and 2+6 orderings for the entire set of perovskites in the work of Armiento et al. [25]. The Goldschmidt tolerance factor t [19,26], defined as where r i are the atomic radii of the A-site cation, the B-site cation or the oxygen anion. We calculate the tolerance factor using the ionic radii from Ref. 27 taken from the database hosted on Ref. 28. In cases where the ionic radii for the relevant oxidation state or coordination number was not available, we used the closest one. In materials where the 4+4 and 2+6 prototypes are significantly lower in energy, the tolerance factor is t > ∼ 1.1, and in materials where only the only the 2+6 prototype significantly lowers the energy, the tolerance factor is t < ∼ 0.9 (Fig. S11). Specifically, Cs-based compounds have t > 1.1; Rb-based compounds have t ≈ 1.1; Ca-based compounds have t ≈ 0.84 − 0.9; Li-based compounds have t ≈ 0.8; Cu-based compounds have t < 0.8; and BiScO 3 has t ≈ 0.7. Note that a number of compounds stabilize the 2+6 and 4+4 prototypes, but the energy gain is quite modest, such as BaTiO 3 or KNbO 3 . In these cases we do not find a general trend of the tolerance factors. Move over, many of the silicates, that have a tolerance factor less than 0.9 or greater than 1.1, do not stabilize neither the 4+4 nor the 2+6 orderings further supporting the fact that we cannot predict what prototypes will be stable or not from the tolerance factor alone. We find that these large energy gains correlate with the change in lattice parameter between the cubic cell without displacements and the one with displacements; however, the energy gain is not due to the change in lattice parameter alone, and a significant energy gain occurs even when the total volume is kept fixed. A better understanding of the behavior of these systems can be obtained by investigating the occurrence of stable displacements as a function of the system volume, see Fig. S12 FIG. S11. Occurrence of cubic-symmetry displacement patterns in 2 × 2 × 2 supercells for the family of 49 perovskites from Armiento et al. [25]. Solid line: 4+4 displacement pattern; dashed line: 2+6 displacement pattern. The lattice parameter is fixed to that of the relaxed 5-atom cubic structure. The tolerance factor is given along the top axis (a) Off-site displacement of the B-site cation. (b) Energy difference per formula unit ∆E and relative lattice parameter change |∆a|, both with respect to the cubic cell without displacements.
Section 9: Stability of 4+4 and 2+6 orderings in the niobates, titanates and zirconates as a function of volume We investigate 12 different perovskites as a function of volume, focusing only on the 4+4 and 2+6 displacement patterns, to demonstrate the generality of the occurrence of local displacements in different perovskite families. The results are shown in Fig. S12 The phonon dispersion of cubic 5-atom KNbO 3 , which also stabilizes the 4+4 and 2+6 displacement patterns, has the same zone-center and zone-boundary instabilities [29] as BaTiO 3 at the equilibrium volume; however, the Nb off-centering as a function of volume does not match that of the titanates. In particular, for all niobates that we investigate here, also a volume compression can trigger the occurrence of local displacements, and there is a very narrow lattice parameter range (or even no range at all) in which displacements do not occur. This occurs as additional phonons become unstable under compression. We also present four zirconates, which have a different B-site displacement vs. volume relationship, where instead a compressive strain is needed to stabilize these two patterns, and at relaxed volume no displacements are stabilized (except for CaZrO 3 in the 2+6 pattern, where also an increasing volume stabilizes them, including at the relaxed volume). In all cases, we confirm that the volume dependence and stability of the structural prototypes in a given material cannot be inferred from equilibrium properties alone (like in Fig. S11), but a detailed understanding requires an analysis as a function of the system volume. The choice of the DFT functional essentially does not affect the magnitude of Ti displacements, as well as the onset as a function of the lattice parameter, for both the 4+4 and 2+6 displacement patterns; however, it does change the predicted equilibrium lattice constant (indicated by vertical arrows). In the figure we compare the PBEsol [3], PBE [30,31], and LDA [32] functionals. The PBEsol equilibrium lattice parameter (at T = 0 K) is the closest to the experimental lattice constant (at finite temperature) of ∼ 4.01Å [6,19] of cubic BaTiO 3 (stable above ∼390 K). FIG. S13. Onset of the cubic-phase instabilities as a function of the lattice parameter for the 4+4 and 2+6 structures of BaTiO3 using the PBEsol, PBE and LDA functionals. The relaxed lattice parameter for each functional in its lowest-energy cubic configuration is indicated by an arrow.