Modeling the high-temperature phase coexistence region of mixed transition metal oxides from ab initio calculations

Suzanne K. Wallace,¹,²,* Ambroise van Roekeghem,¹,² Anton S. Bochkarev,³ Javier Carrasco,⁴
Alexander Shapeev,⁵ and Natalio Mingo¹,²,†

¹CEA, LITEN, 17 Rue des Martyrs, 38054 Grenoble, France

²Université Grenoble Alpes, 621 Avenue Centrale, 38400 Saint-Martin-d'Hères, France

³Interdisciplinary Centre For Advanced Materials Simulation, Universitätsstraße 150, Ruhr-Universität Bochum, 44801 Bochum, Germany

⁴Centre for Cooperative Research on Alternative Energies (CIC energiGUNE), Basque Research and Technology Alliance (BRTA),
Alava Technology Park, Albert Einstein 48, 01510 Vitoria-Gasteiz, Spain

⁵Skolkovo Institute of Science and Technology, Skolkovo Innovation Center, Nobel St. 3, Moscow 143026, Russia

<figure><img src="image_1.png" /></figure>

(Received 28 August 2020; revised 27 October 2020; accepted 1 December 2020; published 12 February 2021)

Accurate knowledge of phase coexistence regions, i.e., solubility gaps (SGs), is key to the development of mixed transition metal oxides for various applications, such as thermochemical energy storage, or catalysis. However, predicting a SG from first principles in these materials is particularly challenging due to the complex interplay between several sources of entropy, the large configuration space, and the computational expense of *ab initio* calculations. We present an approach that yields an accurate prediction of the experimental Hausmannite-spinel SG in the case of (CoₓMn₁₋ₓ)₃O₄. The method uses machine learning to extend an *ab initio* dataset of hundreds of structures, and it includes many different entropic contributions to the free energy. We demonstrate and quantify the crucial roles of phonon and paramagnetic entropy, and the importance of sampling higher-energy configurations, and correcting for finite-size limitations in the *ab initio* supercell configurations.

DOI: 10.1103/PhysRevResearch.3.013139

I. INTRODUCTION

Reversible redox reactions of transition metal (TM) oxides are a promising route to high-density thermochemical solar energy storage at the high operating temperatures of next-generation concentrated solar power plants [1]. Another application envisaged for TM oxides is in catalyzing the oxygen reduction reaction in anion exchange membrane fuel cells [2]. To make a sustainable use of these technologies, the material’s storage capacity and cyclability needs to be balanced against its embodied energy cost [3]. Determining the optimal material for this application in terms of cost, toxicity, cyclability, and energy storage density is an ongoing area of research [1,4–9]. Of the candidate metal oxides that have already been considered for thermochemical energy storage, the Co₃O₄/CoO redox pair stands out for superior energy storage density, reaction kinetics, and cycling stability while the Mn₂O₃/Mn₃O₄ pair is not hindered by the cost and toxicity issues of the Co compound, but suffers from slow oxidation kinetics and poor reversibility. CuO/Cu₂O and Fe₂O₃/Fe₃O₄ have also been considered but are known to suffer from sin-

tering issues, which reduces the cycling stability of the system [1].

Solid solutions of TM oxides have been suggested as a means to mitigate the shortcomings of each binary compound [1,4]. Among the many factors to be taken into account, it has been noted that the reaction reversibility of all these systems is hindered if the cubic spinel structure (S) in the reaction coexists with a tetragonally distorted spinel phase, Hausmannite (H) [1,10]. Thus, to theoretically assess the potential of new mixed TM oxide systems it is important to be able to predict the coexistence region of these two phases in the phase diagram.

Despite the promising qualities of (CoₓMn₁₋ₓ)₃O₄, no first principles calculation of high-temperature phase diagrams for the Co-Mn-O system appears to have been published [1]. Reference [4] has provided zero kelvin formation energy *ab initio* data for hundreds of supercell configurations of (CoₓMn₁₋ₓ)₃O₄. In principle, this data can be used to obtain the free energy and the corresponding phase diagram for this system, including the phase coexistence region, also referred to as the solubility gap (SG). Therefore, the main question we want to answer is whether the SG at high temperatures can be reliably predicted from first principles. And if so, what are the relevant physical mechanisms contributing to it? Often, calculations of phase diagrams employ the Calphad method [11] based on a database of thermochemical information for the relevant phases and minimization of the Gibbs free energy. Other approaches use *ab initio* calculations combined with cluster expansion and Monte Carlo methods for sampling configuration space [12–14], and another work has developed

Published by the American Physical Society under the terms of the
Creative Commons Attribution 4.0 International license. Further
distribution of this work must maintain attribution to the author(s)
and the published article's title, journal citation, and DOI.

<footnote>* suzanne.wallace@cea.fr</footnote>

<footnote>†natalio.mingo@cea.fr</footnote>
a method to augment the standard cluster expansion approach with machine learning [15]. Here we have used the *ab initio* data of Ref. [4] to train a machine learning model to fill in the configuration space of the mixed phase and then consider various other entropic contributions to the free energy.

In what follows, we show that the temperature and composition dependence of the SG can be predicted, but this requires including vibrational entropy and sampling configurations with energies several eV per supercell higher than the actual ground state at each composition (which were absent from Ref. [4] data). We also show that, when employing finite supercells for *ab initio* calculations, it is essential to correct the values of the configurational entropies near the composition edges. Furthermore, we quantify the relative importance of the different contributions to the free energy in the case of $(Co_xMn_{1-x})_3O_4$, which one would need to consider when screening materials for sustainable thermochemical energy storage.

## II. $(Co_xMn_{1-x})_3O_4$

The ground-state crystal structure of the pure end members of the $(Co_xMn_{1-x})_3O_4$ solid solution are cubic spinel in the case of $Co_3O_4$ and a tetragonally distorted spinel for $Mn_3O_4$, where in each case transition metal (TM) ions occupy either tetrahedral (td) or octahedral (oh) lattice sites. Visualizations of the crystal structures are given in the Appendix (Fig. 6). There is a tetragonal (H) to cubic (S) phase transition for $Mn_3O_4$ at a critical temperature $T_C = 1443$ K [16]. Symmetry arguments indicate that the H-S transition of $Mn_3O_4$ is a first-order transition [17]. Many of the TM ions in the end members of the solid solution have nonzero spin magnetic moments. At very low temperatures (below 40 K and 42 K, respectively), $Co_3O_4$ exhibits antiferromagnetic ordering of the atomic spins while $Mn_3O_4$ exhibits ferrimagnetic ordering with a small net magnetization [18]. All $(Co_xMn_{1-x})_3O_4$ compositions have also been found to be ferrimagnetic with measured magnetizations between 0 and $1.5\mu_B$ per formula unit and a paramagnetic transition at very low temperatures, comparable to those of the pure end members [19].

Density-functional theory (DFT) data for $(Co_xMn_{1-x})_3O_4$ from Ref. [4] contain the ground state total energies of relaxed 56-atom (8 formula unit) supercells. This $T = 0$ K data shows a composition-dependent transition from a cubic to tetragonal ground-state structure with increasing Mn content, in good agreement with experimental measurements [18,20]. The plot of formation energy against composition from Ref. [4] shows that all formation energies across the composition range $0 < x < 1$ are greater than zero. Therefore, within this theoretical framework, any concentration would phase separate into a mixture of the two pure oxides. As these solid solutions have been synthesized in various studies [1,4], it is clear that a zero-kelvin *ab initio* model is not sufficient to reproduce the key characteristics of this system. Therefore, we now go on to describe how this data can be used to build a model to predict the SG for $(Co_xMn_{1-x})_3O_4$.

