Prediction of structure and cation ordering in an ordered normal-inverse double spinel

Abstract

Spinels represent an important class of technologically relevant materials, used in diverse applications ranging from dielectrics, sensors and energy materials. While solid solutions combining two “single spinels” have been explored in a number of past studies, no ordered “double” spinels have been reported. Based on our first principles computations, here we predict the existence of such a double spinel compound MgAlGaO₄, formed by an equimolar mixing of MgAl₂O₄ normal and MgGa₂O₄ inverse spinels. After studying the details of its atomic and electronic structure, we use a cluster expansion based effective Hamiltonian approach with Monte Carlo simulations to study the thermodynamic behavior and cation distribution as a function of temperature. Our simulations provide strong evidence for short-ranged cation order in the double spinel structure, even at significantly elevated temperatures. Finally, an attempt was made to synthesize the predicted double spinel compound. Energy Dispersive X-ray Spectrometry and X-ray diffraction Rietveld refinements were performed to characterize the single-phase chemical composition and local configurational environments, which showed a favorable agreement with the theoretical predictions. These findings suggest that a much larger number of compounds can potentially be realized within this chemical space, opening new avenues for the design of spinel-structured materials with tailored functionality.

Introduction

Spinel compounds—typically represented as AB₂O₄, with A and B as metal ions—form an important class of materials^1,2,3,4,5. Spinels are found as naturally occurring minerals all over the globe⁶ and have traditionally been utilized as precious stones, such as red or blue spinels containing either Cr⁺³ or Fe⁺²/Zn⁺² ions, respectively^7,8. The fact that nearly all of the main group and transition metals can be synthesized in a stable spinel structure makes it a relatively large family of compounds with manifold compositions, cation ordering, electron configurations, and valence states⁹. Owing to this large synthetic flexibility and a wide range of atomic and electronic configurations accessible within this chemical space, spinels exhibit interesting magnetic^5,10,11, electronic^12,13, optical¹⁴, and catalytic^15,16 properties useful for a diverse range of applications, including data storage¹⁷, high frequency electronic devices¹⁸, dielectrics¹⁹, transparent conducting oxides^20,21,22,23, lasers^24,25, sensing², energy storage^26,27, superconductivity²⁸, and biotechnology²⁹.

One of the distinguished features of the spinels is the wide range of cation distributions accessible in this system. The “normal” spinal structure was characterized in 1915 by Bragg and Nishikawa, who confirmed that most commonly occurring spinels generally have the A atoms tetrahedrally coordinated, while the B atoms reside in an octahedral coordination. In this structure, $\frac{1}{8}$ of the tetrahedral and $\frac{1}{2}$ of the octahedral voids are occupied by the A and B metal atoms, respectively in a face-centered-cubic close-packed oxygen sublattice. Thus, a normal spinel can be represented by [A]^t[B₂]^oO₄, where the superscripts t and o label the elemental species A and B occupying the tetrahedral and octahedral sites, respectively. In 1935, Barth and Posnjak showed that not all spinels have the normal structure as their ground state configuration and there exist several chemistries with the “inverse” spinel configuration where the tetrahedral sites are occupied by the B atoms and the octahedral sites are shared equally by both A and B atoms, i.e., [B]^t[AB]^oO₄³⁰. While the 0 K ground state configuration for any given spinel is always either the perfectly ordered normal or inverse structure, at a finite temperature a mixing of elemental species within the octahedral lattice or across the octahedral and tetrahedral lattices is often observed, leading to [A_1−xB_x]^t[A_xB_2−x]^oO₄ configurations, with the subscript x characterizing the degree of inversion (often referred to as i in the literature). The inversion parameter can vary from 0 (for a normal spinel) to 1 (for an inverse spinel) and adopts a value of $\frac{2}{3}$ for a completely random distribution of the metal atoms. Prototypical examples of normal and inverse spinels are MgAl₂O₄ (space group $Fd\bar{3}m$) and MgGa₂O₄/MgIn₂O₄ (space group P4₁22), respectively (cf. Fig. 1). While MgIn₂O₄ exhibits a high degree of inversion at room temperature, reported structures for MgGa₂O₄ scatter over a wide range of intermediate degrees of inversion^30,31.

Fig. 1: Normal and inverse spinel structures.

a MgAl₂O₄ normal and b MgGa₂O₄ inverse ground state atomic configurations. In each case the unit cell is shown by solid black lines. Octahedral and tetrahedral atomic coordination environments are also identified by the coordination polyhedra in each case.

The problem of (a) identifying factors that govern the intrinsic tendency of a given binary (i.e., A₃O₄, with A = B and multiple possible charge states of A) or ternary (i.e., AB₂O₄) spinel to adopt either the normal or the inverse configuration and (b) characterizing its cation distribution at a finite temperature has been an active focus of materials research over the last half century^{32,33,34,35,36,37,38,39,40}. In addition to the single spinels, binary solid solutions formed by mixing two single spinels have also been studied^41,42,43. For the binary solid solutions, it is generally considered that normal and inverse single spinel chemistries typically tend not to mix. For instance, a dual-layer oxide formed due to high-temperature corrosion of ferritic-martensitic steels has been attributed to the formation of an inverse magnetite spinel (Fe₃O₄) along with an iron-chromium normal spinel of composition FeCr₂O₄, which tend to phase separate⁴⁴.

