Anion order in oxysulfide perovskites: origins and implications

Abstract

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(O_0.5S_0.5)₃ 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 HfO₃S₃ 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/cm², exhibiting a significant promise for electronic materials applications.

Introduction

Perovskite compounds—with a generic formula ABX₃—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.

Mixed-anion compounds (i.e., containing more than one anion species) add new dimensions to control and tune stability as well as different aspects of functionality, owing to the different elemental characteristics of the anions, such as charge, ionic radii, electronegativity, and polarizability, as well as their relative arrangement in the crystal lattice⁷. Although mixed cation perovskites have been widely studied^{8,9,10,11,12,13}, mixing on the anion sublattice in perovskites remains largely unexplored, with the exception of a few recent studies^{14,15,16,17,18,19,20,21,22}.

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₃ compounds within this class—with A=Ca, Sr, or Ba; B=Ti, Zr, or Hf; X=S or Se—are 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²³. For some ABX₃ sulfides, competing non-perovskite phases have also been reported^{24,27,28,29,30}. More recently, a wide range of ABS₃ perovskites and non-perovskites were screened for relevant properties, such as thermodynamic stability, light absorption, charge mobility, and defect tolerance. The study identified LaYS₃ as a promising candidate³¹. As a first step in the exploration of perovskite oxysulfides, Perera et al¹⁸. recently synthesized a series of BaZr(O_1−xS_x)₃ 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.5S_0.5)₃ 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.5S_0.5)₃, 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₆ octahedron is highly preferred over a trans-arrangement. 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₆ octahedron, which primarily changes the position of the Hf d-states. We then generalize this insight by exploring stability trends for the larger set of AB(O_0.5S_0.5)₃ 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₃, CaHfS₃, and CaHf(O_0.5S_0.5)₃ chemistries

Both CaHfO₃ and CaHfS₃ are known to be stable in a GdFeO₃-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.5S_0.5)₃ chemistry, we simply average the experimentally reported lattice parameters for CaHfO₃ and CaHfS₃ 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.5S_0.5)₃. 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³⁴. Enumeration of all possible symmetry-unique configurations for the equimolar mixed-anion composition CaHf(O_0.5S_0.5)₃ 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.

Table 1 Lattice parameters and bandgaps of the end-point compositions, CaHfO₃ and CaHfS₃ perovskites, forming the mixed CaHf(O_0.5S_0.5)₃ equimolar oxysulfide chemistries studied here.

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₃ and CaHfS₃ perovskites) in units of eV per CaHf(O_0.5S_0.5)₃ 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.5S_0.5)₃ 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?”

Fig. 1: Anion ordering dependent energetics and bandgaps.