## III. MODEL FOR SOLUBILITY GAP PREDICTION

The SG (or phase coexistence region) of a solid solution can be calculated from its Helmholtz free-energy curve,

$F(T, x)$, as a function of the concentration $x$ [21,22]. This process is described in Appendix B. The method to compute free energies is well described in many references [12–14,23,24]. In brief, $F(T, x)$ is obtained from

$$F = -k_B T \ln Q, \quad (1)$$

where $Q$ is the configurational partition function given by

$$Q(T, x) = \sum_{n} e^{-E_{n}(x)/k_{B}T}, \quad (2)$$

where $E_n(x)$ are the total energies of the different microstates of the system, indexed by $n$. In the following sections, we systematically include more possible microstates in our calculation of $F$ from various entropy sources and quantify their impact on the accuracy of the prediction for the H-S SG of $(Co_xMn_{1-x})_3O_4$.

### A. Configurational entropy

When considering configurational entropy (from the various possible atomic arrangements in the system), the ideal-solution approximation is sometimes applied to simplify the computation of $Q$ [22,25]. This approximation disregards energy differences between inequivalent configurations. However, DFT data for $(Co_xMn_{1-x})_3O_4$ from Ref. [4] displays a clear energy difference at most compositions between the two different substitution schemes considered in the study, with the energy difference exceeding 0.4 eV per formula unit at certain compositions. With large energy differences between different configurations in comparison to $k_B T$, the ideal-solution model is not unsuitable in this case (see a comparison of the configurational entropy from our final model to that predicted by the ideal-solution approximation in the Appendix, Fig. 10). To obtain $E_n(x)$ for all possible atomic configurations of the supercell in the computation of the partition function in Eq. (2), we have used the data from Ref. [4] to build a machine-learning model based on moment tensor potentials [26,27] and predict the total energies for the remaining supercell configurations. Details of this procedure will be presented elsewhere [28].

It has been suggested that metastable configurations slightly higher in energy than the ground state may be responsible for the stability of $(Co_xMn_{1-x})_3O_4$ [4]. As $T$ rises, and more of these higher-energy states become accessible, the increasing configurational entropy will lower the free-energy curve in the mid-range of $x$, allowing concentrations other than 0 and 1 to be stable, and progressively closing the SG. $F(T, x)$ calculated considering only additional microstates from different atomic configurations is shown by dashed lines in Fig. 1(a) and the corresponding SG is shown in purple in Fig. 1(b). The figure also shows the experimental phase equilibria diagram for this system taken from Ref. [20] in gray for H, S, and rocksalt (R) $Co_xMn_{1-x}O$ phases. In this work we consider only the H and S phases (and their coexistence region) because our initial dataset did not contain any R structures. It can be seen from Fig. 1(b) that the current model does not reproduce the H + S coexistence region because other important entropy contributions are currently neglected.
<figure><img src="image_2.png" /><figcaption>FIG. 1. (a) Free energy of (CoₓMn₁₋ₓ)₃O₄ vs concentration x of Co at various temperatures computed by considering only configurational entropy contributions (dashed lines) or configurational and paramagnetic entropy (solid lines). (b) Corresponding H (Hausmanite) + S (cubic spinel) solubility gap calculated by including only configurational (purple) or configurational and paramagnetic (blue) entropy with experimental data from Ref. [20] shown in gray for the H and S phases and also for the S + R phase coexistence region, where R is the CoₓMn₁₋ₓO rocksalt phase.</figcaption></figure>

## B. Magnetic entropy

Each atomic configuration considered in the previous section also has associated with it a certain number of microstates due to magnetic degrees of freedom. With the low paramagnetic transition temperatures of (CoₓMn₁₋ₓ)₃O₄, there will be orientational disorder of localized magnetic moments. At temperatures far above the typical magnetic ordering temperatures, such as the 1200–1443 K temperature range of interest in this study, we assume that the system is in a completely paramagnetic phase where all of the atomic spins can assume any possible arrangement. In this scenario, every atom has associated with it (2S + 1) degrees of freedom, where S here denotes the spin quantum number, related to the atomic magnetic moment by $μ = Mμ_B = 2Sμ_B$ [12,29,30]. It is possible that certain frustration effects (especially in defective or inhomogeneous systems) limit the magnetic degrees of freedom even above the critical temperature [31]. However, we assume here that the temperatures are sufficiently far above the typical ordering temperatures that this is unlikely to result in much (if any) overestimation of the magnetic entropy contribution.

Taking the theoretical atomic magnetic moments of the end members of the mixed phase [18], this gives values for M of

$\pm 5μ_B$ for td-coordinated Mn, $\pm 4μ_B$ for oh-coordinated Mn, $\pm 3μ_B$ for td-coordinated Co, and $0μ_B$ for oh-coordinated Co, which to one significant figure are the same as those obtained for the mixed phases in Ref. [4]. Therefore, the total number of magnetic degrees of freedom for a configuration of the mixed phase depends on the number of Co ions, and also strongly depends on the number of them that occupy oh sites, because this results in replacing a high-spin state atom with a zero-spin state atom. Many works have postulated a possible thermally induced partial low-to-high spin-state transition of oh-coordinated Co in Co₃O₄ [32–35]. However, here we have neglected this possibility and consider the spin state of this species to always be zero. First, because we are interested in the lower Co-concentration side of the phase diagram for the H + S SG and second because it is unclear if the same effect should be expected in the mixed phases.

The magnetic entropy $S_{mag}$ in the paramagnetic state can be obtained for the supercells from all of the atoms with nonzero M by the mean-field term [12,29,36]

$$S_{\text{mag}} = k_B \ln \left( \prod_i (M_i + 1)^{i_{\text{tot}}} \right), \quad (3)$$

where i is each TM species [distinguishing when they are on an oh or td site as this impacts the value of $M_i$ in the case of (CoₓMn₁₋ₓ)₃O₄] and $i_{tot}$ is the total number of a particular species on a particular type of crystallographic site. To account for all possible magnetic states associated with each configuration, we make the assumption that the spin arrangements of the same configuration are close enough in energy that we can multiply each configuration by its number of magnetic states to determine the total number of states for this configuration at this particular energy.

Including magnetic entropy in the free-energy calculation modifies the free-energy curve very noticeably, as can be seen by comparing the dashed and solid curves in Fig. 1(a). However, it does not make the SG resemble the experimental one any more closely, as can be seen from the blue shaded region in Fig. 1(b). There are clearly important components of the model that are still missing, with other possible entropic contributions including electronic and vibrational entropy [12].

## C. Electronic entropy

For an insulator or sufficiently-wide-gap semiconductor, we could expect the contributions to entropy from thermally induced electronic excitations to be negligible. There are measurements in the literature for two different compositions of (CoₓMn₁₋ₓ)₃O₄ which indicate that the compounds are semiconductors with band gaps of around 1.5 e V [37], which is a large reduction relative to Co₃O₄ and Mn₃O₄, but still sufficiently large that we could justify the neglect of electronic excitations. However, first principles calculations of the electronic density of states (EDOS) across the full composition range of (CoₓMn₁₋ₓ)₃O₄ in Ref. [18] have shown the appearance of intermediate states in the band gap and much-reduced fundamental band gaps, which were also very variable depending on the level of theory used.