Note that a number of binary spinel solid solutions with a considerable miscibility at finite temperatures have been known for a long time now^41,42,43. For instance, cation mixing in high temperature (in 973–1473 K range) solid solutions of the MgAl_2−xGa_xO₄ quaternary system (with x = 0.38, 0.76, 0.96, 1.52) have been investigated by Ito et al.⁴²; however, no low temperature ordered double spinel structures formed by combining two single spinels (in this case, MgAlO₄ and MgGaO₄) have been reported until now. Given that such ordered double structures are well known to exist in perovskite and perovskite-related crystal systems^{45,46,47,48,49} with strong theoretical evidence for their existence in other more complex structures⁵⁰, it is natural to ask if, analogous to perovskites, double spinel chemistries exist. Further, if there is such a possibility, what specific type(s) of ordering(s) is(are) preferred by cations shared by the two parent single spinels.

In this contribution, we use density functional theory (DFT)^51,52 computations in combination with a cluster expansion^53,54 based effective Hamiltonian approach to show that, in contrast to the general chemical intuition-based understanding, ordered double spinel ground states formed by combining a normal and an inverse spinel can indeed be possible at low temperatures in specific spinel chemistries. Towards this end, we explore all quaternary double spinels possible via pairwise combining the three single spinels MgAl₂O₄, MgGa₂O₄, and MgIn₂O₄ over the entire range of compositions and a wide range of ordered cation configurations (nearly 350 and 160,000 configurations explored with DFT and cluster expansion, respectively). Our results show that, while the MgAl₂O₄–MgIn₂O₄ and MgGa₂O₄–MgIn₂O₄ binary systems do not like to mix at low temperatures (i.e., do not have any ordered configuration with a negative enthalpy of mixing), there exists a specific equimolar ordered configuration for the MgAl₂O₄–MgGa₂O₄ (more specifically [Ga]^t[MgAl]^oO₄) system that is distinctly favored by a negative enthalpy of mixing. After studying the details of the underlying electronic structure and cation distribution as a function of temperature for this newly-found double spinel structure, we synthesize the compound via a conventional ball-milling route and measure X-ray diffraction (XRD) patterns on samples after annealing at 200 °C and 500 °C over a range of annealing times. A direct comparison of the measured XRD patterns with those simulated for the ordered double spinel structure and a random distribution of the cations presents strong evidence for the existence of the identified double spinel ordering. Rather than being associated with the specific binary spinels studied here, if the discovered novel double spinel ground state configuration turns out to be a general phenomenon for the entire set of known spinels, the known spinel chemical space can be extended many-fold, thereby opening new avenues to rationally design and tune material functionality for a plethora of applications^{55,56,57,58,59,60}.

Results

DFT and cluster expansion-based effective Hamiltonians

We start by looking at ground states of bulk MgM₂O₄ (with M = Al, Ga, and In) spinel chemistries. The perfectly ordered normal and inverse spinel configurations for each of the three chemistries were simulated in cubic $Fd\bar{3}m$ and tetragonal P4₁22 space groups. In agreement with previous theoretical and experimental studies, we confirm that MgAl₂O₄ has a normal spinel ground state, while both MgGa₂O₄ and MgIn₂O₄ prefer an inverse spinel ground state. Table 1 compares the computed 0 K ground state symmetry and lattice parameters for these spinels with previously reported computational (0 K) and experimental results. Our DFT-computed ΔE = (E_Inverse − E_Normal) values of 0.192, −0.077, and −0.092 eV per formula unit for MgAl₂O₄, MgGa₂O₄ and MgIn₂O₄, respectively, are again in good agreement with those reported in previous calculations^62,63. Further, based on the magnitude of the ΔE, it can be inferred that MgAl₂O₄ has a strong tendency for a normal spinel structure and MgIn₂O₄ exhibits a relatively stronger tendency for the inverse structure as compared to MgGa₂O₄. These predictions are qualitatively supported by the experimentally observed room temperature structures for these chemistries, as indicated in Table 1. We note, however, that previous work has shown that the difference in inversion tendencies in MgGa₂O₄ and MgIn₂O₄ is due, in large part, to differences in vibrational entropy contributions⁴⁰.

Table 1 Predicted ground states, relative energetics, and lattice constants of the three single spinels along with previously reported experimental and theoretical DFT data.

As a next step, to efficiently examine the relative energetics of various ordered double spinel configurations reached by pairwise-combining any of the two single spinels, we resort to the cluster expansion-based effective Hamiltonian approach described in the Methods section. Toward this aim, three different cluster expansion models were developed for the three possible double spinel chemistries, i.e., MgAl₂O₄-MgGa₂O₄, MgAl₂O₄-MgIn₂O₄, and MgGa₂O₄-MgIn₂O₄. In each case, we enumerate all allowed 53,298 configurations of type [Mg_1−zB_u B′_w]^t[Mg_zB_2−u−vB′_v−w]^oO₄ (where u, v, w, and z are real numbers all ≥0 such that u + w = z, z ≤ 1, 0 ≤ u + v ≤ 2 and v ≥ w; B ≠  B′ and B, B′ ∈ {Al, Ga, In}). A subset of these compounds is then systematically selected to perform DFT computations and to iteratively augment the training dataset used for fitting each cluster expansion model, as discussed in the Methods section. Figures 2a–c present results of the five-fold cross-validated cluster expansion models for the three double spinels, with the insets in each panel comparing the fitted cluster expansion-predicted mixing enthalpies (defined in the Methods section) with those computed using DFT calculations. We note that in each case, the fitted cluster expansion models are able to correctly reproduce the ground state configurations and relative energetics of the three single spinel chemistries shown in Table 1.

