Anion order in oxysulfide perovskites: origins and implications

Heteroanionic oxysulfide perovskite compounds represent an emerging class of new materials allowing for a wide range of tunability in the electronic structure that could lead to a diverse spectrum of novel and improved functionalities. Unlike cation ordered double perovskites—where the origins and design rules of various experimentally observed cation orderings are well known and understood—anion ordering in heteroanionic perovskites remains a largely uncharted territory. In this contribution, we present and discuss insights that have emerged from our first-principles-based electronic structure analysis of a prototypical anion-ordered SrHf(O0.5S0.5)3 oxysulfide chemistry, studied in all possible anion configurations allowed within a finite size supercell. We demonstrate that the preferred anion ordering is always an all-cis arrangement of anions around an HfO3S3 octahedron. As a general finding beyond the specific chemistry, the origins of this ordering tendency are traced back to a combined stabilization effect stemming from electronic, elastic, and electrostatic contributions. These qualitative notions are also quantified using state-of-the-art machine learning models. We further study the relative stability of the identified ordering as a function of A (Ca, Sr, Ba) and B (Ti, Zr, Hf) site chemistries and probe chemistry-dependent trends in the electronic structure and functionality of the material. Most remarkably, we find that the identified ground-state anion ordering breaks the inversion symmetry to create a family of oxysulfide ferroelectrics with a macroscopic polarization >30 μC/cm2, exhibiting a significant promise for electronic materials applications.


INTRODUCTION
Perovskite compounds-with a generic formula ABX 3 -represent a versatile class of materials exhibiting tremendous synthetic flexibility. For instance, >90% of the metal ions in the Periodic Table are known to be stable in a perovskite-type structure 1,2 . Moreover, the anion sublattice X can be populated by several members of halide (X=F, Cl, Br, and I), chalcogenide (X=O, S, Se) or pnictide (X=N) families. A further possibility of partial substitutions at A, B, and X sites accompanied with a rich set of possible cation-and/or anion-orderings make this class a truly staggering collection of compounds. This diverse perovskite composition and ordering naturally leads to a rich portfolio of chemical, electronic and magnetic properties for applications ranging from catalysis to electronic devices and from sensors to energy storage [3][4][5][6] .
Little is yet known about the structure and properties of perovskite oxysulfides, despite the potential advantages of mixing these two anions. Although metal oxides are a mainstay of many devices owing to their stability and good dielectric properties, the increased covalency of metal-sulfide bonds offers the chance to narrow the bandgap and improve electronic transport properties. Indeed, sulfide perovskites are under active investigation as environment-friendly solar energy harvesters. Recent studies focused on chalcogenides have shown that ABX 3 compounds within this class-with A=Ca, Sr, or Ba; B=Ti, Zr, or Hf; X=S or Seare stable in a distorted perovskite structure with the bandgaps ranging from 0.3 eV to 2.3 eV [23][24][25][26] . Many of these semiconducting compounds possess direct bandgaps with high optical absorption coefficients-comparable to those of traditional optoelectronic semiconductors such as GaAs-with large valence and conduction band dispersions, suggesting good carrier mobility 23 . For some ABX 3 sulfides, competing non-perovskite phases have also been reported 24,[27][28][29][30] . More recently, a wide range of ABS 3 perovskites and non-perovskites were screened for relevant properties, such as thermodynamic stability, light absorption, charge mobility, and defect tolerance. The study identified LaYS 3 as a promising candidate 31 . As a first step in the exploration of perovskite oxysulfides, Perera et al 18 . recently synthesized a series of BaZr(O 1−x S x ) 3 oxysulfides over the wide range of x ∈ [0, 0.54] to show that all the compounds within this composition range are stable in a perovskite crystal structure with a systematic bandgap variation of >1 eV.
In general, the properties of perovskites are sensitive to small structural changes, such as rotations, distortions, and tilting of the metal-anion octahedra. In the case of oxysulfide perovskites (or more broadly in heteroanionic perovskites), the B-site octahedral coordination is controlled by two or more anions. Different anion configurations around a given metal ion could, in principle, lead to different levels of crystal field strength, bonding characteristics, and local polarization, which in turn dictate the ground-state structure and functionalities. In spite of this potential for rich structure-property relationships, investigations of anion order in mixed chalcogenide perovskites are missing.
Here, we report a first-principles investigation of anion order in AB(O 0.5 S 0.5 ) 3 perovskites, with A=Ca, Sr, Ba and B=Ti, Zr, Hf. We first analyze all 246 symmetry-unique configurations accessible within a 20-atom supercell of CaHf(O 0.5 S 0.5 ) 3 , and distill down various factors that control their relative energetics and bandgap variations. Specifically, we find that a cis-arrangement of O and S around the BX 6 octahedron is highly preferred over a transarrangement. We show that a combination of electronic, local strain, and electrostatic effects largely dictates the anion ordering in this class of compounds. The bandgap is controlled by the crystal field splitting imposed by the surrounding anion in the BX 6 octahedron, which primarily changes the position of the Hf dstates. We then generalize this insight by exploring stability trends for the larger set of AB(O 0.5 S 0.5 ) 3 oxysulfides with A ∈ {Ca, Sr, Ba} and B ∈ {Ti, Zr, Hf} to demonstrate that the trends can be rationalized using geometrical notions and ionic radii of the constituent ions. In particular, our computations predict that zirconates of Ca and Sr are thermodynamically stable and, therefore, should be amenable to experimental synthesis. At last, we demonstrate that the asymmetric and stable cis-arrangement of O and S anions in the perovskite structures can support ferroelectricity.