To determine the impact of electronic entropy on the F of (CoₓMn₁₋ₓ)₃O₄, we have calculated the EDOS of the lowest
and second lowest energy structure at each composition in the dataset from Ref. [4] and checked for cases with significant DOS close to the Fermi level. We use the calculated EDOS to obtain an over-estimated upper limit for the number of accessible electronic states for all compositions of (CoₓMn₁₋ₓ)₃O₄ which allows us to determine a maximum range for the electronic entropy we could expect across the full composition range. Further details are included in the Appendix, but, in brief, we obtain a maximum range of accessible electronic states of 0 to 207 per supercell. To put this range into context, the minimum and maximum number of magnetic states per supercell across all possible configurations are 5 764 801 and 3.972 × 10²³, respectively. The number of accessible electronic states (and any variation between different configurations) will clearly be a much smaller contribution compared with magnetic states in the case of (CoₓMn₁₋ₓ)₃O₄. However, for other mixed compounds without magnetic atomic moments or with a stronger metallic character, electronic entropy might be a more dominant contribution that should not be neglected.

## D. Vibrational entropy

In the absence of contributions from electronic entropy for even the mixed phases with the narrowest band gaps, it is necessary to look to contributions from vibrational entropy. Phonon entropy effects are sometimes disregarded in phase stability calculations because they are difficult to calculate from first principles. However, it has been demonstrated that their effect can be large [38–41]. The configurations sampled by Zaki *et al.* in Ref. [4] contain both cubic and tetragonal structures, with formation energies in the range 0–2 eV per formula unit. In principle, one would not think it necessary to include any configurations lying higher in energy, since their occupation probability may seem negligible at the temperatures of interest. However, consideration of vibrational degrees of freedom, or phonon entropy, reveal the need to sample over cubic systems with considerably higher formation energies. This is due to the larger vibrational entropy of the cubic phases, as compared with the tetragonal ones, in the Mn-rich side of the phase diagram. Because of this, the original cubic data in Ref. [4] are not sufficient to describe the thermodynamically stable phases of the system at high temperature.

Insights for the vibrational entropy of (CoₓMn₁₋ₓ)₃O₄ can be gained from the behavior of the pure Mn₃O₄ compound. It is sometimes the case that the spectral differences across compositions for a given crystallographic phase are smaller than the differences between two crystallographic phases of the same composition [42]. At the cubic-tetragonal phase transition temperature T_C of pure Mn₃O₄, the free energy of the cubic phase should equal that of the tetragonal phase and cross below it upon further temperature increase. Since there is no species disorder in these pure phases and electronic entropy can be considered negligible for electrical insulators, this crossover can only be driven by differences in vibrational and/or magnetic entropy between the two phases. To see how this happens, let us consider the different vibrational spectra of the cubic and tetragonal phases of Mn₃O₄.

As we show below, the tetragonal-cubic transition of Mn₃O₄ is most likely caused by the different vibrational entropy of the two phases, although we cannot discard a possible additional influence from the different magnetic entropy of the two phases. To quantify the vibrational free energy of the two phases, we need the phonon density of states for each of them. This is done via the calculation of phonon frequencies. A standard approach is to compute the interatomic force constants (IFCs) of a large-enough supercell using finite displacements, generate the dynamical matrix, and solve the eigenvalue problem for a grid of phonon wave vectors in the Brillouin zone. The density of states is then obtained from the collection of frequencies in this grid by using a variety of approaches such as Gaussian smearing or the tetrahedron method. We used phonopy [43] to perform this task for the tetragonal phase, using the IFCs as calculated for the magnetic system. It is not strictly correct to use the 0 K ordered magnetic solution to represent a paramagnetic system at ≈1400 K. However, proper inclusion of magnetism at high temperatures in DFT is a highly specialized task and is beyond the scope of this article [36]. The tetragonal phase vibrational DOS is given in Sec. 1 of the Supplemental Material (SM) [44]. Its calculated vibrational free energy at 1440 K is −2046.3 k J − mol, and the vibrational entropy is 2838.8 J K⁻¹ mol⁻¹.

When trying to perform the same finite-displacement calculation for the cubic phase, which is unstable at 0 K, the resulting spectrum yields unphysical imaginary frequencies. To calculate the high-temperature vibrational DOS of a phase that is not stable at 0 K, we first need to obtain its effective high-temperature force constants using a self-consistent stochastic approach via the QSCAILD code (for more details see Refs. [45–47]). In brief, the method consists of calculating the energies and forces of randomly displaced atomic configurations u⁽ⁿ⁾ generated according to their thermal probability distribution at that temperature. Then the effective IFCs of the system are chosen so as to minimize the sum of the squared differences between the DFT forces and those predicted from the adjusted IFCs Φ,

$$ \sum_{n \in \text{figs}} \left| F_{i}^{(n), \text{DFT}} - \sum_{j} \Phi_{i,j} u_{j}^{(n)} \right|^{2}, \quad (4) $$

where i and j run over the degrees of freedom in the supercell. A new set of displaced configurations is then generated by using the IFC-predicted frequencies to obtain the new thermal probability distribution. The forces of the new configurations are again computed by DFT and the approach is repeated in a loop with a mixing algorithm, until the effective IFCs no longer change.

Using this approach we are able to obtain a converged vibrational DOS for the cubic phase at 1500 K. This calculation is orders of magnitude more expensive than those using finite displacements. Therefore, we take these self-consistent high-temperature IFCs to be approximately valid for high temperatures around the transition temperature. The corresponding vibrational free energy at 1400 K is −2415.0 k J − mol, and the vibrational entropy is 3089.7 J K⁻¹ mol⁻¹. We also tried to perform a self-consistent high-temperature calculation for the tetragonal phase at 1500 K, but this systematically yielded imaginary frequencies
<figure><img src="image_3.png" /><figcaption>FIG. 2. Vibrational free energy, $F_{vib}$, versus temperature for the tetragonal (tet) and cubic (cub) phases (solid lines) and their entropic contributions, $TS_{vib}$ (dashed lines), as a function of temperature.</figcaption></figure>

and did not reach a converged solution. The tetragonal phase is unstable at high temperature, so there is in principle no guarantee that the self-consistent approach should reach convergence, because the system may be trying to transition to the more stable cubic phase but it is not able to reach it due to the unit cell being forced to remain tetragonal.

Figure 2 shows the vibrational free energies of the cubic, $F_{vib}^{cub}$, and tetragonal, $F_{vib}^{tet}$, phases of Mn$_3$O$_4$, with that of the cubic being always lower than that of the tetragonal and this difference increasing with temperature. From the calculated values of the vibrational free energies at the experimental transition temperature, we obtain a vibrational energy difference between the two phases of $F_{vib}^{cub}(T_C) - F_{vib}^{tet}(T_C) \simeq -368.6$ kJ/mol $\simeq -3.82$ eV/supercell. Our DFT-calculated difference in total formation energy between the cubic and tetragonal phases is 3.80 eV/supercell. So this amount of vibrational free energy is able to bring down the total free energy of the cubic phase to coincide with the tetragonal one.

Such good agreement would indicate that any effect coming from differences in the magnetic entropy of the two systems must be small. However, this is difficult to ascertain since the typical errors of phase-transition temperatures predicted by self-consistent ab initio lattice dynamics calculations can be as large as 100 K [48,49]. Looking at the calculated temperature-dependent free-energy difference, a 100 K miscalculation of the transition temperature would be associated with an additional 0.14 eV/supercell that would have to be supplemented by differences in magnetic free energy. This is still less than 4% of the total free-energy difference leading the phase transition, which is therefore overwhelmingly due to phonon entropy.