Fig. 2: DFT and cluster expansion energetics.

Composition-dependent mixing enthalpies computed using DFT calculations and fitted using the cluster expansion-based effective Hamiltonian approach for a Mg(Al_xIn_1−x)₂O₄, b Mg(Ga_xIn_1−x)₂O₄, and c Mg(Al_xGa_1−x)₂O₄ chemistries. The insets in a–c provide a comparison of the DFT-computed and the cluster expansion-fitted mixing enthalpies. d Atomic structure of the only stable MgAlGaO₄ structure (space group P4₁22) identified in c with a negative enthalpy of mixing. Octahedral and tetrahedral atomic coordination environments are identified via the polyhedra.

Interestingly, the predictions from the cluster expansion models suggest that there is not a single ordered double spinel compound that is thermodynamically stable (i.e., ΔE_mix < 0 eV) for the MgAl₂O₄-MgIn₂O₄ (Fig. 2a) and MgGa₂O₄-MgIn₂O₄ (Fig. 2b) chemistries. In both cases, all the double spinels have a higher energy relative to the respective two single spinel end points, indicating that in these cases any ordered double spinel structure is actually thermodynamically preferred to decompose into two single spinels, regardless of their relative mixing compositions. This is somewhat surprising for MgGa₂O₄-MgIn₂O₄ as conventional wisdom would suggest that two inverse spinels could mix but, in this case, if they do, it would be driven by entropy, not enthalpy. For the MgAl₂O₄-MgGa₂O₄ case shown in Fig. 2c, however, we predict a special configuration, reached by equimolar mixing of the two single spinel compositions, with a significantly negative enthalpy of mixing (a DFT computed ΔE_mix of 0.088 eV/f.u.), which is on the order of the stability of the inverse structure for MgGa₂O₄ relative to the normal spinel structure. Interestingly, this is the only compound that appears on the convex hull for the entire range of compositional mixing in Fig. 2c. To the best of our knowledge, we are not aware of any previous reports of such an ordered double spinel over the entire spinel chemical space. Therefore, as a natural next step, we take a closer look at the cation ordering and electronic structure of the identified ordered double spinel configuration.

Predicted ordered ground state for MgAlGaO₄

The specific cation ordering pattern in the double spinel can easily be understood and rationalized in terms of the ground state cation ordering configurations of the two parent single spinels. As discussed earlier, MgAl₂O₄ and MgGa₂O₄ exhibit normal and inverse spinel ground states with [Mg]^t[Al₂]^oO₄ and [Ga]^t[Mg Ga]^oO₄ configurations, respectively. In the inverse spinel, Mg and Ga exhibit an ordering on the octahedral sublattice which is symmetric with respect to swapping of the two cation chemistries. In other words, if we distinguish the two sets of octahedral sites by o_I (occupied by Mg) and o_II (occupied by Ga), then a symmetry transformation going from [Ga]^t[Mg${\,}^{{o}_{I}}$Ga${\,}^{{o}_{II}}$]O₄ to [Ga]^t[Ga${\,}^{{o}_{I}}$Mg${\,}^{{o}_{II}}$]O₄ renders the single inverse spinel symmetry invariant. We find that the cation ordering of the double spinel MgAlGaO₄ (cf. Fig. 2d) is closely related to the parent inverse spinel with a [Ga]^t[Mg${\,}^{{o}_{I}}$Al${\,}^{{o}_{II}}$]O₄ configuration. The atomic configuration of the identified stable ordered double spinel is exactly identical to the inverse MgGa₂O₄ structure (with the same P4₁22 space group) where all the Ga atoms appearing on the octahedral sublattice have been replaced with Al atoms. Therefore, the double spinel local configurational environments preserve the specific coordinations of their atomic constituents coming from the two end point single spinels (i.e., Al in a normal spinel and Mg and Ga in an inverse spinel configuration).

In Fig. 3, we compare the electronic band structure of the MgAlGaO₄ double spinel with the two parent single spinels. The relaxed geometries for each of the single spinels and the double spinel structure are provided as Supplementary Data 1 accompanying the manuscript. From the band structures in Fig. 3a, b, it can be seen that both the normal MgAl₂O₄ ($Fd\bar{3}m$) and the inverse MgGa₂O₄ (P4₁22) compounds are direct bandgap materials (minimum occurring at the Γ point) with bandgaps of 5.11 and 2.63 eV, respectively (computed using the generalized gradient approximation (GGA) Perdew, Burke, and Ernzerhof (PBE) exchange correlational functional). Further, while the conduction bands have a significant dispersion, the valence bands are relatively flat for both materials. Figure 3c shows that going to the double spinel structure, all the electronic structure features of the parent materials are preserved with an intermediate bandgap of 3.33 eV. To further understand the nature of the specific atomic orbital contributions to the formation of the valence and conduction bands, in Fig. 3d we plot the orbital decomposed density of states for these chemistries. It can be seen that O 2p-states consistently and predominantly contribute to the valence band edge formation in all three cases; Ga 4s-states are largely responsible for the conduction bands in both MgGa₂O₄ and MgAlGaO₄, while the conduction bands in MgAl₂O₄ derive major contributions from Mg 3s- and Al 3p-states. Thus, the addition of Ga adds new electronic states in the double spinel that significantly reduces the bandgap.