RESULTS
Bulk CaHfO 3 , CaHfS 3 , and CaHf(O 0.5 S 0.5 ) 3 chemistries Both CaHfO 3 and CaHfS 3 are known to be stable in a GdFeO 3 -type crystal structure with a Pnma space group and can be represented within a 20-atom orthorhombic supercell. Using density-functional theory (DFT) calculations, we computed the lattice parameters and bandgaps for these two chemistries, and we compare them with past studies from the literature in Table 1. Going to equimolar mixed-anion CaHf(O 0.5 S 0.5 ) 3 chemistry, we simply average the experimentally reported lattice parameters for CaHfO 3 and CaHfS 3 in the same space group as a starting point for further geometry optimization using DFT computations. Our theoretical investigation starts with an explicit enumeration of all possible configurations spanned by different relative arrangements of S and O atoms within a 20-atom supercell of CaHf(O 0.5 S 0.5 ) 3 . As both the oxide and sulfide parent perovskite chemistries stabilize in a Pnma space group, it is fair to assume that this size of supercell is enough to enumerate all low energy anion orderings and octahedral distortion patterns possible for the equimolar mixed oxysulfide chemistry and going to larger cells will not reveal any new orderings. To enumerate the symmetry-unique configurations exhaustively for the target composition, we rely on an algorithm for generating derivative structures 32,33 , using the implementation made available in enumlib 34 . Enumeration of all possible symmetry-unique configurations for the equimolar mixed-anion composition CaHf(O 0.5 S 0.5 ) 3 within this supercell leads to a total of 246 possible candidate structures. Further details of our DFT computations are provided in the Methods section towards the end of the manuscript. Figure 1 captures the DFT-computed relative energetics and Kohn-Sham bandgaps for the entire set of configurations. The relative energy is defined by the enthalpy of mixing with respect to the two end-point compositions (i.e., CaHfO 3 and CaHfS 3 perovskites) in units of eV per CaHf(O 0.5 S 0.5 ) 3 5-atom formula unit. It is interesting to note that the different anion orderings have a remarkable effect on the stability and bandgap. For a fixed composition, just by changing the relative arrangement of S and O atoms on the anion sublattice results in a large variation in the mixing enthalpies (17.6-1.135 eV) and bandgaps (1.52-2.75 eV). Note that, for this particular chemistry, it turns out that all of the mixing energies are positive, such that it is thermodynamically favorable (at least at 0 K, i.e., considering enthalpic contributions alone) for the oxysulfide to phase separate into the respective oxide and sulfide components. We later show that by systematically changing A and B-site chemistries, this mixing enthalpy can be made negative. However, first we are interested in understanding configuration-dependent relative trends across the CaHf(O 0.5 S 0.5 ) 3 data set. With this aim, natural next questions are: "What different factors are responsible for the observed trends and variations in the energetics and electronic structure?" and "Is it possible to develop a more intuitive understanding of the same at a microscopic level?" Relative stability of anion-ordered configurations In the quest of addressing the above questions, we first focus on the mixing enthalpies. A visualization of the relaxed geometries for the entire set of 246 configurations reveals that the moststable arrangement is a fac-type HfO 3 S 3 octahedron, presented by configurations shown in Fig. 1a, b (i.e., all S-Hf-S and O-Hf-O bonds are cis with a 90 ∘ interaction), followed by a mer-type HfO 3 S 3 octahedron (each anion type forms a "T" shape having two cis interactions, and one trans). The least-stable configuration is mixed HfO 4 S 2 /HfO 2 S 4 octahedra with the maximum number of trans interactions (cf., configurations shown in Fig. 1c, d). A deeper analysis reveals that there are at least three different factorsnamely, electronic, elastic strain, and electrostatic effects-that are responsible for energetically favoring the cis-configurations and penalizing the trans-configurations.
Electronic effects are known to play an important role in stabilizing a cis-configuration for strongly bonded ligands in octahedral complexes containing high valence d 0 transition metal ions. For instance, Mo(NR) 2þ 2 (with R = alkyl or aryl group) or MoO 2þ 2 complexes with a cis-configuration of ligands tend to maximize the metal(d π )-ligand(p π ) covalency 35,36 . In addition to octahedral complexes, covalency effects are known to favor the formation of cisoctahedra over the trans-arrangement in a number of oxynitride perovskites 20,37 . Similar effects also play a role in stabilization of configurations with a cis-arrangement of anions around the transition metal ion octahedra in oxysulfides. In fact, a partial density-of-states (DOS) analysis of representative cis-versus transconfigurations in the BX 6 octahedra (shown in Fig. 2a) results in very different valence band DOS for the p-states of O and S-with O pstates pushed significantly lower in the cis-configuration as compared with the trans-configuration, as depicted in Fig. 2bconfirming that the electronic effects indeed have an important role in the relative stabilization of the anion cis-configurations.
The role of elastic strain can be understood on the basis of simple and intuitive arguments by considering different arrangements of metal-anion-metal chains reached by the cis-and transanion configurations. As graphically illustrated in Fig. 2c, whereas a cis-configuration of the O and S ions in an octahedral Experimental bandgap from ref. 68 . f Bandgap using the hybrid Heyd-Scuseria-Ernzerhof functional (HSE06) from ref. 39 .  Fig. 2d). In the latter case, the two different chains will have different "relaxed" lengths, owing to the significantly different ionic radii of the O and S anions (1.40 and 1.84 Å, respectively) 38 . Therefore, to be able to accommodate these chains along a given lattice vector, the crystal periodicity would require the longer -S-Hf-S-Hf-Schains to be compressed and the shorter -O-Hf-O-Hf-Ochains to be stretched. As a result, the trans-configurations lead to an elastic strain along one or more directions in the mixed-anion crystal, and the effect eventually leads to the relative stabilization of the cisconfigurations.
In brief, alluded to above, in oxysulfides the two anionic species -the O 2− and S 2− ions-have a large size mismatch and therefore the elastic strain effects naturally become important. However, the two anions share the same nominal charge states and, at least for a fixed anion sublattice, rearranging the anions at a given composition cannot alter the net electrostatic energy of the system. An Ewald summation based on the nominal charge states of the constituents for a perfect (unrelaxed) perovskite crystal results in a constant electrostatic energy for the entire set of configurations. This might lead to an erroneous conclusion that, unlike oxynitrides (where the two anion species have different nominal charge states: O 2− and N 3− ), for oxysulfides containing isovalent O 2− and S 2− anions the electrostatic effects are not important. However, our analysis shows that relative energetics of the different configurations is directly governed by lowering of the overall electrostatic energy during geometry relaxation, as discussed below.
To understand the role played by electrostatics in stabilizing certain anion orderings in a mixed-anion oxysulfide, it is instructive to look at the A-and/or B-site cation ordered perovskites where such effects are well understood 9,10,12 . Figure  2e-g show local coordination environments of an anion X in different A-and B-site cation ordered configurations of a double perovskite. In each case, the anion is coordinated by two B-site cations and four A-site cations. In case of a B-site rocksalt-ordered perovskite A 2 [B 1 B 2 ]X 6 , the two different B-site cations are on opposite sides from one another, as shown in Fig. 2e. As a result of the asymmetric local configuration caused by the B-site sublattice, the anion may easily shift towards the smaller cation and away from the larger cation to not only relieve lattice strains arising from the size mismatch, but also minimize the overall electrostatic energy of the system during the course of relaxation. Figure 2f, g compare local anion configurations for two different A-site orderings, namely rocksalt and layer orderings, in a [A 1 A 2 ]B 2 X 6type double perovskite. Unlike the B-site ordered case, the rocksalt ordering of the A 1 and A 2 cations creates a symmetric local configuration for the anion whereas the layer ordering breaks this symmetry, allowing for a larger degree of relaxation for the anion. Similar arguments can also be used to rationalize energetics of simultaneous ordering at the A-and B-sites of a [A 1 A 2 ][B 1 B 2 ]X 6 double perovskite. Indeed, for ordered perovskite systems where electrostatics has a dominant role (i.e., compounds that demonstrate large size and charge disparity with respect to the ordered cation pair at a given cation sublattice), a rocksalt ordering at the B-sublattice is often accompanied by a layer ordering at the Asublattice.
The aforementioned symmetry arguments can be readily extended to an anion-ordered perovskite. Although cis-ordering of the anions leads to a highly asymmetric octahedron, a transordered octahedron is symmetric. As a result, an oxysulfide configuration with cis-ordering of the O 2− and S 2− anions around the B-site cations can lower its strain and electrostatic energies more easily than those with the trans-ordering of the anions. If the local-symmetry arguments above are valid, relative energetics of mixing enthalpies presented in Fig. 1 must exhibit a strong correlation with the corresponding lowering of the overall electrostatic energy during the course of geometry relaxation. To confirm this, we computed the Madelung energy for each DFToptimized configuration with respect to a common unrelaxed reference state. Indeed, the results presented in Fig. 2h exhibit a strong correlation between the two quantities, indicating the validity of the arguments discussed above to understand and rationalize the relative energetics of anion ordering in oxysulfide perovskites.
Configuration-dependent trends in the bandgap After rationalizing the trends in the mixing enthalpy as a function of different anion configurations, next we focus our attention on the bandgap. As can be seen from Fig. 1, there is no obvious correlation between the plotted mixing enthalpies and bandgaps for the different configurations. In fact, for different configurations with a fixed value of mixing enthalpy, the DFT-computed bandgap can vary up to >1.5 eV. To systematically understand the factors responsible for this variation, in Fig. 3a-d, we plot orbitaldecomposed DOS for four different configurations, all having a mixing enthalpy value of 0.60 ± 0.02 eV per formula unit, with significantly different bandgaps. In particular, the DOS associated with Ca s-states, Hf d-states and O and S p-states are shown in each panel. In agreement with previous reports in literature 39,40 , the valence band maximum and conduction band minimum in this class of ABX 3 compounds are primarily formed by anion pstates and B-site transition metal d-states. From the results plotted in Fig. 3, it can be seen that the bandgap variation across different configurations is largely controlled by the extent of lowering of the B-site transition metal ion d-states, relative to the valence band maximum. This lowering for a particular B-site ion is mainly governed by the crystal field imposed by a specific arrangement of different anions surrounding the octahedrally coordinated ion. The overall bandgap is then a collective effect of these local dband lowerings at the individual sites. Therefore, the bandgap trends in the target chemical space are expected to be determined by not only local arrangements of anion species around the B-site octahedra but also global arrangement of these octahedral motifs with respect to each other in a crystal.
Configurational machine learning Although the analysis presented thus far provides a physics-based qualitative understanding of relative energetics and bandgap variations as a function of anion configurations, next we focus on learning these trends in a quantitative manner using machine learning (ML). Further, a ML-based analysis allows us to quantify the impact of local relaxation effects on the two properties, by training the model on either the relaxed or unrelaxed geometries. Regardless of specific aims, materials ML strategies generally consist of two distinct steps. The first step is to come up with a physically meaningful numerical representation for each material in the data set. As a result of this step, each input material gets reduced to a string of numbers (frequently referred to as a fingerprint or descriptor). This is an enormously important step, requiring significant expertize and knowledge of the materials class and the application, i.e., "domain expertize". The second step establishes a mapping between the fingerprinted input and the target property, and is entirely numerical in nature, largely devoid of the need for domain knowledge.
For the problem at hand, as the relative arrangement of O 2− and S 2− species in different systems is deemed most relevant, we naturally resort to fingerprinting schemes that explicitly account for atomic configurations. Specifically, we employ (1) bag of bonds (BoB) fingerprint based on a sine Coulomb matrix (SCM) and (2) simulated X-ray diffraction (XRD) powder pattern spectra. The SCM represents a simple extension of the molecular Coulomb matrix to periodic crystals 41 . For the XRD spectra based fingerprint, we used a one-dimensional array representing powder diffraction pattern of a crystal structure 42 . A fixed 2θ range (from 0 ∘ to 127 ∘ ) with Cu K α radiation was consistently used for the entire set of configurations and a Gaussian broadening of 0.05 was used for smoothening the XRD spectra. To map the configurations on to the target properties, namely the mixing enthalpy and the bandgap, we employ kernel ridge regression (KRR) as a statistical learning algorithm [43][44][45] . A detailed description of both the . e-g Representative anion-site local environments in rocksalt-ordered A 2 B 1 B 2 X 6 -and A 1 A 2 B 2 X 6 -type, and layered-ordered A 1 A 2 B 2 X 6 -type double perovskites, respectively. h The mixing enthalpy of the oxysulfide configurations, which combines the effects of electronic effects, strain relaxation and local-symmetry-breaking, is well-correlated with the increase of electrostatic cohesive energy during relaxation. Figure 1. For the SCM-based BoB and XRD powder spectra descriptors implementations available in the matminer library 42 were used and KRR implementation with a Gaussian kernel available in scikitlearn was used for ML 46 .