Configuration-dependent mixing enthalpies and bandgaps of 246 different anion-ordered CaHf(O_0.5S_0.5)₃ configurations possible within a 20-atom $\sqrt{2}\times \sqrt{2}\times 2$ pseudo-cubic perovskite supercell. Representative octahedral configurations of anions around Hf are highlighted in a few selected most-stable a, b and least-stable c, d supercell geometries.

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 most-stable arrangement is a fac-type HfO₃S₃ 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₃S₃ octahedron (each anion type forms a "T” shape having two cis interactions, and one trans). The least-stable configuration is mixed HfO₄S₂/HfO₂S₄ 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 factors—namely, 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⁰ 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 cis-octahedra 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 trans-configurations in the BX₆ octahedra (shown in Fig. 2a) results in very different valence band DOS for the p-states of O and S—with O p-states pushed significantly lower in the cis-configuration as compared with the trans-configuration, as depicted in Fig. 2b—confirming that the electronic effects indeed have an important role in the relative stabilization of the anion cis-configurations.

Fig. 2: Factors dicatating anion ordering in oxysulfides.

A graphical representation of the electronic, elastic strain, and electrostatic effects that stabilize cis- ordering while penalizing the trans- order in oxysulfide perovskites. a A schematic representing electronic stabilization of cis-configuration for strongly bonded ligands in octahedral complexes containing d⁰ transition metal ions, as compared with the trans-configuration. b Partial density-of-states (DOS) plots showing that the cis-configured octahedra in a have oxygen p-states at lower energies. The cis- and trans-configurations of anions also lead to different anion-orderings in anion-metal-anion chains. For instance, in a cis-configured CaHf(O_0.5S_0.5)₃ structure c all the chains have the same –S–Hf–O–Hf–S–Hf–O– arrangement, as shown in top panel, which can fit in a crystal without introducing any additional lattice strain (bottom panel). However, the trans-configuration shown in d has mismatched chain lengths, as O and S have significantly different radii (top panel), which have to be either stretched or compressed to fit in the periodic boundary conditions posed by a crystal lattice (bottom panel). e–g Representative anion-site local environments in rocksalt-ordered A₂B¹B²X₆- and A¹A²B₂X₆-type, and layered-ordered A¹A²B₂X₆-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.

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 trans-anion configurations. As graphically illustrated in Fig. 2c, whereas a cis-configuration of the O and S ions in an octahedral environment leads to equi-length zigzag chains of only one type (i.e.,–O–Hf–S–Hf–O–) in all the three directions within a crystal, a trans-configuration can result in both –O–Hf–O–Hf–O– and –S–Hf–S–Hf–S– chains running parallel to each other in a given crystal direction (cf. 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)³⁸. Therefore, to be able to accommodate these chains along a given lattice vector, the crystal periodicity would require the longer –S–Hf–S–Hf–S– chains to be compressed and the shorter –O–Hf–O–Hf–O– chains 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 cis-configurations.

In brief, alluded to above, in oxysulfides the two anionic species—the O²⁻ and S²⁻ 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²⁻ and N³⁻), for oxysulfides containing isovalent O²⁻ and S²⁻ 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₂[B¹B²]X₆, 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¹A²]B₂X₆-type double perovskite. Unlike the B-site ordered case, the rocksalt ordering of the A¹ and A² 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¹A²][B¹B²]X₆ 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 A-sublattice.

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 trans-ordered octahedron is symmetric. As a result, an oxysulfide configuration with cis-ordering of the O²⁻ and S²⁻ 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 DFT-optimized 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 orbital-decomposed 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₃ compounds are primarily formed by anion p-states 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 d-band 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.

Fig. 3: Anion ordering dependent electronic structure.

Orbital-decomposed partial DOS for four different CaHf(O_0.5S_0.5)₃ configurations a–d, selected to have a similar range of the mixing enthalpies (within 0.60 ± 0.02 eV per formula unit), but with significantly different GGA-PBE bandgaps. Ca s-states, O and S p-states, and Hf d-states are shown in blue, red, black, and green, respectively.

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²⁻ and S²⁻ 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⁴¹. For the XRD spectra based fingerprint, we used a one-dimensional array representing powder diffraction pattern of a crystal structure⁴². 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 employed descriptors and the KKR ML model is provided in Supplementary Discussion and Supplementary Figure 1. For the SCM-based BoB and XRD powder spectra descriptors implementations available in the matminer library⁴² were used and KRR implementation with a Gaussian kernel available in scikit-learn was used for ML⁴⁶.

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₃-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.

Table 2 Performance of KRR ML models built on different sets of configurational fingerprints or descriptors in predicting the computed mixing enthalpies and DFT-PBE bandgaps of various CaHf(O_0.5S_0.5)₃ configurations considered in the present study.

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 B-site 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.5S_0.5)₃ chemistry employed here, with exceedingly lower amounts of required training data as new chemistries are added.

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.5S_0.5)₃ 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.5S_0.5)₃ 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₆ 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⁴⁷, with a tolerance of 0.01 Å for the internal atomic positions as well as the lattice, resulted in the Pmc2₁ 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.5S_0.5)₃ chemistry is likely to extrapolate to the other related cation chemistries.

Fig. 4: Relative chemical trends across chemistries.

Systematic trends in the computed mixing enthalpies and DFT bandgaps in AB(O_0.5S_0.5)₃ chemistries with A ∈ {Ca, Sr, Ba} and B ∈ {Ti, Zr, Hf}. Each chemistry is studied in two different cis-ordered configurations (both with a Pmc2₁ space group) identified as the lowest-energy configurations for CaHf(O_0.5S_0.5)₃ from within the entire set of 246 configurations studied and shown in Figs. 1a, b.

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.5S_0.5)₃ 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.5S_0.5)₃. 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 trends in mixing enthalpies and GGA-PBE bandgaps identified in Figs. 4a, b, respectively, using the GGA-PBE functional do not change when using higher level meta-GGA and hybrid functionals. In particular, for a selected set of cation orderings, we recomputed the mixing enthalpies and GGA-PBE bandgaps using Strongly Constrained and Appropriately Normed (SCAN)⁴⁸ and hybrid Heyd–Scuseria-Ernzerhof HSE06 functionals⁴⁹, respectively, for each of the nine chemistries shown in Fig. 4. Our results, presented in the Supplementary Discussion, Supplementary Figs. 2 and 3 indicate that these functionals possibly improve the quantitative nature of our results, albeit at a significantly higher computational cost, but do not change any of the physical trends or conclusions.