From Fig. 3, the temperature dependence of the free-energy difference that we include in our model is approximately

$$ \Delta F(T) \simeq \Delta E \left( 1 - \frac{T}{T_C} \right), \qquad (5) $$

where $\Delta E$ is the energy difference between the cubic and tetragonal phases of Mn$_3$O$_4$. The same relationship can also be derived from an analytical approximation which is included in Sec. 2 of the SM [44].

A main obstacle to obtaining the SG from the dataset of Ref. [4] is its lack of data for the cubic configurations with

<figure><img src="image_4.png" /><figcaption>FIG. 3. Vibrational free-energy difference ($\Delta F_{vib}$) between the tetragonal and cubic phases, compared with its entropic contribution ($T \Delta S_{vib}$), as a function of temperature. The 3.8 eV/supercell marked by the horizontal dashed line corresponds to the DFT-calculated formation energy difference between the two phases at 0 K. This energy is almost perfectly compensated by the vibrational free-energy difference at the phase-transition temperature (1443 K, vertical dashed line).</figcaption></figure>

$x < 1/3$. These are missing because standard total-energy minimization algorithms in DFT codes aim at finding a global minimum, so high-energy structures that are metastable at 0 K are typically not found. But without these high-energy cubic structures it is not possible to predict the SG. As an approximation, we can assume that every calculated tetragonal configuration has a corresponding cubic one, with a higher formation energy, but which is still correlated with the formation energy of the tetragonal configuration. So, as the simplest approximation, we consider all the energies of the cubic systems with $x < 1/3$ to be shifted upwards with respect to the tetragonal ones by $\Delta E$, plus an additional shift depending linearly on the configuration's Co concentration only. In other words,

$$ E_n^{cub} = E_n^{tet} + \Delta E + cx_n, \qquad (6) $$

where $E_n^{cub}$ is the formation energy of a specific configuration $n$ in the cubic phase, $E_n^{tet}$ is the formation energy of the corresponding configuration in the tetragonal phase, and $c$ is the only adjustable parameter in the model. This parameter would not be necessary if sufficient ab initio data were available for the cubic phase. Also, it is possible that a slight composition-dependence of the vibrational free energy may be accounted for by this term. We find that $c = 0.021$ eV per formula unit per Co ion in the supercell is a suitable value. Since the formation energy of the cubic structure becomes zero on the Co$_3$O$_4$ side of the phase diagram, its variation across the supercell's Co-composition range of 0-24 Co atoms is expected to be of the order of $\Delta E/24$. Using $\Delta E = 3.8$ eV/supercell $= 0.475$ eV/(formula unit), implies a rough estimate of $c \approx 0.02$ eV per formula unit per Co. In the Appendix we show the sensitivity of the fit on the value of this parameter.

Therefore, from Eqs. (5) and (6),

$$ F^{cub}(T, x) \simeq F^{tet}(T, x) + \Delta E \left\{ 1 - \frac{T}{T_C} \right\} - cx. \qquad (7) $$
<figure><img src="image_5.png" /><figcaption>FIG. 4. (a) Free energy of (Co$_x$Mn$_{1-x}$)$_3$O$_4$ vs Co concentration $x$ calculated when including configurational, paramagnetic, and vibrational entropy after correcting the configurational entropy with tangents indicated by black lines. (b) Corresponding solubility gap before correcting the configurational entropy (dashed lines) and after (solid lines), with experimental data indicated by the shaded region.</figcaption></figure>

The resulting SG in the $0 < x < 1/3$ composition range (close to the H + S phase coexistence region), is shown in Fig. 4(b) (dashed lines, labeled “Uncorrected model”). The predicted SG is now beginning to resemble much more closely the experimental SG, demonstrating the importance of phonon contributions when calculating phase diagrams, which has also been noted in other works [50,51]. However, we still see substantial disagreement as $x \to 0$. Due to the incorrect temperature dependence, as the left-hand side tends to $x = 0$ too early, the algorithm fails to find a solution of Eq. (B3) for the right-hand side of the gap above $T = 1420$ K.

## IV. CORRECTION FOR THE LOW-CONCENTRATION EXTREME

The free energy calculated directly by the current model yields an unphysical SG behavior that does not end in a single point at $x = 0$, $T = T_C$, as in the experimental data also shown in Fig. 4(b). This is due to the use of finite supercells in the model, because it is not possible to represent systems with these very low concentrations of the substituting species or to represent all possible configurations at a given low concentration. This problem is inherent to the direct calculation of the partition function from ab initio configuration energies, and it would not appear if the configurational entropy is approximated as an analytical nonideal solution model, as often done in Calphad methods and some ab initio-based approaches [52,53]. We therefore proceed to develop an ana-

lytical correction to the configurational free energy from our finite-supercell model, $S_{\text{config}}^{\text{model}} T$, in the dilute limit.

To determine an analytical correction we return to the ideal-solution model, which we initially disregarded as unsuitable for this system, just as a starting point. The advantage of using the ideal-solution model is that we can easily consider very large system sizes, much more easily than performing ab initio calculations for infeasibly large supercells. In the dilute limit we can consider substituting species as being sufficiently separated such that they do not interact with each other and that the configurations are energetically equivalent, the only difference being if the substituting species occupies an oh- or td-coordinated lattice site.

To begin with, we write the general expression for the configurational entropy from the ideal-solution approximation in terms of the concentration of Co, $x = \frac{n_{\text{Co}}}{N_{\text{tot}}}$, for the 56-atom supercell,

$$S_{\text{supercell}}^{\text{ideal}} = k_{\text{B}} \ln \left( \frac{N_{\text{tot}}!}{(xN_{\text{tot}})!(N_{\text{tot}} - xN_{\text{tot}})!} \right), \quad (8)$$

where $N_{\text{tot}}$ is the total number of TM lattice sites that the substituting species can occupy. In the 56-atom supercell, $N_{\text{tot}}$ is 24 sites, of which 8 are td- and 16 are oh-coordinated. For the case of large $N_{\text{tot}}$, Eq. (8) can be rewritten using the Stirling approximation as [22]

$$S_{\text{Stirling}}^{\text{ideal}} = -N_{\text{tot}}k_{\text{B}}[(1-x)\ln(1-x) + x\ln x]. \quad (9)$$

When considering $S_{\text{Stirling}}^{\text{ideal}}$ per site, the expression is independent of $N_{\text{tot}}$ and we show in the SM (Sec. 3) that, when increasing $N_{\text{tot}}$ in Eq. (8), this tends towards the value obtained with Eq. (9). We therefore use Eq. (9) in our analytical correction to obtain an approximation to the configurational entropy for an infinite supercell [44].

The ideal-solution model expression, as it stands, assumes that all $N_{\text{tot}}$ TM sites are energetically equivalent and hence all are accessible. However, from the ab initio data of Ref. [4], we know that, at $T = 0$ K, the ground state in the low-concentration limit is to occupy the td sites with Co first, because these are the lower-energy structures. This would suggest that the number of accessible sites in the supercell should be eight, but at finite temperatures there will be some occupation probability of the oh sites. Across the temperature range of interest for this study (1200–1443 K), the number of accessible sites in the supercell should lie somewhere between these two extreme cases, i.e., between 8 and 24.