employed descriptors and the KKR ML model is provided in Supplementary Discussion and Supplementary
Results of our ML model's training and prediction performance for the two targeted properties are presented in Table 2. In each case, the two configurational descriptors (SCM-based BoB and XRD powder spectra) are computed on either the DFT-relaxed ground-state geometry or the unrelaxed geometry where both lattice parameters and the internal atomic coordinates were fixed to ideal GdFeO 3 -type structure for all different configurations. A set of 90% of randomly selected data was used for the KRR ML models' training and the remaining 10% unseen data set was employed to test prediction performance of the trained models. To avoid overfitting, the optimal KRR model hyperparameters were selected during model training stage via minimizing a fivefold cross-validation root mean squared error (RMSE). To account for slight variability with respect to different training test data splits, in Table 2 the reported mean absolute errors (MAEs) and RMSEs are averaged over five different training/test splits.
It can be seen from the results presented in Table 2 that while reasonable estimates of the two properties on unseen data are possible using either of the two descriptors, descriptors computed on top of the relaxed geometries result in significantly better estimates of the target properties. This is not entirely unexpected and points towards a sensitive dependence of the target properties on the details of local structural relaxations. Although the ML models built on the descriptors computed using the relaxed geometries are more accurate, the models employing unrelaxed geometries are practically more useful as they are able to provide reasonably accurate trends of the properties without requiring an actual DFT computation. Furthermore, owing to systematic trends in the two properties across different A and Bsite chemistries (as discussed in detail in the following section), it is expected that such ML models can easily be extended beyond the specific CaHf(O 0.5 S 0.5 ) 3 chemistry employed here, with General trends across chemistries After developing an understanding of the trends in energetics and bandgaps as a function of anion configuration in CaHf(O 0.5 S 0.5 ) 3 perovskite, next we focus on how these trends might be affected in different (but systematically varying) cation chemistries. In particular, we consider nine different chemistries with A ∈ {Ca, Sr, Ba} and B ∈ {Ti, Zr, Hf}. Each of these nine chemistries were studied in two different structures, which are close in energy and exhibit the two lowest enthalpies of mixing for CaHf(O 0.5 S 0.5 ) 3 within the entire set of 246 configurations studied. Although the relaxed geometries for all the structures are provided as Supplementary Data 1, we note that both the structures show a cis-configuration of S and O atoms around HfO 6 octahedra and are closely related to each other. The only difference lies in the direction in which the zigzag -S-Hf-S-Hf-S-cis-bonded chains run. Henceforth these two configurations are referred to as Configs. A and B. A space group determination using FINDSYM 47 , with a tolerance of 0.01 Å for the internal atomic positions as well as the lattice, resulted in the Pmc2 1 space group for both the structures considered in the entire set of chemistries with A ∈ {Ca, Sr, Ba} and B ∈ {Ti, Zr, Hf}.
The computed mixing enthalpies and generalized gradient approximation (GGA)-PBE bandgaps across this chemical space for the two configurations are presented in Fig. 4a, b, respectively. Note that for a given chemistry, both configurations always show very similar mixing enthalpies and bandgaps. Further, it can be seen from the figure that there are systematic trends going from one cation chemistry to the other, on both the A-and B-site sublattices, which suggests that the analysis made thus far for the specific CaHf(O 0.5 S 0.5 ) 3 chemistry is likely to extrapolate to the other related cation chemistries.
Besides the systematic variation of the two properties across chemistries, Figs. 4a, b reveal several important points. First, the ordering of chemical trends in the two properties differ with respect to the B-site chemistries. Although the mixing enthalpies decrease in an order of Ti > Hf > Zr, the bandgaps generally increase in an order of Ti < Zr < Hf, irrespective of the A-site chemistry. In other words, zirconates and hafnates have different ordering for the two properties. This can be rationalized by noting the fact that, for different B-site cation chemistries, the size effects dictate the variation in mixing enthalpies, whereas the electronegativity of the cation plays a dominant role in determining the bandgap. Although the B cation electronegativities follow the expected trend and decrease going down the Periodic Table from Ti to Hf, the Shannon's ionic radii (for a VI-fold coordinated tetravalent oxidation state) of the cations increase in an order of Ti (0.605 Å) < Hf (0.71 Å) < Zr (0.72 Å). The different ionic radii ordering is due to the well-known Lanthanide contraction effect -as a result of the filling of the 4f shell in Hf, an increased nuclear charge leads to a much stronger electrostatic interaction between the nuclei and the electron cloud resulting in a slightly smaller ionic radii for Hf (0.71 Å) as compared with that of Zr (0.72 Å). Second, the DFT-computed bandgap for BaTi(O 0.5 S 0.5 ) 3 is anomalously large, as can be seen from Fig. 4b. A closer look at the relaxed structure reveals that, despite being dynamically stable (i.e., no soft mode phonon instabilities), a combination of a large A-site cation and a small B-site cation leads to extremely stretched Ti-S bonds in this system. As a consequence of this local strain, the bandgap further opens up, breaking the common trend exhibited by the remaining chemistries. This local strain also leads to a relatively higher positive enthalpy of mixing for BaTi(O 0.5 S 0.5 ) 3 . Based on our results, only zirconates of Ca and Sr are predicted to be thermodynamically stable based on the mixing enthalpies. Finally, we note that we have also confirmed that the chemical   Figs. 5c, d, respectively). Given that the anion p-states largely contribute to the valence band maximum, along the axis of anion ordering (i.e., along the Γ-X and Γ-Y paths), the valence band dispersion resembles that of the pure sulfide perovskite, whereas in the other two directions, oxygen blocks electron hopping and the bands formed present details that are characteristic of the oxide. The band structure of SrZr (O 0.5 S 0.5 ) 3 exhibits a qualitatively similar behavior.
To understand why only Ca and Sr zirconates exhibit negative mixing enthalpy and to further investigate thermodynamic formability in these oxysulfide chemistries, we resort to an empirical analysis based on classical structure factors. The problem of perovskite formability has traditionally been addressed using a pair of key geometric factors, namely the Goldschmidt tolerance factor (t f ) and the octahedral factor (o f ). Although the former directly captures the ability of a given ABX 3 -type perovskite chemistry to achieve a packing arrangement of different atoms in alternate AX and BX 2 {001} family of planes as close to that of an ideal cubic perovskite crystal structure, the latter is only concerned with the stability window of the B-site octahedral coordination.
The conventional t f and o f are defined as: ðr A þ r X Þ= ffiffi ffi 2 p ðr B þ r X Þ and r B /r X , respectively, where r i represents Shannon's coordination-dependent ionic radius for the atom i. A closer look at the cis-configured AB(O 0.5 S 0.5 ) 3 structure, shown in Fig. 6a, readily points out the problem associated with the employed definition of the compositionally averaged modified t f . From Fig. 6a, it can be seen that for the B(O 0.5 S 0.5 ) 2 {001} planes, the alternating arrangement of O and S does justify use of an averaged anion ionic radius in both the o f and denominator of the t f . However, the {100} planes with A-site cations exhibit atomic layers containing either 100% O or 100% S as anionic species. Therefore, given the much larger atomic size of S as compared with O, only AS planes can be close packed. This insight indicates that the modified tolerance factor should rather be defined as ðr A þ r S Þ= ffiffi ffi 2 p ½r B þ 0:5ðr O þ r S Þ for the present case. Using this revised modified t f and an anion-averaged o f , a structure map presented in Fig. 6b correctly explains the DFT-computed relative stabilities across the considered set of chemistries. In particular, titanates are unstable because of a much smaller B-site cation radius (average o f < 0.414) and Ba-containing chemistries are unfavorable because A-site radius is larger than desired (modified t f > 1.04). The figure also shows that the energetics of the oxysulfides is very sensitive to the size of B-site cations (as compared with that of the A-site cations) and a very small size difference between the Zr and Hf is actually sufficient to switch the trends of thermodynamic stabilities between some of these chemistries. At last, we note that, in addition to being helpful in rationalization of stability trends, the structure map presented in Fig. 6b can, in principle, be used to screen other potentially stable equimolar cis-ordered oxysulfide chemistries. If a given chemistry falls near and above the highlighted region (shown by a dashed rectangle) of stable zirconates, the geometric arguments favor stability in a perovskite crystal structure, though many other factors may also need to be addressed before any firm conclusions can be made.
At last, we note that based on the computed mixing enthalpies both SrZr(O 0.5 S 0.5 ) 3 and CaZr(O 0.5 S 0.5 ) 3 are predicted to be thermodynamically stable. Since synthesis of BaZr(O 0.5 S 0.5 ) 3 has already been demonstrated in the past 18 , we anticipate that experimental realization of relatively stable (although thermodynamically metastable) 50 CaHf(O 0.5 S 0.5 ) 3 and SrHf(O 0.5 S 0.5 ) 3 should also be highly likely. Anion-order-driven ferroelectricity This work not only identifies and analyzes the factors that control anion ordering in perovskite oxysulfides, but also has considerable implications towards functionality exhibited by heteroanionic perovskites. For instance, considering the two stable zirconate chemistries shown in Fig. 4a, we note that while the two end members AZrO 3 and AZrS 3 (with A=Ca and Sr) are both centrosymmetric with a Pnma space group, the most-stable cistype anion ordering breaks mirror symmetry along the c axis owing to the alternating AO and AS planes, lowering the symmetry to a polar Pmc2 1 space group. Therefore, the cisordering could also potentially lead to the phenomena of hybrid improper ferroelectricity [51][52][53][54] , if the preferred anion ordering can manifest in large single domains. To examine the potential ferroelectricity in these polar materials, we compute the macroscopic polarization and the associated switching barrier using the Berry phase approach 55 made available within the framework of the modern theory of polarization [55][56][57] . As a non-polar reference, we take the centrosymmetric Pmma structure obtained by undoing the octahedral rotation around the c axis using PSEUDO 58 , and we relax this with constrained symmetry. To fit the polarization branch 55 , we use ISODISTORT 59 to interpolate three intermediate images along the switching pathway between Pmma and Pmc2 1 . The calculated polarizations and switching barriers for the two stable zirconate chemistries are shown in Fig.  7a and the structural details of the Pmc2 1 and Pmma phases are depicted in Figs. 7b, c, respectively. We find that both chemistries are predicted to be ferroelectric with a computed macroscopic polarization of >30 μC/cm 2 . Further, quantitative details for our computations are provided in Table 3.
Although the ferroelectric mechanism here is similar to that of the known hybrid improper ferroelectrics 52-54 , here we note an important distinction. In a Pnma ABX 3 -type perovskite, the A-site cations in each {001} plane along the unique axis of tilting displace in alternating directions. The perovskite is non-polar because these dipoles cancel exactly. However, typically when symmetry is lowered to the Pmc2 1 space group by breaking the mirror symmetry along the c axis (for instance in a double perovskite containing two A-site cations arranged in a layered-type ordering) 54 , one of the dipoles becomes larger than the other, resulting in a relatively small net polarization. However, in these ordered oxysulfides, as a consequence of the cis-type anion ordering, all the A-site cations get displace in the same direction, resulting in a significantly larger macroscopic polarization than that of typical improper ferroelectrics reported in the literature 54 . Specifically, the cations in the AS layers behave as expected, whereas those in the AO layer reverse direction. This anion-order-driven polarization is  Owing to the much larger size of the S atom, only AS planes can be close-packed, forming the basis for a modified tolerance factor t f proposed in this work. b The proposed modified tolerance factor t f and anion-averaged octahedral factor o f readily cluster the two stable zirconate chemistries (shown within a dashed rectangle) and allow us to rationalize the DFT-computed trends in relative energetics. Note that the Cartesian axes in the figure refer to a five-atom cubic perovskite unit cell.
predicted to result in a macroscopic polarization that is larger than that of the well-known ferroelectrics such as BaTiO 3 60 and KNbO 3