Fig. 3: Electronic band structure and density of states.

DFT-computed (GGA-PBE) electronic band structure for a MgAl₂O₄ normal (space group $Fd\bar{3}m$), b MgGa₂O₄ inverse (space group P4₁22), and c MgAlGaO₄ inverse double spinel (space group P4₁22) compounds. d Orbital decomposed density of states for these chemistries shown in a–c.

Finite temperature cation ordering from Monte Carlo simulations

Our DFT and cluster expansion results indicate the existence of a new double spinel structure at 0 K. However, to understand the stability of this structure under realistic conditions, we examine the temperature dependence of the cation structure both theoretically and experimentally. First, to understand the finite temperature cation distribution in the identified 0 K ordered double spinel structure, we present results of our Monte Carlo simulations utilizing the fitted cluster expansion model. While the potential energy or enthalpy E is directly accessible in the Monte Carlo simulations, the configurational entropy S is accessed via thermodynamic integration following E = ∂(βΦ)/∂β and S = −∂Φ/∂T, where β  = 1/k_BT, Φ is the characteristic potential, k_B is Boltzmannʼs constant and T represents the temperature. We note that within this approach vibrational contributions to free energy are not included. However, this neglect is not expected to change our results and the vibrational contributions are expected to be negligible due to the presence of very similar local atomic configurational environments in the single reference spinels and the double spinels⁶⁴.

The Monte Carlo-computed enthalpic and entropic contributions to the free energies as a function of temperature are shown in Fig. 4. While below 300 K the ordered double spinel configuration is preferred, mixing between Mg and Al on the octahedral sublattice starts above this temperature, leading to a random-inverse-like cation ordering. This mixing results in relatively small entropic gains, which are reflected in the free energy profile (cf., inset in Fig. 4). Similarly, mixing of Mg and Al on the octahedral sublattice has little impact on the enthalpy of the material. Further increasing the temperature beyond 500 K sets off intermixing between the disordered octahedral and the tetrahedral sublattices, allowing Ga to appear on the octahedral sublattice. As a result of this inter-sublattice mixing, the induced local strain effects lead to a sharp increase in the mixing enthalpy, while configurational entropy also rises abruptly. As a result, the fully ordered compound is not stable once temperature is introduced, but some level of ordering is still maintained at higher temperatures.

Fig. 4: Finite temperature energetics from Monte Carlo simulations.

Temperature-dependent mixing enthalpy, configurational entropy, and free energy for MgAlGaO₄ computed using Monte Carlo simulations. The inset resolves the free energy in the low-temperature limit and the three regimes that correspond to perfect ordering (yellow), mixing only on octahedral sublattice (green), and mixing on both the sublattices (blue), respectively, with increasing temperature, are highlighted.

The octahedral-tetrahedral site cation mixing in single inverse spinel compounds has been studied in the past both experimentally and theoretically. For instance, for MgAl₂O₄ the tetrahedral-octahedral mixing temperature (${T}_{c}^{t-o}$) as been reported to be 750⁶⁵, 950⁶⁶, and 930 K⁶⁷ in different experimental studies, which agree with the value of 860 K estimated using Monte Carlo simulation³¹. For the inverse spinels MgGa₂O₄ and MgIn₂O₄, the corresponding Monte Carlo-estimated ${T}_{c}^{t-o}$ values are relatively higher (1100 and 1300 K, respectively)³¹. It is interesting to note that the MgAlGaO₄ double spinel exhibits a ${T}_{c}^{t-o}$ which is significantly lower than that reported for the closely related single inverse spinels (MgGa₂O₄ and MgIn₂O₄). This lowering of the critical mixing temperature can be traced back to the availability of unique cation mixing options (such as Ga^t–Al^o) not available in the single spinels.

While disordering between tetrahedral and octahedral sites can be probed using a variety of experimental techniques, including neutron^68,69,70 and X-ray diffraction⁷¹, high-resolution nuclear magnetic resonance⁷², electron spin resonance⁷³, and Raman spectroscopy⁷⁴, the phase transition involving intra-octahedral site mixing is rather subtle and is yet to be studied directly in experiments, to the best of our knowledge. Based on the results of their Monte Carlo simulations, Seko et al.  have reported a critical temperature for octahedral-octahedral mixing (${T}_{c}^{o-o}$) of 270 K and 310 K for MgGa₂O₄ and MgIn₂O₄, respectively, which are close to the value of  ~300 K obtained in the present study³¹. Given that the elemental species participating in the mixing in these single and double spinels are different (Mg with Ga, In, and Al, respectively), these results point towards a very similar qualitative nature of the underlying potential energy surface involved in the octahedral-coordinated site mixing across these chemistries.

Lastly, from Fig. 4, it can be seen that with increasing temperature the computed entropic and enthalpic contributions gradually taper off. However, even at temperatures as high as 3000 K the two curves are not completely flat. This could point to the fact that even at higher temperatures, the cation distributions have not completely randomized in this double spinel and exhibit some amount of short-range ordering. Indeed such a behavior (i.e., a lower than the completely randomized configurational entropy) has also been suggested for MgAl₂O₄ at high temperatures^70,75.