To obtain a suitable “modified-ideal-solution” expression for the model, we replace $N_{\text{tot}}$ in Eq. (8) by $\alpha N_{\text{tot}}$, where $\alpha$ is an adjustable parameter between $\frac{1}{3}$ and 1, which would scale the total number of accessible sites in the supercell to be between 8 and 24. To determine a suitable value of $\alpha$, we compare $S_{\text{config}}^{\text{ideal}}$ to the true expression for $S_{\text{config}}^{\text{true}}$ which we obtained by explicitly calculating the partition function $Q$ and using Eq. (1) to compute $F_{\text{config}}$. We then calculate the configurational entropy term via [24],

$$S_{\text{config}} = - \left( \frac{\partial F_{\text{config}}}{\partial T} \right)_V . \quad (10)$$
From comparing the curves for $S_{\text{config}}^{\text{ideal}}$ T and $S_{\text{config}}^{\text{true}}$ T against T, we find that a value of $\alpha = 0.6$ (which is the equivalent of 14.4 accessible sites in the 24-atom supercell) results in good agreement over the temperature range of interest in this study. Further details of this procedure are provided in Secs. 3 of the SM along with a plot showing the sensitivity of the calculated SG on the $\alpha$ parameter [44].

Thus, the correction to $S_{\text{config}}^{\text{model}}$, $\Delta S(T, x)$ is

$$ \Delta S(T, x) = \left( \frac{S_{\text{Stirling}}^{\text{ideal}}(x, \alpha)}{S_{\text{supercell}}^{\text{ideal}}(x, \alpha, N_{\text{tot}})} - 1 \right) S_{\text{config}}^{\text{model}}(T, x), \quad (11) $$

where $N_{\text{tot}}$ has been replaced by $\alpha N_{\text{tot}}$ in $S_{\text{supercell}}^{\text{ideal}}$. We then correct the total free energy from our model, $F_{\text{tot}}^{\text{model}}$ (which considered not only configurational but also magnetic and vibrational entropy) as

$$ F_{\text{tot}}^{\text{corr}}(T, x) = F_{\text{tot}}^{\text{model}}(T, x) - T \Delta S(T, x). \quad (12) $$

A plot of $F_{\text{tot}}^{\text{corr}}$ as a function of $x$ at various temperatures is shown in Fig. 4(a). This corrected free energy yields the SG shown by the solid line in Fig. 4(b). The extremal dependence has been altered the most noticeably with the intermediate $x$ range less impacted by the correction (see details in Sec. 3 of the SM). After applying this correction the predicted SG now resembles much more closely the experimental data, in particular, the spindle shape as $x \to 0$. This shows that it is possible to predictively calculate the SG for the mixed TM oxides.

## V. IMPACT OF MAGNETISM REVISITED

After having developed the methods to account for vibrational entropy in the model and for correcting the configurational entropy as $x \to 0$, it is worth reexamining the role of paramagnetic entropy in the final model. As shown in Fig. 1, the magnetic degrees of freedom had a dramatic effect on the free-energy curve but the predicted SG was still far from that measured experimentally. Now, if we only take into account configurational entropy, vibrational entropy and the correction to the configurational entropy, but completely neglect magnetism, we obtain the F curve shown in Fig. 5(a), which results in the calculated SG shown in Fig. 5(b). The c parameter used in the vibrational entropy calculation here was chosen to be the one that minimized the root mean squared error (RMSE) of the predicted SG relative to the experimental data, which in this case was 0.014. The sensitivity of the model prediction on the c parameter is discussed in the Appendix.

The minimum RMSE for the SG calculated without magnetism is 0.0338, compared with 0.0096 for the corrected curve shown in Fig. 4(b). Furthermore, without including paramagnetic entropy in the model we are only able to find solutions for the SG in the temperature range of 1320 to 1440 K. Upon comparing Fig. 4 with Fig. 5, it is clear that the shape of the F curve has changed dramatically by the inclusion of paramagnetic entropy. The contribution of paramagnetic entropy lowers the free energy across the full composition range but lowers it the most dramatically at lower Co concentrations due to the larger atomic magnetic moments of Mn ions. The influence of paramagnetic entropy produces a

<figure><img src="image_6.png" /><figcaption>FIG. 5. (a) The Helmholtz free energy F calculated with a model neglecting paramagnetic entropy but including configurational entropy, vibrational entropy, and a correction to the configurational entropy at the composition extreme at temperatures ranging from 1320 to 1440 K. (b) The solubility gap (SG) calculated with F from panel (a).</figcaption></figure>

much more convex F curve. The crucial role of paramagnetic entropy was less clear in Fig. 1, but with our final model we see that inclusion of this contribution to the entropy is vital to accurately predict the SG for this system.

## VI. CONCLUSIONS

In summary, we have shown how to accurately predict the SG of TM oxides, and systematically demonstrated the relative importance of different contributions to the free energy in the case of $(\text{Co}_x\text{Mn}_{1-x})_3\text{O}_4$ when calculating its SG from *ab initio* data. First, vibrational entropy plays a central role in determining the SG. If it is neglected, the calculated SG is unable to achieve even qualitative agreement with experimental measurements. Second, for a spindle-shaped SG, it is necessary to correct for finite-size limitations when the model is derived from *ab initio* calculations employing finite supercells. Not doing this leads to a wrong prediction of the temperature dependence of the SG for compositions close to those of the pure phases. Lastly, for this system the inclusion of paramagnetic entropy is crucial to obtain the SG in the right range of temperatures and compositions. The neglect of this contribution produces a SG in an incorrect position and with an incorrect width. The effect of paramagnetic entropy is particularly pronounced for $(\text{Co}_x\text{Mn}_{1-x})_3\text{O}_4$. This is due to the composition dependence of the number of magnetic degrees of freedom in this system, where the type of TM
occupying a particular type of crystallographic site influences the atomic magnetic moment.

Software developed in this work for calculating the solubility gap of (Co$_x$Mn$_{1-x}$)$_3$O$_4$ is available under a BSD 3-Clause license [55].

ACKNOWLEDGMENTS

We thank Stefano Curtarolo and Jesús Carrete for useful comments on the paper. Work at CEA-Grenoble was supported by Institut Carnot, through projects MAPPE and PREDICT. This work was performed using HPC resources from GENCI-TGCC Grant A0060910765. J.C. acknowledges the computer resources at SGI/IZO-SGIker UPV/EHU (Arina cluster), the i2BASQUE academic network, and MareNostrum and the technical support provided by the Barcelona Supercomputer Center (QS-2020-2-0002). Work at CIC energiGUNE was supported by the Ministerio de Ciencia e Innovación of Spain through the project ION-SELF (No. PID2019-106519RB-I00). A.S. was supported by the Russian Foundation for Basic Research under Grant No. 20-53-12012

APPENDIX A: Crystal STRUCTURES OF Co$_3$O$_4$ AND Mn$_3$O$_4$

Visualizations of the crystal structures are given in Fig. 6.

APPENDIX B: CALCULATION OF A SOLUBILITY GAP

For solids at atmospheric pressure the product PV is generally negligible [56], therefore the Gibbs free energy G of each phase can be approximated by the Helmholtz free energy F at its minimizing volume [42]. At a given temperature T,