Comparing the electronic band structure of the two configurations (configs. A and B) in each of these stable chemistries reveals that details of anisotropic distributions of anions are reflected in the band structures. For instance, as shown in Figs. 5a, b for CaZr(O_0.5S_0.5)₃, config. A showcases continuous zigzag –S–Zr–S–Zr– chains running along the y axis, whereas in config. B, these chains run along the x axis (for comparison, the band structures of CaZrO₃ and CaZrS₃ are shown in 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.5S_0.5)₃ exhibits a qualitatively similar behavior.

Fig. 5: Electronic band structures of homo and heteroanion chemistries.

DFT-computed band structures for a CaZr(O_0.5S_0.5)₃ config. A (space group Pmc2₁), b CaZr(O_0.5S_0.5)₃ config. B (space group Pmc2₁), c CaZrO₃ (space group Pnma), and d CaZrS₃ (space group Pnma). The two configurations for the mixed heteroanion chemistry differ in the manner the –S–Zr–S–Zr– zigzag chains run along the supercell, as identified by black circles in the first two panels.

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₃-type perovskite chemistry to achieve a packing arrangement of different atoms in alternate AX and BX₂ {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})/\sqrt{2}({r}_{B}+{r}_{X})$ and r_B/r_X, respectively, where r_i represents Shannon’s coordination-dependent ionic radius for the atom i. Based on simple geometrical arguments, t_f should be as close to unity as possible and o_f should fall within [$\sqrt{2}$−1, $\sqrt{3}$−1] to form a stable cubic perovskite. Extensions of t_f and o_f to a mixed-anion chemistry AB(O_0.5S_0.5)₃ by taking compositional averages of the anion ionic radii as $[{r}_{A}+0.5({r}_{O}+{r}_{S})]/[\sqrt{2}\left({r}_{B}+0.5({r}_{O}+{r}_{S})\right.]$ and r_B/0.5(r_O + r_S), however, fail to capture the stability trends presented in Fig. 4a. In particular, such a definition of a modified t_f results in a metric which is near the formability borderline (on the lower side) for CaZr(O_0.5S_0.5)₃ and SrZr(O_0.5S_0.5)₃, whereas BaZr(O_0.5S_0.5)₃ is predicted to lie in a highly favorable region of the structure map.

A closer look at the cis-configured AB(O_0.5S_0.5)₃ 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.5S_0.5)₂ {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})/\sqrt{2}[{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.

Fig. 6: Geometric factors influencing formability in a cis-ordered oxysulfide.

a A schematic showing local atomic arrangements in a cis-ordered AB(O_0.5S_0.5)₃ perovskite oxysulfide. The AS and BOS close-packed planes along the [010] direction are explicitly shown. Note that alternate {001} AX atomic planes feature either all-oxide or all-sulfide stackings. 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.

At last, we note that based on the computed mixing enthalpies both SrZr(O_0.5S_0.5)₃ and CaZr(O_0.5S_0.5)₃ are predicted to be thermodynamically stable. Since synthesis of BaZr(O_0.5S_0.5)₃ has already been demonstrated in the past¹⁸, we anticipate that experimental realization of relatively stable (although thermodynamically metastable)⁵⁰ CaHf(O_0.5S_0.5)₃ and SrHf(O_0.5S_0.5)₃ 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₃ and AZrS₃ (with A=Ca and Sr) are both centrosymmetric with a Pnma space group, the most-stable cis-type anion ordering breaks mirror symmetry along the c axis owing to the alternating AO and AS planes, lowering the symmetry to a polar Pmc2₁ space group. Therefore, the cis-ordering 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⁵⁵ 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⁵⁸, and we relax this with constrained symmetry. To fit the polarization branch⁵⁵, we use ISODISTORT⁵⁹ to interpolate three intermediate images along the switching pathway between Pmma and Pmc2₁. The calculated polarizations and switching barriers for the two stable zirconate chemistries are shown in Fig. 7a and the structural details of the Pmc2₁ 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². Further, quantitative details for our computations are provided in Table 3.

Fig. 7: Anion ordering induced ferroelectric polarization.

a The DFT-computed barrier for ferroelectric switching in the two thermodynamically stable zirconate chemistries. The results for both config. A and config. B are shown in each case. b, c present the top and side views of the Pmc2₁ ground state and Pmma saddle point geometric structures in the config. A.

Table 3 The DFT-computed macroscopic polarizations and ferroelectric switching barriers in the two thermodynamically stable zirconate chemistries.

Although the ferroelectric mechanism here is similar to that of the known hybrid improper ferroelectrics^52,53,54, here we note an important distinction. In a Pnma ABX₃-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₁ 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)⁵⁴, 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⁵⁴. Specifically, the cations in the AS layers behave as expected, whereas those in the AO layer reverse direction. This anion-order-driven polarization is predicted to result in a macroscopic polarization that is larger than that of the well-known ferroelectrics such as BaTiO₃⁶⁰ and KNbO₃ ⁶¹.

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.5S_0.5)₃ chemistry. Subsequently, the most-stable anion ordering was studied and compared in a range of different AB(O_0.5S_0.5)₃ 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.5S_0.5)₃ oxysulfide perovskites with a B-site cation in a d⁰ electronic configuration exhibit a strong preference for a cis-type 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 anion - configuration-dependent variations in the mixing enthalpies and the bandgaps, with both descriptors leading to a comparable performance.