To further understand the temperature evolution of the cation distributions on different sites in the structure, we carried out a detailed analysis of the Monte Carlo-simulated equilibrium structures at different temperatures from 0 K to 3000 K, in discrete steps of 100 K. For these simulations, a 12 × 12 × 12 supercell of the original 28-atom double spinel unit cell containing a total of 48,384 atoms (or 20,736 cations) was used after performing a careful convergence study of the cation distribution profiles with respect to the supercell size employed in the Monte Carlo simulations.

For the coordination-dependent cation distribution analysis, we compute the fractional occupation of sites on the octahedral or the tetrahedral sublattices with each cation type as a function of temperature. Further, to be able to systematically analyze the number of average first, second, and third cation nearest neighbors (i.e., 1NN, 2NN, and 3NN, respectively) at a given temperature, we first compute the radial distribution function (RDF) to determine the cutoff radii for the 1NN, 2NN, and 3NN shells. From the RDF, the 1NN, 2NN and 3NN cut-off radii were determined to be ≤4.5, 4.5–5.5, and 5.5–6.5 Å, respectively. Subsequently, for each cation type (C = Al, Mg, Ga), and coordination site (s = t or o), we calculated the total number of coordination-dependent sites in the nth NN shell (for n = 1, 2, and 3) as ${N}_{{\mathrm{tot}}}^{s}$(nNN), as well as the number of each type of cation for a given coordination and NN shell, ${N}_{C}^{s}$(nNN). Defined this way ${N}_{tot}^{s}$(nNN) = ∑_C${N}_{C}^{s}$(nNN). To quantify the ordering or randomness of the cation distribution at a given temperature, we calculate the ratio ${\Re }_{{C}^{s}}$ = ${N}_{C}^{s}$/$\frac{1}{3}$(${N}_{{\mathrm{tot}}}^{s}$) for the 1NN, 2NN and 3NN shells. Furthermore, we always compute the ratio ${\Re }_{{C}^{s}}$ with respect to a center species $C{^{\prime} }^{s^{\prime} }$ (with $C^{\prime}$ = Al, Mg or Ga and $s^{\prime}$ = o or t) looking radially outward to count specific cations in different NN shells ${\Re }_{{C}^{s}}$($C{^{\prime} }^{s^{\prime} }$) and this procedure is repeated over the entire structure to compute the average statistics. Note that the factor of 3 in ${\Re }_{{C}^{s}}$ comes from the fact that in each shell there are three different types of cations that can occupy either the tetrahedral or octahedral sublattices and for a perfectly random structure the fraction ${\Re }_{{C}^{s}}$ would converge to unity.

The results of our temperature-dependent on-site and NN cation distribution analysis are captured in Fig. 5. The on-site distribution of each atom type at either of the two sublattices is presented in Fig. 5a. More specifically, for each temperature, we plot the total number of atoms for a given cation species found on a specific sublattice divided by the total number of cations of the same elemental species in the entire structure. In a perfectly random distribution, these fractional occupations should converge to $\frac{1}{3}$ and $\frac{2}{3}$ for the tetrahedral and octahedral sites, respectively. The results plotted in Fig. 5a clearly demonstrate that the fractional occupations, computed at a temperature as high as 3000 K, are still far away from those expected in a structure harboring a completely randomized distribution of cations. These deviations in the cation distributions from an ideal random state can potentially be attributed to a short range ordering in this system.

Fig. 5: Temperature-dependent cation distributions in the MgAlGaO₄ double spinel.

a Temperature and coordination environment-dependent fractional occupation of different cation species. b Average fractional occupations of different species computed in the 1NN, 2NN, and 3NN shells of selected cations at specific coordination environments. c Average fractional occupations of different species computed in the 1NN shell of all six configurationally unique cation sites. See text for details.

To quantify the potential short range ordering, in Fig. 5b we graphically represent the ratios ${\Re }_{{C}^{s}}$($C{^{\prime} }^{s^{\prime} }$) computed for 1NN, 2NN and 3NN shells and centered around a selected set of elemental species on a given sublattice type, namely, $C{^{\prime} }^{s^{\prime} }$ = Al^o, Mg^o or Ga^t, over the entire range of temperatures considered here. In particular, for a reference center species (Al^o, Mg^o or Ga^t) the ratio captures an averaged count of specific cation pair types occurring over the specific sublattices in 1NN, 2NN or 3NN shells over the entire structure and converges to unity for a completely randomized structure. Furthermore, for completeness, in Fig. 5c we present all six possible 1NN interactions where the ratios ${\Re }_{{C}^{s}}$ are computed with reference to each of the three cations (Al, Mg, Ga) for the two different coordination sites (o or t). The results from our NN analysis clearly demonstrate strong evidence for short range ordering in this structure. For instance, focusing on the Mg(o)–Mg(o) 1NN, 2NN and 3NN interactions shown in Fig. 5b (and highlighted in the bottom-left inset), it can be seen that while ${\Re }_{M{g}^{o}}$(Mg^o) computed for the 2NN and 3NN shells closely approaches unity at ~3000 K—a temperature well above the melting temperature of MgAl₂O₄—the 1NN ratio is still quite far from the value expected in a structure with completely random cation distribution. This indicates that a complete random distribution of cations is unlikely to occur in this compound even at temperatures near the melting temperature. This entire analysis, based on the Monte Carlo simulations and the on-site and NN local cation distributions, provides strong theoretical evidence for short range ordering in the double spinel at elevated temperatures and is consistent with recent neutron diffraction experiments showing short-range order even after irradiation⁷⁰.

Discussion