<figure><img src="image_7.png" /><figcaption>FIG. 6. (a), (c) Mn$_3$O$_4$ and (b), (d) Co$_3$O$_4$ with O in red, Co in blue, and Mn in purple. Tetrahedrally coordinated (td) and octahedrally coordinated (oh) sites are indicated in panel (c). Figures generated using VESTA [54].</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>FIG. 7. Illustration of the tangent rule to determine the phase coexistence region between Hausmannite (H) and cubic spinel (S) phases. The dashed blue line is tangent to the free-energy curves of the cubic, F$_{cub}$, and tetragonal, F$_{tet}$, phases at concentrations x$_l(T)$ and x$_r(T)$, respectively. These concentrations mark the limits of the phase coexistence region as a function of T (1260 K in the figure). As T increases, the cubic free energy (pink line) gets lower with respect to the tetragonal one (orange line), eventually making the coexistence region vanish at T = 1443 K.</figcaption></figure>

a homogeneous system of single concentration x$_0$ can lower its free energy by separating into domains of two different concentrations if there exist two temperature-dependent compositions x$_l$ < x$_0$ and x$_r$ > x$_0$, where subscripts l and r denote left and right sides of the SG, respectively, such that

$$F(T, x_0) > c_r F(T, x_r) + c_l F(T, x_l), \quad (B1)$$

where

$$c_r = \frac{x_0 - x_l}{x_r - x_l} \quad \text{and} \quad c_l = \frac{x_r - x_0}{x_r - x_l}. \quad (B2)$$

are the relative fractions of material with concentration x$_r$ and x$_l$ respectively. As demonstrated graphically in Fig. 7, this means that the line tangent to F(T, x) at x$_l$ and x$_r$ is below F(T, x) for x = x$_0$. The SG at a given T is the range of concentrations between x$_l$ and x$_r$. These two points can be obtained by numerically solving the two coupled differential equations,

$$\left. \frac{\partial F}{\partial x} \right|_{x_l} = \left. \frac{\partial F}{\partial x} \right|_{x_r} = \frac{F(x_l) - F(x_r)}{x_l - x_r}. \quad (B3)$$
<figure><img src="image_9.png" /><figcaption>FIG. 8. (a) Solubility gap (SG) calculated by models including configurational, paramagnetic, and vibrational entropy with a correction to the configurational entropy at the composition extreme with different values for the *c* parameter when calculating the vibrational entropy. (b)–(d) Corresponding Helmholtz free energies calculated with different values of *c* across a temperature range of 1200–1440 K at every 1 K increment.</figcaption></figure>

## APPENDIX C: SENSITIVITY OF MODEL PREDICTIONS TO THE *c* PARAMETER

As discussed in the main text, we approximate the free energy of the cubic phase of $(Co_xMn_{1-x})_3O_4$, $F^{cub}$, using the difference in the formation energy of the cubic and tetragonal phases of $Mn_3O_4$, $\Delta E$, and the free energy of the tetragonal phase of the corresponding configuration of $(Co_xMn_{1-x})_3O_4$, $F^{tet}$, using Eq. (7). Here we discuss the sensitivity of our model on the choice of the parameter *c* in Eq. (7).

We found that $c = 0.021$ reproduces well the experimental solubility gap (SG) for $(Co_xMn_{1-x})_3O_4$. Figure 8(a) shows the SGs calculated using models accounting for configurational, paramagnetic, and vibrational entropy with a correction to the configurational entropy at the composition extreme but with different values for the *c* parameter when calculating the vibrational entropy. Figure 8(a) shows that a larger *c* shifts the SG to the left and broadens the gap. Figures 8(b)–8(d) show the corresponding free-energy curves calculated at temperatures 1200 to 1440 K at every 1 K increment. The dark regions in each plot indicate the tangent of each curve used to calculate the SG. Here it can also be seen that the dark region shifts to the left and broadens as *c* increases. We can also see in the free-energy curves that the kink in the curve, due to the transition between the tetragonal and cubic phase, is more pronounced with a larger *c*.

The observations discussed above can be understood using Fig. 9. The figure shows that the crossover point between the tetragonal and cubic phase occurs furthest to the left (at the lowest number of Co atoms in the supercell) when *c* is the largest value (crossovers between the solid and dotted lines). We can also see that the kink is more pronounced and

narrow at this crossover than that when *c* is the lowest value (crossovers between the solid and dashed lines).

## APPENDIX D: CONFIGURATIONAL ENTROPY EXTRACTED FROM THE MODEL VS IDEAL-SOLUTION APPROXIMATION

Figure 10 shows the configurational entropy for the 56-atom (24 substitutional TM sites) supercell from the ideal-solution model against that extracted from the DFT-based model at various temperatures. At the intermediate concentrations the ideal-solution approximation overestimates the

<figure><img src="image_10.png" /><figcaption>FIG. 9. Helmholtz free energy *F* of $(Co_xMn_{1-x})_3O_4$ calculated at *T* = 1300 K (blue lines) and *T* = 1400 K (orange lines) where solid lines are for the tetragonal phase, and the dotted and dashed lines are for the cubic phase with free energy calculated using different values of the parameter *c*.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>FIG. 10. Configurational entropy as a function of concentration of Co in (CoₓMn₁₋ₓ)₃O₄ extracted from our model at various temperatures plotted against the ideal-solution approximation for an equivalent supercell.</figcaption></figure>

configurational entropy relative to the DFT-based model and that the overestimation is not systematic. There are regions where configurations are closer in energy (the temperature dependence of the configurational entropy is reduced) and regions where this is more pronounced. This is likely to be due to composition-dependent energetic variations in the configurations which are likely to be related to site occupations (oh or td sites) and to chemical interactions, both of which are not accounted for in the ideal-solution approximation.

## APPENDIX E: JUSTIFICATION FOR NEGLECTING THE ELECTRONIC ENTROPY OF (CoₓMn₁₋ₓ)₃O₄

The EDOS for various relaxed configurations of (CoₓMn₁₋ₓ)₃O₄ from the dataset of Ref. [4] are calculated using VASP [57,58], a 12 × 12 × 12 k grid, and the SCAN meta-GGA functional [59], which has been demonstrated to give better band-gap predictions than GGA and overcomes the issue of subjectivity when selecting the value of U for GGA + U. However, this method does still typically underestimate experimental band gaps [60]. The calculated EDOS shows variable fundamental band gaps ranging from clearly semiconducting or insulating (with band gaps in excess of 1.5 e V) to semimetallic. Noting that we can still expect some underestimation of the band gap with the SCAN meta-GGA functional and further reduction in electronic conductivity may arise from electronic correlations, which are not accounted for in our predictions. Plots of all calculated EDOS are included in Sec. 4 of the SI [44].

In all cases the DOS close to the Fermi level was typically low. Of all the calculated EDOS, Fig. 11 shows the case with the highest number of states close the Fermi level. The peak number of states close to the Fermi level can be read off from Fig. 11 as approximately 80 states per eV per supercell. This value is used to approximate an overestimated value for the electronic entropy for all compositions of the mixed phase.

The electronic entropy $S_{elec}$ defined from the occupation probabilities of single electron states [12],

$$f(E) = \frac{1}{1 + e^{\frac{E-E_F}{k_B T}}}, \qquad (E1)$$

<figure><img src="image_12.png" /><figcaption>FIG. 11. Calculated electronic density of states for (CoₓMn₁₋ₓ)₃O₄. In this case, the second minimum-energy structure from the training data of all structures with 20 Co in the supercell.</figcaption></figure>

is given by

$$S_{elec} = -k_B \int D(E) \{f(E) \ln f(E) + [1 - f(E)] \ln[1 - f(E)] \} dE, \qquad (E2)$$

where $D(E)$ is the electronic density of states and $E_F$ is the Fermi level. The function inside the curly parentheses in Eq. (E2) is only significantly greater than zero in an interval of less than $E_F - 10k_B T < E < E_F + 10k_B T$. A plot of this function is shown in Fig. 12.