DISCUSSION
Oxysulfide perovskite compounds represent an emerging class of materials with a potential use in a diverse range of technological applications, including solar energy harvesting, catalysis, and sensing. In contrast to heterocationic double perovskite compounds, where the underlying design rules and cation ordering tendencies are relatively well understood and documented, the factors that govern the anion ordering phenomena within this class of heteroanionic compounds remain poorly understood. As an important step in this direction, the present manuscript employed a selected set of first-principles-based DFT computations to first study a wide range of anion orderings in a prototypical equimolar oxysulfide perovskite CaHf(O 0.5 S 0.5 ) 3 chemistry. Subsequently, the most-stable anion ordering was studied and compared in a range of different AB(O 0.5 S 0.5 ) 3 chemistries with A ∈ {Ca, Sr, Ba} and B ∈ {Ti, Zr, Hf}. The salient findings of our work can be summarized as follows: • AB(O 0.5 S 0.5 ) 3 oxysulfide perovskites with a B-site cation in a d 0 electronic configuration exhibit a strong preference for a cistype anion ordering. The driving force for this specific ordering tendency stems from a combined interplay of electronic, elastic, and electrostatic factors.
• Descriptors based on either BoB-SCM or simulated XRD spectra are shown to be effective in learning the anionconfiguration-dependent variations in the mixing enthalpies and the bandgaps, with both descriptors leading to a comparable performance.  Table, mixing enthalpies exhibit a relatively more complex behavior.

•
The trends in mixing enthalpies of these oxysulfides (computed relative to the respective pure oxide and sulfide chemistries) can be understood by using the compositionally averaged octahedral factor and a modification of the classical tolerance factor that accounts for presence of AS planes, as the larger sulfide ions are the limiting factor in forming a close-packed structure. The results for both the config. A and config. B are shown for each of the chemistries. Net polarization direction refers to the Cartesian axes in a GdFeO 3 -type Pnma crystal structure.