Our theoretical calculations predict the existence of a novel inverse double spinel structure for mixed normal MgAl₂O₄ and inverse MgGa₂O₄. This structure is not stable for other mixtures of Mg-bearing spinels, including inverse MgGa₂O₄ and inverse MgIn₂O₄. Our results indicate that the formation of ordered double spinel structures depends sensitively on the nature of the cations in the material. However, our results also suggest that, even in the stable case of MgAlGaO₄, the fully ordered structure is only stable at low temperatures and that, by temperatures of 300 K, mixing will begin on the octahedral sublattice.

To validate the theoretical predictions, we have synthesized the MgAlGaO₄ compound via ball milling, as described in the Methods section. Energy Dispersive X-ray Spectrometry (EDS) (FEI Quanta 400F FEG with a Thermo Noran Systems NSS System 7 EDS detector) was used to confirm the homogeneous dispersion of cations within the material. The EDS characterization results shown in Fig. 6a, b indicate that the material exhibits micron-scale homogeneity of elemental distributions. Furthermore, the synthesized single phase compound reveals a chemical composition of the individual elemental species aligning with that of MgAlGaO₄ (cation percentages of 27% Mg, 36% Al, 37% Ga), within the error bars expected from the technique. The deviation from the anticipated 33% for each element is likely due to the marked topological variation of the dispersed powder, as well as having used a standardless fitting procedure⁷⁶.

Fig. 6: Experimental characterization of synthesized MgAlGaO₄ double spinel.

EDS-based elemental maps a and spectra b showing a uniform distribution of Mg, Al and Ga cations in the as synthesized MgAlGaO₄ double spinel. c XRD patterns for MgAlGaO₄ for anneals at (left) 473 K over 6 weeks and (right) 773 K over 8 weeks in 2 week intervals, starting with material synthesized at 1773 K. (center) The 220 to 331 peak ratios, as extracted from the XRD patterns. The 220 peak is more sensitive to the cation ordering and this ratio gives some measure of the degree to which the cations are ordered in the structure. For reference, the peak ratios as extracted from simulated patterns for the fully ordered double spinel and the random spinel are also shown.

The XRD results highlighting cation ordering evolution as a function of annealing time are presented in Fig. 6c, where XRD patterns for the compound are shown as a function of anneal time for two different temperatures: 473 and 773 K. These results immediately highlight a challenge in validating the theoretical predictions: at the low temperatures at which the ordered MgAlGaO₄ structure is predicted to be stable, the kinetics are simply too sluggish to form that structure. That is, at 473 K, we see no evolution of the structure of the compound over a time scale of 6 weeks.

To overcome kinetic limitations, we annealed the material at 773 K. At this temperature, we do not expect to be able to form the fully ordered double spinel structure, but we do expect that, as shown in Fig. 4, the material to exhibit significant order, and to be more ordered than the starting material synthesized at 1773 K. This is born out by the evolution of the XRD pattern shown in Fig. 6. As time progresses, the peak ratios are much closer to the ordered case than they are at the beginning of the anneal. We do not claim that the fully ordered structure has formed (as might be inferred from Fig. 6) but rather that the material is indeed evolving toward the low-temperature ordered structure. While not serving as a direct confirmation of the theoretical predictions, the XRD analysis of the annealed material does indicate the strong ordering tendency of this compound.

While the XRD results discussed above based on the ratio of the 220–331 peak heights are only semi-quantitative in nature, further analysis showing a more detailed and direct quantitative comparison of the XRD patterns measured for the synthesized compound and fitted to the simulations starting from the predicted computational structure is presented in the Supplementary Figs. 1 and 2 and Supplementary Tables 1 and 2 accompanying the manuscript. Our results from the fitting analysis show that the XRD pattern measured for both the as-synthesized and the 8-weeks annealed samples are in an excellent agreement with the corresponding patterns fitted using the computationally-predicted structure. Further, the best fits are achieved by allowing the occupancies of the three cations on the two sublattices to vary. Fractional site occupancies extracted from the XRD patterns refinements, presented in Table 2, indicate an increase of Ga content from nearly 75–85% on the tetrahedral sublattice, which is again in line with the theoretical predictions.

Table 2 Fractional site occupancies for as synthesized and 8-week annealed at 773 K sample extracted from the XRD patterns refinements.

Together, our theoretical and experimental results indicate that, similarly to the perovskite family of materials, ordered double spinels are likely to form in at least some spinel chemistries. This opens a new avenue for designing materials. As has been shown in previous work, cation ordering in complex oxides (containing multiple cations) has a significant impact on diffusion⁷⁷, magnetic structure⁷⁸, and bandgap⁷⁹. To the extent that even more stable compounds can be identified that can be synthesized under realistic conditions, this provides new possibilities for tailoring functionality by mixing chemistries that do order.

However, as also revealed by our calculations, not all single spinels will mix. While the case we have identified, MgAlGaO₄, might be counter intuitive as it is the result of mixing a normal and an inverse spinel, another case in which two inverse spinels, MgGa₂O₄ and MgIn₂O₄, were mixed did not lead to a stable compound. A third case, mixing MgAl₂O₄ with MgIn₂O₄, was also unstable. Thus, clearly, not all spinels can be mixed to form double spinels but the rules governing which can and cannot mix cannot be simply distilled to mixing of normal and/or inverse structures. Further work is needed to identify the factors that dictate the formation of the double spinel structure. Our results suggest that ionic radius, the biggest difference between Ga and In, is one such factor. A high throughput exploration of a large number of double spinel chemistries is necessary to quantify stability prediction rules for these new compounds and identify the chemical space in which these compounds can form.