Taking this as the relevant energy interval around $E_F$, an upper bound for the number of accessible electronic states is calculated at the highest temperature considered in this work ($T = 1400$ K) and by assuming an infinitely dense $D(E)$ as a rectangle of height 80 states per eV per supercell (the maximum number close to $E_F$ for all calculated EDOS) and width of $2 \times 10k_B T$ eV. This gives a maximum range of electronic states that could contribute to electronic entropy as 0 to 207

<figure><img src="image_13.png" /><figcaption>FIG. 12. Function taken from the expression for electronic entropy $S_{elec}$ in Eq. (E2) (y axis) to show the range of energies above or below the Fermi energy $E_F$ in units of $k_B T$ where the value of the function is significantly larger than zero.</figcaption></figure>
per supercell. To put this range into context, the minimum and maximum number of magnetic states per supercell [obtained using Eq. (3)] across all possible configurations are 5 764 801 and 3.972 × 10²³, respectively. Expressing entropy in terms of accessible states of the system, w,

$$S = k_B \ln(w), \qquad (E3)$$

and considering all possible accessible states contributing to the total w, the electronic states (and any variation in the number between different configurations) will clearly be a much smaller contribution compared with magnetic states and variations in the number of magnetic states between different compositions in the case of (CoₓMn₁₋ₓ)₃O₄.

[1] L. André, S. Abanades, and L. Cassayre, Mixed metal oxide systems applied to thermochemical storage of solar energy: Benefits of secondary metal addition in Co and Mn oxides and contribution of thermodynamics, Appl. Sci. 8, 2618 (2018).

[2] Y. Wang, J. Li, and Z. Wei, Transition-metal-oxide-based catalysts for the oxygen reduction reaction, J. Mater. Chem. A 6, 8194 (2018).

[3] R. Jacob, M. Belusko, A. I. Fernández, L. F. Cabeza, W. Saman, and F. Bruno, Embodied energy and cost of high temperature thermal energy storage systems for use with concentrated solar power plants, Appl. Energy 180, 586 (2016).

[4] A. Zaki, J. Carrasco, D. Bielsa, and A. Faik, Tunable redox temperature of a Co₃₋ₓMnₓO₄ (0 ≤ x ≤ 3) continuous solid solution for thermochemical energy storage, ACS Appl. Mater. Interfaces 12, 7010 (2020).

[5] C. Agrafiotis, M. Roeb, M. Schmücker, and C. Sattler, Exploitation of thermochemical cycles based on solid oxide redox systems for thermochemical storage of solar heat. Part 1: Testing of cobalt oxide-based powders, Sol. Energy 102, 189 (2014).

[6] E. Alonso, A. Gallo, C. Pérez-Rábago, and E. Fuentealba, Thermodynamic study of CuO/Cu₂O and Co₃O₄/CoO redox pairs for solar energy thermochemical storage, in SOLARPACES 2015: International Conference on Concentrating Solar Power and Chemical Energy Systems, edited by V. Rajpaul and C. Richter, AIP Conf. Proc. No. 1734 (AIP, New York, 2015), p. 050004.

[7] C. Agrafiotis, M. Roeb, and C. Sattler, Exploitation of thermochemical cycles based on solid oxide redox systems for thermochemical storage of solar heat. Part 4: Screening of oxides for use in cascaded thermochemical storage concepts, Sol. Energy 139, 695 (2016).

[8] A. J. Carrillo, D. P. Serrano, P. Pizarro, and J. M. Coronado, Thermochemical heat storage based on the Mn₂O₃/Mn₃O₄ redox couple: Influence of the initial particle size on the morphological evolution and cyclability, J. Mater. Chem. A 2, 19435 (2014).

[9] A. Bayon, R. Bader, M. Jafarian, L. Fedunik-Hofman, Y. Sun, J. Hinkley, S. Miller, and W. Lipiński, Techno-economic assessment of solid-gas thermochemical energy storage systems for solar thermal power applications, Energy (Oxford, UK) 149, 473 (2018).

[10] L. André, S. Abanades, and L. Cassayre, Experimental investigation of Co–Cu, Mn–Co, and Mn–Cu redox materials applied to solar thermochemical energy storage, ACS Appl. Energy Mater. 1, 3385 (2018).

[11] P. Spencer, A brief history of CALPHAD, CALPHAD: Comput. Coupling Phase Diagrams Thermochem. 32, 1 (2008).

[12] M. Widom, Modeling the structure and thermodynamics of high-entropy alloys, J. Mater. Res. 33, 2881 (2018).

[13] A. van de Walle, Methods for first-principles alloy thermodynamics, JOM 65, 1523 (2013).

[14] A. Van der Ven, J. Thomas, B. Puchala, and A. Natarajan, First-principles statistical mechanics of multicomponent crystals, Annu. Rev. Mater. Res. 48, 27 (2018).

[15] A. R. Natarajan and A. Van der Ven, Machine-learning the configurational energy of multicomponent crystalline solids, npj Comput. Mater. 4, 56 (2018).

[16] M. C. Kemei, J. K. Harada, R. Seshadri, and M. R. Suchomel, Structural change and phase coexistence upon magnetic ordering in the magnetodielectric spinel Mn₃O₄, Phys. Rev. B 90, 064418 (2014).

[17] C. Haas, Phase transitions in crystals with the spinel structure, J. Phys. Chem. Solids 26, 1225 (1965).

[18] R. Arras, T. L. Le, S. Guillemet-Fritsch, P. Dufour, and C. Tenailleau, First-principles electronic structure calculations for the whole spinel oxide solid solution range MnₓCo₃₋ₓO₄ (0 ≤ x ≤ 3) and their comparison with experimental data, Phys. Chem. Chem. Phys. 18, 26166 (2016).

[19] S. Guillemet-Fritsch, C. Tenailleau, H. Bordeneuve, and A. Rousset, Magnetic properties of cobalt and manganese oxide spinel ceramics, Adv. Sci. Technol. (Durnten-Zurich, Switz.) 67, 143 (2010).

[20] Y. Golikov, S. Tubin, V. Barkhatov, and V. Balakirev, Phase diagrams of the Co-Mn-O system in air, J. Phys. Chem. Solids 46, 539 (1985).

[21] A. I. Kitaigorodsky, Mixed Crystals, edited by M. Cardona, P. Fulde, and H.-J. Queisser (Springer, Berlin, Heidelberg, 1984), Vol. 33.

[22] A.-B. Chen and A. Sher, Semiconductor Alloys (Plenum Press, New York, 1995).

[23] A. van de Walle, G. Ghosh, and M. Asta, Ab initio modeling of alloy phase equilibria, in Applied Computational Materials Modeling: Theory, Simulation and Experiment, edited by G. Bozzolo, R. D. Noebe, P. B. Abel, and D. Vij (Springer US, Boston, 2007), pp. 1–34.

[24] L. D. Landau and E. M. Lifshitz, Statistical Physics, 3rd ed. (Butterworth-Heinemann, Amsterdam, London, 1980), Vol. 5, Part 1.

[25] C. M. Rost, E. Sachet, T. Borman, A. Moballegh, E. C. Dickey, D. Hou, J. L. Jones, S. Curtarolo, and J.-P. Maria, Entropy-stabilized oxides, Nat. Commun. 6, 8485 (2015).

[26] A. V. Shapeev, Moment tensor potentials: A class of systematically improvable interatomic potentials, Multiscale Model. Simul. 14, 1153 (2016).

