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

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 MgAlGaO4, formed by an equimolar mixing of MgAl2O4 normal and MgGa2O4 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. Materials with a spinel structure are used in various applications, including in the nuclear industry and as dielectrics. Here, first principle calculations and Monte Carlo simulations predict that an ordered double spinel structure is stable, supported by preliminary experimental data.

S pinel compounds-typically represented as AB 2 O 4 , 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 6 and have traditionally been utilized as precious stones, such as red or blue spinels containing either Cr +3 or Fe +2 / Zn +2 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 9 . 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 14 , and catalytic 15,16 properties useful for a diverse range of applications, including data storage 17 , high frequency electronic devices 18 , dielectrics 19 , transparent conducting oxides [20][21][22][23] , lasers 24,25 , sensing 2 , energy storage 26,27 , superconductivity 28 , and biotechnology 29 .
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, 1 8 of the tetrahedral and 1 2 of the octahedral voids are occupied by the A and B metal atoms, respectively in a face-centeredcubic close-packed oxygen sublattice. Thus, a normal spinel can be represented by [ 4 , 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] 4 30 . 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 [ 4 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 2 3 for a completely random distribution of the  30,31 .
The problem of (a) identifying factors that govern the intrinsic tendency of a given binary (i.e., A 3 O 4 , with A = B and multiple possible charge states of A) or ternary (i.e., AB 2 O 4 ) 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 hightemperature corrosion of ferritic-martensitic steels has been attributed to the formation of an inverse magnetite spinel (Fe 3 O 4 ) along with an iron-chromium normal spinel of composition FeCr 2 O 4 , which tend to phase separate 44 .
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−x Ga x O 4 quaternary system (with x = 0.38, 0.76, 0.96, 1.52) have been investigated by Ito et al. 42 ; however, no low temperature ordered double spinel structures formed by combining two single spinels (in this case, MgAlO 4 and MgGaO 4 ) 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 50 , 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 O 4 ) 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 2 O 4 (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 3m and tetragonal P4 1 62,63 . Further, based on the magnitude of the ΔE, it can be inferred that MgAl 2 O 4 has a strong tendency for a normal spinel structure and MgIn 2 O 4 exhibits a relatively stronger tendency for the inverse structure as compared to MgGa 2 O 4 . 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 2 O 4 and MgIn 2 O 4 is due, in large part, to differences in vibrational entropy contributions 40 .
As a next step, to efficiently examine the relative energetics of various ordered double spinel configurations reached by pairwisecombining 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 , 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.
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 2 O 4 -MgIn 2 O 4 ( Fig. 2a) and MgGa 2 O 4 -MgIn 2 O 4 ( 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 2 O 4 -MgIn 2 O 4 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 2 O 4 -MgGa 2 O 4 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 2 O 4 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 4 . The specific cation ordering pattern in the double spinel can easily be understood and rationalized in terms of the ground state cation Table 1 Predicted ground states, relative energetics, and lattice constants of the three single spinels along with previously reported experimental and theoretical DFT data.
This study (DFT-PBE) Past theory (DFT-LDA) Experiments The experimental structures and lattice constants are from ref. 61 , while the theoretical results, computed with the local density approximation (LDA) to the exchange correlations functional within DFT, are from ref. 31 . A positive (negative) ΔE represents that the normal (inverse) spinel structure is more stable.
ordering configurations of the two parent single spinels. As discussed earlier, 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 4 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 2 O 4 (Fd 3m) and the inverse MgGa 2 O 4 (P4 1 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 2pstates 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 2 O 4 and MgAlGaO 4 , while the conduction bands in MgAl 2 O 4 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.
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 B T, Φ 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 64 .
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.
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 2 O 4 the tetrahedraloctahedral mixing temperature (T tÀo c ) as been reported to be 750 65 31 . It is interesting to note that the MgAlGaO 4 double spinel exhibits a T tÀo c which is significantly lower than that reported for the closely related single inverse spinels (MgGa 2 O 4 and MgIn 2 O 4 ). 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-70 and X-ray diffraction 71 , high-resolution nuclear magnetic resonance 72 , electron spin resonance 73 , and Raman spectroscopy 74 , 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 oÀo c ) of 270 K and 310 K for MgGa 2 O 4 and MgIn 2 O 4 , respectively, which are close to the value of~300 K obtained in the present study 31 . 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 2 O 4 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 s tot (nNN), as well as the number of each type of cation for a given coordination and NN shell, N s C (nNN). Defined this way N s tot (nNN) = ∑ C N s C (nNN). To quantify the ordering or randomness of the cation distribution at a given temperature, we calculate the ratio < C s = N s C / 1 3 (N s tot ) for the 1NN, 2NN and 3NN shells. Furthermore, we always compute the ratio < C s with respect to a center species C 0 s 0 (with C 0 = Al, Mg or Ga and s 0 = o or t) looking radially outward to count specific cations in different NN shells < C s (C 0 s 0 ) and this procedure is repeated over the entire structure to compute the average statistics. Note that the factor of 3 in < 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 < 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 1 3 and 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.
To quantify the potential short range ordering, in Fig. 5b we graphically represent the ratios < C s (C 0 s 0 ) computed for 1NN, 2NN and 3NN shells and centered around a selected set of elemental species on a given sublattice type, namely, C 0 s 0 = 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 < 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 < Mg 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 2 O 4 -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 70 .

Discussion
Our theoretical calculations predict the existence of a novel inverse double spinel structure for mixed normal MgAl 2 O 4 and inverse MgGa 2 O 4 . This structure is not stable for other mixtures of Mg-bearing spinels, including inverse MgGa 2 O 4 and inverse MgIn 2 O 4 . 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 4 , 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 4 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 4 (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 76 .
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 4 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 8weeks annealed samples are in an excellent agreement with the corresponding patterns fitted using the computationallypredicted 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.
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 77 , magnetic structure 78 , and bandgap 79 . 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 4 , 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 2 O 4 and MgIn 2 O 4 , were mixed did not lead to a stable compound. A third case, mixing MgAl 2 O 4 with MgIn 2 O 4 , 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 2 O 4 +MgGa 2 O 4 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 4 Table 2 Fractional site occupancies for as synthesized and 8-week annealed at 773 K sample extracted from the XRD patterns refinements.

Site type
Site multiplicity and symmetry Atom type As synthesized occupancy Occupancy after annealing Thermal parameter (B eq ) Al 4 b  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.
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 2 O 4 , MgGa 2 O 4 , and MgIn 2 O 4 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 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 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 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 ddimensional 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) 82 and employed the Perdew, Burke, and Ernzerhof (PBE) 83 generalized gradient approximation (GGA) exchangecorrelation functional, which is known to provide a good description of the ground states and relative energetics for the spinel chemistries studied here 39 . 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 84 explicitly included the following valence electronic configurations for different elemental species in the relevant spinel chemistries-Mg: 2p 6 3s 2 , Al: 3s 2 3p 1 , Ga: 3d 10 4s 2 4p 1 , In: 4d 10 5s 2 5p 1 and O: 2s 2 2p 4 . A Gamma-centered automaticallygenerated 9 × 9 × 9 Monkhorst-Pack k-point mesh 85 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 −2 GPa, respectively.
The intrinsic thermodynamic tendency for two single spinels MgB 2 O 4 and MgB' 2 O 4 , 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, where E MgBB 0 O 4 and E MgB 2 O 4 , E MgB 0 2 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 2 O 3 (99.999% purity) and Ga 2 O 3 (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 2 O 4 :MgGa 2 O 4 . 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 3m. 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 2 O 3 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 90 . 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 4 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.

Data availability
The data that support the findings of this study are available from the corresponding authors upon reasonable request.

Code availability
The VASP code used in this study is a commercial electronic structure modeling software, available from https://www.vasp.at. The CASM code used for Monte Carlo simulations is openly available at https://github.com/prisms-center/CASMcode. Received: 13 May 2020; Accepted: 12 October 2020;