Conclusions

Using density functional calculations, we have identified a novel double spinel structure that is thermodynamically favorable for the mixed MgAl₂O₄+MgGa₂O₄ spinel system. This structure is not stable for two other Mg-bearing spinel mixtures. Cluster expansion techniques, combined with canonical Monte Carlo, indicate that, even in the stable case, mixing of cations on the various sublattices will occur at relatively low temperature, suggesting that synthesis of the double spinel would be challenging at best. However, X-ray diffraction of synthesized MgAlGaO₄ highlights the tendency of the compound to order, partially validating the computational results. Together, these results demonstrate the theoretical possibility and experimental reality of double spinel compounds and open new avenues for the design of spinel-structured materials with tailored functionality.

Methods

Cluster expansion and DFT calculations

We are interested in identifying any double spinel configurations that can be formed by pairwise combining the single spinels MgAl₂O₄, MgGa₂O₄, and MgIn₂O₄ over the entire range of composition. Since the three cations forming such a double spinel can occupy either or both the tetrahedral and octahedral lattice sites in different ratios, this leads to a staggering number of possible configurations. For instance, for a 28 atom supercell, the total number of possible unique configurations is nearly 160,000. Since calculating configurational energies for all possible structures using DFT computations is clearly impractical in this case, we resort to a cluster expansion-based effective Hamiltonian approach^{53,54,62,80,81}.

Within this formalism, the energy E of a configuration is expanded as a linear combination of averaged cluster functions ϕ_α as

$$E=\sum_{\alpha }{v}_{\alpha }\cdot \left\langle {\phi }_{\alpha }\right\rangle ,$$

(1)

where the coefficients v_α represent effective cluster interactions (ECIs) for clusters labeled with index α. The cluster function ϕ_α is simply defined as the product of discrete pseudo-spin configuration variables σ_i which form a given cluster α, for respective lattice sites i, as

$${\phi }_{\alpha }=\prod_{i\in \alpha }{\sigma }_{i}.$$

(2)

Since any function of the discrete site-specific configuration variables σ can be expanded in cluster functions, for any given configuration, its energy can be obtained using Eq. (1), provided the ECIs involved in the expansion can be determined. Although in principle, the summation in Eq. (1) runs over an infinite number of clusters that are theoretically necessary to reconstruct the physical system under consideration, in practice, the sum is truncated to include only a finite number of dominant n-body cluster contributions, with the higher term interactions ignored. For this truncated series, next we address the problem of how to select and fit an optimal set of ECIs that can appropriately describe the underlying potential energy surface of the system of interest. To identify and fit the ECIs for a given double spinel structure, we resort to a genetic algorithm with a cross-validated (CV) least square fitting procedure that minimizes the root mean square error (RMSE) between the cluster expansion-predicted and DFT-computed configuration energies. While the use of a genetic algorithm allows us to efficiently explore different combinations of ECIs in a high dimensional space of all possible ECIs, the CV procedure ensures generalizability of the fitted cluster expansion effective Hamiltonian on unseen configurations.

For a given cluster expansion fit, the RMSE ϵ_rms is defined over the differences between the DFT-computed and the cluster expansion-predicted energies for the set of N configurations as follows

$${\epsilon }_{rms}=\sqrt{\frac{1}{N}\mathop{\sum }\nolimits_{i = 1}^{N}{\left({E}_{DFT}^{i}-{E}_{CE}^{i}(\{{v}_{\alpha }\})\right)}^{2}}.$$

(3)

Note that the cluster expansion-predicted energy for a configuration i is a function of a set of ECIs {v_α}. The primary aim of the genetic algorithm-based search is to identify an optimal set of {v_α} that minimizes a CV ϵ_rms (i.e., the error computed on unseen data) as a fitness score. In the present study, a 5-fold CV ϵ_rms was used as the fitness metric. Next, in the genetic algorithm, we define an individual using a d-dimensional binary vector with its different components as either “0” or “1”. Here d indicates the number of maximum possible ECIs allowed in the truncated summation in Eq. (1) and a “0” or “1” appearing in a given component of the binary vector refers to the absence or presence, respectively, of a specific ECI in a cluster expansion fit represented by the particular individual. We start with an initial population of 100 randomly generated individuals and in subsequent generations this population is iteratively evolved by means of standard genetic algorithm operations such as crossovers and mutation. In particular, a size three tournament selection was employed for selecting parents in crossovers while the mutation rate was set to 1.0%. During this evolution process a hall-of-fame catalog of up to the 20 best individuals ever encountered by the algorithm was maintained and, in the end, the best individual in this pool was selected as the final cluster expansion model.