• Heteroanionic cis-ordered AB(O_0.5S_0.5)₃ chemistries with A ∈ {Ca, Sr, Ba} and B ∈ {Ti, Zr, Hf} exhibit well-behaved trends in both mixing enthalpies and bandgaps. While the bandgap variations are well aligned with trends in the metallic nature of the A- and B-site cations in the Periodic 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.

• Based on the computed mixing enthalpies, we predict a relative stability ordering of the following oxysulfide perovskites as: SrZr(O_0.5S_0.5)₃ > CaZr(O_0.5S_0.5)₃ > CaHf(O_0.5S_0.5)₃ > SrHf(O_0.5S_0.5)₃ > BaZr(O_0.5S_0.5)₃, with a negative enthalpy of mixing for the first two chemistries. Given that a number of BaZr(O_1−xS_x)₃ perovskites spanning over the wide range of x ∈ [0, 0.54] have already been synthesized and reported¹⁸, all of the above chemistries should be amenable to experimental synthesis. In particular, zirconates of Ca and Sr should exhibit long term stability.

• The strongly favored cis-type anion ordering breaks inversion symmetry locally, and can potentially lead to a large macroscopic polarization >30 μC/cm², if large anion-ordered domains are sustained.

Methods

Details of DFT computations

The DFT calculations were performed using the Vienna Ab initio Simulation Package (VASP)^62,63 employing the Perdew, Burke, and Ernzerhof (PBE)⁶⁴ GGA exchange-correlation functional. The electronic wave functions were expanded in plane waves up to a cutoff energy of 500 eV. We used potentials based on the projector-augmented wave method⁶⁵ explicitly including the following valence electronic configurations for different elemental species; Ca: 3s²3p⁶4s², Sr: 4s²4p⁶5s², Ba: 5s²5p⁶6s², Zr: 4p⁶4d²5s², Hf: 5p⁶5d²6s², O: 2s²2p⁴, and S: 3s²3p⁴. A Gamma-centered 5 × 5 × 3 Monkhorst–Pack k-point mesh⁶⁶ was used for Brillouin-zone integrations for an orthorhombic supercell with Pnma space group containing 20 atoms (i.e., a $\sqrt{2}\times \sqrt{2}\times 2$ pseudo-cubic supercell). Spin-unpolarized calculations were employed as none of the simulated systems contained any unpaired electrons. To obtain a geometry optimized equilibrium structure, atomic positions and the supercell lattice parameters were fully relaxed using the conjugate-gradient method until all the Hellmann–Feynman forces and the stress components were >0.05 V/Å and 1.0 × 10⁻² GPa, respectively.