[27] I. S. Novikov, K. Gubaev, E. V. Podryabinkin, and A. V. Shapeev, The MLIP package: Moment tensor potentials with
MPI and active learning, Mach. Learn.: Sci. Technol. **2**, 025002 (2021).

[28] S. K. Wallace, A. S. Bochkarev, A. van Roekeghem, J. Carrasco, A. Shapeev, and N. Mingo, Free energy of CoₓMn₃₋ₓO₄ mixed phases from machine-learning-enhanced *ab initio* calculations (unpublished).

[29] N. Ashcroft and N. Mermin, Diamagnetism and paramagnetism, Solid State Physics (Saunders College Publishing, 1976).

[30] Y. Ikeda, B. Grabowski, and F. Körmann, Ab initio phase stabilities and mechanical properties of multicomponent alloys: A comprehensive review for high entropy alloys and compositionally complex alloys, Mater. Charact. **147**, 464 (2019).

[31] U. G. Nielsen, I. Heinmaa, A. Samoson, J. Majzlan, and C. P. Grey, Insight into the local magnetic environments and deuteron mobility in jarosite (AFe₃(SO₄)₂(OD, OD₂)₆, A = K, Na, D₃O) and hydronium alunite ((D₃O)Al₃(SO₄)₂(OD)₆), from variable-temperature ²H MAS NMR spectroscopy, Chem. Mater. **23**, 3176 (2011).

[32] V. Brabers and A. Broemme, Low-spin-high-spin transition in the Co₃O₄ spinel, J. Magn. Magn. Mater. **104-107**, 405 (1992).

[33] H. C. O'Neill, Thermodynamics of Co₃O₄: A possible electron spin unpairing transition in Co₃, Phys. Chem. Miner. **12**, 149 (1985).

[34] K. Mocala, A. Navrotsky, and D. Sherman, High-temperature heat capacity of Co₃O₄ spinel: Thermally induced spin unpairing transition, Phys. Chem. Miner. **19**, 88 (1992).

[35] M. Chen, B. Hallstedt, and L. J. Gauckler, Thermodynamic assessment of the Co-O system, J. Phase Equilib. **24**, 212 (2003).

[36] N. Shulumba, B. Alling, O. Hellman, E. Mozafari, P. Steneteg, M. Odén, and I. A. Abrikosov, Vibrational free energy and phase stability of paramagnetic and antiferromagnetic CrN from *ab initio* molecular dynamics, Phys. Rev. B **89**, 174108 (2014).

[37] G. Salek, P. Dufour, S. Guillemet-Fritsch, and C. Tenailleau, Sustainable low temperature preparation of Mn₃₋ₓCoₓO₄ (0 ≤ x < 3) spinel oxide colloidal dispersions used for solar absorber thin films Mater. Chem. Phys. **162**, 252 (2015).

[38] N. Shulumba, O. Hellman, Z. Raza, B. Alling, J. Barrirero, F. Mücklich, I. A. Abrikosov, and M. Odén, Lattice Vibrations Change the Solid Solubility of an Alloy at High Temperatures, Phys. Rev. Lett. **117**, 205502 (2016).

[39] B. P. Burton and A. van de Walle, First-principles phase diagram calculations for the system NaCl-KCl: The role of excess vibrational entropy, Chem. Geol. **225**, 222 (2006).

[40] C. Wolverton and V. Ozoliņš, Entropically Favored Ordering: The Metallurgy of Al₂Cu Revisited, Phys. Rev. Lett. **86**, 5518 (2001).

[41] C. Ravi, C. Wolverton, and V. Ozoliņš, Predicting metastable phase boundaries in Al–Cu alloys from first-principles calculations of free energies: The role of atomic vibrations, Europhys. Lett. **73**, 719 (2006).

[42] R. Feng, P. K. Liaw, M. C. Gao, and M. Widom, First-principles prediction of high-entropy-alloy stability, npj Comput. Mater. **3**, 50 (2017).

[43] A. Togo and I. Tanaka, First principles phonon calculations in materials science, *Scr. Mater.* **108**, 1 (2015).

[44] See Supplemental Material at http://link.aps.org/supplemental/10.1103/PhysRevResearch.3.013139 for details on the calculations and additional graphs.

[45] A. van Roekeghem, J. Carrete, and N. Mingo, Anomalous thermal conductivity and suppression of negative thermal expansion in ScF₃, Phys. Rev. B **94**, 020303(R) (2016).

[46] A. van Roekeghem, J. Carrete, C. Oses, S. Curtarolo, and N. Mingo, High-throughput computation of thermal conductivity of high-temperature solid phases: The case of oxide and fluoride perovskites, Phys. Rev. X **6**, 041061 (2016).

[47] A. van Roekeghem, J. Carrete, and N. Mingo, Quantum self-consistent ab-initio lattice dynamics, arXiv:2006.12867.

[48] P. Souvatzis, O. Eriksson, M. I. Katsnelson, and S. P. Rudin, Entropy Driven Stabilization of Energetically Unstable Crystal Structures Explained from First Principles Theory, Phys. Rev. Lett. **100**, 095901 (2008).

[49] T. Tadano and S. Tsuneyuki, Self-consistent phonon calculations of lattice dynamical properties in cubic SrTiO₃ with first-principles anharmonic force constants, Phys. Rev. B **92**, 054301 (2015).

[50] R. Chinnappan, B. Panigrahi, and A. van de Walle, First-principles study of phase equilibrium in Ti–V, Ti–Nb, and Ti–Ta alloys, CALPHAD: Comput. Coupling Phase Diagrams Thermochem. **54**, 125 (2016).

[51] V. Ozoliņš and M. Asta, Large Vibrational Effects Upon Calculated Phase Boundaries in Al-Sc, Phys. Rev. Lett. **86**, 448 (2001).

[52] M. Y. Lavrentiev, N. L. Allan, G. D. Barrera, and J. A. Purton, Ab initio calculation of phase diagrams of oxides, J. Phys. Chem. B **105**, 3594 (2001).

[53] N. Saunders, Phase diagram calculations for high-temperature structural materials, in *High-Temperature Structural Materials*, edited by R. W. Cahn, A. G. Evans, and M. McLean (Springer, Dordrecht, 1996), pp. 125–143.

[54] K. Momma and F. Izumi, VESTA: A three-dimensional visualization system for electronic and structural analysis, J. Appl. Crystallogr. **41**, 653 (2008).

[55] N. Mingo and S. K. Wallace, Zenodo, version v.1, 2020, doi: 10.5281/zenodo.4133862.

[56] A. van de Walle and G. Ceder, The effect of lattice vibrations on substitutional alloy thermodynamics, Rev. Mod. Phys. **74**, 11 (2002).

[57] G. Kresse and J. Furthmüller, Efficiency of ab-initio total energy calculations for metals and semiconductors using a plane-wave basis set, Comput. Mater. Sci. **6**, 15 (1996).

[58] G. Kresse and J. Furthmüller, Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set, Phys. Rev. B **54**, 11169 (1996).

[59] J. Sun, A. Ruzsinszky, and J. P. Perdew, Strongly Constrained and Appropriately Normed Semilocal Density Functional, Phys. Rev. Lett. **115**, 036402 (2015).

[60] Z. H. Yang, H. Peng, J. Sun, and J. P. Perdew, More realistic band gaps from meta-generalized gradient approximations: Only in a generalized Kohn-Sham scheme, Phys. Rev. B **93**, 205205 (2016).