To fit the cluster expansion models, DFT calculations were performed using the Vienna ab initio Simulation Package (VASP)⁸² and employed the Perdew, Burke, and Ernzerhof (PBE)⁸³ generalized gradient approximation (GGA) exchange-correlation functional, which is known to provide a good description of the ground states and relative energetics for the spinel chemistries studied here³⁹. The electronic wave functions were expanded in plane waves up to a cut-off energy of 500 eV. The pseudopotentials based on the projector augmented wave method⁸⁴ explicitly included the following valence electronic configurations for different elemental species in the relevant spinel chemistries—Mg: 2p⁶3s², Al: 3s²3p¹, Ga: 3d¹⁰4s²4p¹, In: 4d¹⁰5s²5p¹ and O: 2s²2p⁴. A Gamma-centered automatically-generated 9 ×  9  ×  9 Monkhorst-Pack k-point mesh⁸⁵ was used for Brillouin-zone integrations for a 14-atom-containing primitive unit cell and appropriately scaled for larger supercells to ensure a nearly similar reciprocal space k-point density. Spin-unpolarized calculations were employed. To obtain a geometry optimized equilibrium structure, atomic positions as well as the supercell lattice parameters were fully relaxed using the conjugate gradient method until all the Hellmann-Feynman forces and the stress component were >0.02 eV/Å and 1.0 × 10⁻² GPa, respectively.

The intrinsic thermodynamic tendency for two single spinels MgB₂O₄ and MgB’₂O₄, with B, B′ ∈ {Al, Ga, In} and B ≠  B′, to form a specific ordered double spinel structure was quantified using the enthalpy of mixing, defined as follows,

$$\Delta {E}_{{\mathrm{mix}}}={E}_{MgB{B}^{\prime}{O}_{4}}-\frac{1}{2}({E}_{Mg{B}_{2}{O}_{4}}+{E}_{Mg{B}_{2}^{\prime}{O}_{4}}),$$

(4)

where ${E}_{MgB{B}^{\prime}{O}_{4}}$ and ${E}_{Mg{B}_{2}{O}_{4}}$, ${E}_{Mg{B}_{2}^{\prime}{O}_{4}}$ represent the DFT-computed enthalpies per formula unit (f.u.) for a given double spinel configuration and for the two single spinel ground state configurations, respectively. For each of the three possible double spinel chemistries, we considered a set of 53,298 configurations possible within a 28 atom supercell. For DFT computations, starting with a randomly selected set of 20 configurations for each double spinel chemistry, additional configurations were added following an iterative process where the next set of configurations were selected based on the predictions of a cluster expansion-model fitted on the present set. At every step, we included and retained all the configurations predicted to have negative enthalpies of mixing relative to the respective single spinel ground state configurations.

Finite temperature thermodynamic properties and temperature-dependent cation distributions were evaluated using canonical Monte Carlo simulations by the Metropolis algorithm^86,87. Supercells for the Monte Carlo simulations were constructed by a 12 × 12 × 12 expansion of the primitive unit cell and contained 24,192 atoms (10,368 cations). The simulations were performed on 1000 trial steps per cation for calculating thermodynamic averages of energy and cluster functions after equilibration over 1000 Monte Carlo steps per cation. The temperature intervals of the Monte Carlo simulations were set to 100 K from 0 K to 3000 K. The cluster expansion-model development and Monte Carlo simulations were performed using the CASM software package [https://github.com/prisms-center/CASMcode]^88,89.

Experimental methods

Commercial powders of MgO (99.999% purity), Al₂O₃ (99.999% purity) and Ga₂O₃ (99.999% purity) were calcined at 1000 °C for 24 h. While the powders were at 500 °C, they were weighed to yield the chemistry of 50:50 atomic percentage of MgAl₂O₄:MgGa₂O₄. The three end-member powders were mixed into a slurry with spectroscopic grade ethanol and placed in zirconia vials with two zirconia balls. The vials were put into a Spex 8000M ball-mill and milled for a total of 8 h in 30 min intervals with 15 min cooling cycles to prevent the vials from heating. Following the ball-milling sequence, the powders were removed from the vials and placed into a drying oven at 200 °C overnight to evaporate the ethanol dispersant. The dried powder was ground in a mortar and pestle into a fine powder. The powder was then pressed in a 13 mm diameter die at 500 MPa. The pellets were ejected and placed into an ambient environment furnace and ramped at a rate of 10 °C/min to 1500 °C and held for 36 h, then allowed to furnace cool. The resulting material, after the entire preparation procedure discussed above was complete, was measured using XRD to determine the resulting phase. The XRD measurements were done on a Bruker Nano D8 Advance with a Cu K-α source using locked couple θ-2θ scans. The XRD measurement revealed a near phase pure spinel-type phase belonging to space group $Fd\bar{3}m$. A minor, almost undetectable, secondary phase of zirconia was also observed which was most likely from the zirconia vial/balls of similar hardness to the MgO and Al₂O₃ powders.

XRD Rietveld refinements were performed on XRD patterns obtained from the two specimens, namely the as-synthesized and the sample annealed at 500 °C for 8 weeks. The structural refinements were performed using the Rietveld refinement code TOPAS⁹⁰. Fitting was performed exclusively using the fundamental parameters convolution based model. Specifically for the calculated predicted phase of interest, in addition to the lattice parameter and scaling factor the oxygen positions and the cations occupancy for each specific cation site were also refined from the measured patterns. The fractional occupancies for the three cation sites were mathematically constrained to have their summation equal unity to maintain the original stoichiometry of the starting fabricated material.

Energy Dispersive X-ray Spectrometry (EDS) (FEI Quanta 400F FEG with a Thermo Noran Systems NSS System 7 EDS detector) was used to confirm the homogeneous dispersion of cations within the material. MgAlGaO₄ powder was dispersed in acetone and applied to a carbon stub using a dropper. EDS spectra were obtained with an incident beam voltage of 5 KV, as this was 2–3 times higher than the highest energy Ga L lines probed. A standardless quantification using the ZAF correction method was employed for measuring cation percentages.