Towards Ultra Low Cobalt Cathodes: A High Fidelity Computational Phase Search of Layered Li-Ni-Mn-Co Oxides

Layered Li(Ni,Mn,Co,)O$_2$ (NMC) presents an intriguing ternary alloy design space for optimization as a cathode material in Li-ion batteries. Recently, the high cost and resource limitations of Co have added a new design constraint and high Ni-containing NMC alloys have gained enormous attention despite possible performance trade-offs. It is not fully understood if this material space is a disordered solid solution at room temperature and any arbitrary combination can be used or if there exist distinct transition metal orderings to which meta-stable solid solutions will decay during cycling and affect performance. Here, we present a high fidelity computational search of the ternary phase diagram with an emphasis on high-Ni, and thus low Co, containing compositional phases to understand the room temperature stability of the ordered and disordered solid solution phases. This is done through the use of density functional theory training data fed into a reduced order model Hamiltonian that accounts for effective electronic and spin interactions of neighboring transition metal atoms at various lengths in a background of fixed composition and position lithium and oxygen atoms. This model can then be solved to include finite temperature thermodynamics into a convex hull analysis to understand the regions of ordered and disordered solid solution as well the transition metal orderings within the ordered region of the phase diagram. We find that for the majority of transition metal compositions of the layered material, specifically medium to high-Ni content, prefer transition metal ordering and predict the collection of preferred compositions in the ordered region.

at room temperature and any arbitrary combination can be used or if there exist distinct transition metal orderings to which meta-stable solid solutions will decay during cycling and affect performance. Here, we present a high fidelity computational search of the ternary phase diagram with an emphasis on high-Ni, and thus low Co, containing compositional phases to understand the room temperature stability of the ordered and disordered solid solution phases. This is done through the use of density functional theory training data fed into a reduced order model Hamiltonian that accounts for effective electronic and spin interactions of neighboring transition metal atoms at various lengths in a background of fixed composition and position lithium and oxygen atoms. This model can then be solved to include finite temperature thermodynamics into a convex hull analysis to understand the regions of ordered and disordered solid solution as well the transition metal orderings within the ordered region of the phase diagram. We also provide a method to propagate the uncertainty at every level of the analysis to the final prediction of thermodynamically favorable compositional phases thus providing a quantitative measure of confidence for each prediction made. Due to

Introduction
The first rechargeable Li-ion battery cathode material, layered LiCoO 2 , 1 revolutionized the world of portable electronics. The high cost and low operational capacity of this battery material was improved with the addition of Ni and Mn to create LiNi x Mn y Co (1−x−y) O 2 (NMC).
This was first successfully done by Lui et al. 2 and later the most widespread phase, x=y=1/3 (NMC111) was synthesized by Ohzuku et al. 3 Since then a small collection of phases have been used in various applications with some of the most promising candidates being x=0.5, y=0.3 (NMC532), 4-6 and x=0.8,y=0.1 (NMC811). 7,8 Additionally, NMC cathode materials are now present in a majority of electric vehicle battery chemistries. 9 Many studies have explored trends of operational performance as a function of Ni, Mn, and Co content. Julien et al. demonstrated that with increasing Co content, the specific capacity increases. 10 Manthiram et al on the other hand demonstrated a capacity increase with increasing Ni content. 11 The discrepancy between these studies can be attributed to their distinct, one-dimensional sweeps of the full two-dimensional phase space. Another consideration regarding transition metal composition is that of thermal stability. Recent work has shown that thermal stability decreases with increasing Ni content. 12 In this study, however, no phases were tested for Ni content between 60% and 80% and only one possible combination of Co and Mn content was used for each Ni content. Only with a full understanding of the entire compositional phase space, can trends such as these be understood and possibly broken. Additionally within all of these studies, the material space is assumed to be a disordered solid solution. While this is likely true at synthesis temperatures given the synthesizability of NMC at arbitrary compositions, this is likely not true at ambient temperatures at which the material is electrochemically cycled within a battery. In fact a tendency for ordering has been seen even for NMC materials that were rapidly quenched after synthesis. [13][14][15][16][17] Additionally an influence of this ordering on aspects such as the voltage profile during cycling has also been seen. 14 Thus the proper ordering of the transition metals is necessary for proper simulation of the material. Many previous computational works have provided reasonable guesses for various other phases with, however, no exploration for the stability with respect to other cation orderings at the same composition. 5,18,19 Within this work, we also specifically attempt to address the need for new low-Co content cathode materials as Co production presents major price volatility and supply concerns in the near future. 20 Amid these concerns of material availability, the United States Department of Energy has identified the need for low-Co cathode materials as imperative for future battery technologies, setting a target of 50 mg W h of Co. 21 This target however is far from realization in materials that can meet the performance needs currently satisfied by high-Co cathodes.
We attempt to understand if the region of ultra low Co is a disordered solid solution where any composition can be expected to have stable transition metal orderings, or if there are preferred compositions and orderings that will emerge within low Co cathodes ultimately cause inhomogeneous electrochemical activity and cycling.
Density functional theory (DFT) alongside thermodynamic modeling has been successfully used as a method for understanding and predicting the stable phases of many binary alloys 22,23 and battery materials. [24][25][26] Density functional theory calculations train a computationally efficient model and with the use of statistical simulations using the model Hamiltonian, the thermodynamically stable phases with respect to composition can be predicted.
One of the largest factors in the accuracy of density functional theory calculations is the choice of the exchange correlation functional. Different exchange correlation functionals are shown to have drastically different levels of accuracy for different material classes with certain functionals performing better for certain applications. 27,28 Recently, a class of functionals known as Bayesian Error Estimation Functionals have been developed to incorporate a collection of exchange correlation functionals and therefore provide an estimate of uncertainty related to the calculation of exchange correlation energetics. Previous work has applied this uncertainty estimation to the prediction of magnetic ground states, bulk properties, and activity of catalysts. [29][30][31][32] Given the complexity of the problem, it is necessary to fully understand the uncertainty from systematic errors in the DFT training data as well as random errors from the statistical nature of Monte Carlo simulations. We focus most heavily on understanding the uncertainty derived from density functional theory. Reduced order models such as cluster expansion can be fit to high accuracy given enough terms, and the statistical error of Monte Carlo can be controlled and understood given a sufficient number of sweeps in your simulation as well as repeated simulations. The systematic error in density functional theory, however, is impossible to control as there is no way to find an exchange correlation functional for an arbitrary system with arbitrary accuracy. We therefore attempt to understand the sensitivity of our results to variation of the exchange correlation functional within a systematically chosen subspace of exchange-correlation functional space. Additionally, we attempt to show that without considering uncertainty propagation, the goal of identifying promising new low-Co phases could not be done with trustable results. We do however find that for the majority of the computational phase space, the energetic difference between the ordered and disorder phase is much larger than the the uncertainty of our calculations. Ultimately, within this work, we present an extensive computational search of the ternary phase diagram of NMC cathodes to understand the compositional regions of order versus disorder, as well as the most favorable compositions in the ordered region.

Computational Details
In this work, we use a total of 118 density functional theory calculations for the formation energy of various NMC compositions, focusing on high Ni-content structures, to train a reduced order model in order to more efficiently evaluate the formation energy with only a minor decrease in accuracy compared to input. The reduced order model is then solved in a high fidelity search of the composition space using Monte Carlo simulations that incorporate finite temperature effects. From these simulations, the change in Gibbs free energy at 298 K from pure layered metal oxide end member states to a NMC compositions at a given composition is fed into a convex hull analysis where we are able to predict all thermodynamically stable NMC phases. Here we describe briefly each step of this process providing the full details in the Supplementary Materials.
For all of the 118 density functional calculations used as training data, the Bayesian Error Estimation Functional with van der Waals (BEEF-vdW) 33 was used to treat the exchangecorrelation energy at the level of the generalized gradient approximation. Previous work has shown that systematic error related to oxide materials can be reduced if the prediction of energy differences is made with the proper reference state without the use of Hubbard U, even in the case of transition metal oxides where a U may commonly be used. 34,35 This was done by referencing the oxygen energy in the oxide to water rather than O 2 in order to match the oxidation state. This idea was further extended to show that the error in DFT predicted formation enthalpies of two similar reactions are correlated. That is the difference in formation enthalpies is constant and independent of exchange correlation function allowing the difference in these formation enthalpies to be compared accurately without the use of a Hubbard U correction. 36 Again this is achieved by comparing two similar chemical environments to one another. Within this work, we therefore chose to predict and compare the formation energy of the mixed NMC phases with respect to the pure layered oxide end members and do not include the Hubbard correction within our calculations. As this reference compares oxygens in nearly identical chemical environments, we expect the errors to be reduced drastically versus the prediction of the formation energy from elements.
The BEEF-vdW functional has the ability to estimate the error of density functional theory with respect to experiment through the generation of an ensemble of energetic predictions and was chosen for this reason. The ensemble of energies is generated by feeding the electron density of a fully self consistent calculation of the empirically fit functional from BEEF-vdW to an ensemble of exchange correlation functionals. The spread of the energetic predictions has been pretuned to recreate the error of the DFT calculation with respect to the experimental data it was trained on. In this way, BEEF-vdW links the precision of the ensemble of predictions, to the accuracy of the prediction of the self consistent calculations. With this empirical error estimation, we can provide a quantification of uncertainty within each DFT calculation and propagate that uncertainty to the model parameters that are trained with respect to the DFT input. Beyond the result that the spread of energies estimates the uncertainty of the calculation, the ensemble of functionals provides a way to probe exchangecorrelation space to understand the confidence of a prediction at the level of the generalized gradient approximations. This aspect can be used to ask if another exchange-correlation functional was used, would the prediction qualitatively change. The spread of the BEEF functional estimates has been been shown to bound the prediction of other general gradient approximation functionals for mechanical properties, magnetic ground states including, and reaction enthalpies for hydrocarbons. 29,30,36,37 As there is no way to ensure accuracy or even the estimation of accuracy, we argue this aspect of uncertainty quantification strengthens the interpretability of DFT prediction and predictions made by models trained on DFT.
Once the DFT data is generated, a reduced order model for the formation enthalpy with respect to the layered transition metal oxide end-members is chosen. This reference is chosen as it has been shown that for Ni, Mn, and Co can all be synthesized in the layered structure. 1,38,39 The change in volume of the lattice from pure states to the mixed state simulated by DFT is on the order of 1Å 3 and therefore the enthalpy is assumed to be the same as the energy predicted from DFT. The goal of this reduced order model is to recreate the change in enthalpy with respect to the end members of homogeneous layered lithium metal oxides, rather than the enthalpy from constituent elements. Therefore, our model focuses terms for the species that occupies a lattice site, and N-body interactions over various length scales as well as two body spin terms over various length scales given by: It should be noted that for the DFT calculation, the Co atom converged to the a very small spin. The inclusion of Co spin interaction versus ignoring Co spin interactions was explored. In order to do this as well as select which of the interactions considered gave the best reduced order model, four model selection techniques were implemented. The goal of this model selection was to determine the best model that balances predictive power with model complexity to prevent over-fitting. Figure 2 shows the results of these tests as a function of number of terms in the model. were also used. In each of the formulas, k is the number of parameters in the model,L is the maximized likelihood function, and N is the number of training points. We assume a Gaussian likelihood function for distribution of residual errors, x i , and therefore the maximum likelihood estimators for the mean, µ, and variance, σ 2 , are the mean absolute error and the sample standard deviation given by For each of these model selection techniques, the model with the lowest value of the metric is predicted to be the most ideal model. In our case, the same 42 parameter model which includes occupation, nearest neighbor spin not including Co spin, 4-body in plane, and 2-body out of plane next nearest neighbor interactions, was selected by leave one out cross validation, AICc, and BIC. The AIC selected a model that was slightly more complex than the 42 parameter model used which is a known bias of the AIC for smaller number of training points. 44 The final model gave a root mean squared error from leave one out cross validation of 5.86 meV/atom.
Metropolis Monte Carlo simulations were then performed using all supercells ranging from 5x5 to 10x10 to sample a large collection of possible compositions and ordering to ultimately find the lowest energy states. Each of these supercells also contained 3 layers in the z direction to reflect the A-B-C layer stacking structure of these layered oxide materials. In order to aid in converging the lowest energy ground state, a simulated annealing scheme starting at 1500K and stepping down to 0K at 150K steps was performed. At each temperature in this scheme, the simulation was allowed to equilibrate a total number of steps equal to 100 times the number of lattice sites in the simulation box. Once the annealing is done, the simulation is then equilibrated at room temperature and then sampled for another total number of steps equal to 100 times the number of lattice sites in the simulation box. The average energy is then reported for that composition and for the purposes of understanding the ordered solid solution phases, the configurational entropy is ignored. This is justified as once the simulation was equilibrated at room temperature, there were little to no accepted transition metal swaps during any of the thermodynamic samplings and thus we treat the ordered phases as fully ordered with respect to the placement of transition metals. There was however significant spin configurational entropy as the simulations were performed above the magnetic phase transition temperature of NMC. 46 When the entropy was computed via thermodynamic integration of the heat capacity, it was within 10% of the theoretical full spin entropy. The inclusion of full spin configurational entropy however, did not change the predictions for the ordered phases and therefore entropy was not explicitly computed for the ordered phases. The possibility of a disordered arrangement of transition metals was also considered. To simulate the disordered phase, a simulation was performed where all swaps were accepted and the energy of that configuration was sampled. This was done for 5000 steps and the average energy was reported. For the disordered phase, the full configurational entropy was included in the Gibbs energy and the results are compared to the ordered phase in Figure 3. The blue region shows the range of compositions predicted to be a disordered solid solution at that temperature. It is seen, however, that most compositions at room temperature favor an ordered arrangement of transition metals, with the region of disorder growing as temperatures reach that of synthesis.

Results
To further understand the energetic differences between the ordered and disordered phases Within this range the ordered and disordered hulls were predicted for comparison. At approximately x=0.2 and above the energy difference between order and disorder is well above the errors associated with the calculations and therefore we predict with high confidence the materials with lower amounts of Ni to be ordered. In the high-Ni region, however, the energy difference is much smaller and could arguably be within error. As most of the phases a room temperature were predicted to be ordered, we then perform a hull analysis in an attempt to find the lowest-lying ordered compositional states and determine what compositional states will appear as pure states within a battery at room temperature. To do this, the points were fed into a 3D convex hull algorithm and all points on the hull with formation energy less than or equal to 0 eV were plotted as points on a ternary phase diagram. Any composition that does not appear as a vertex of the convex hull is thus predicted to exist as a phase separated mixture of neighboring compositions that do appear as a vertex. Additionally, we simulate the stability of these phases at room temperature, the expected operational temperature, rather than elevated synthesis temperatures as states that are stable at high temperature but metastable at room temperature will likely end up in the stable state as the cathode is cycled accompanied by volume change and capacity loss. The effect of higher temperature with regard to ordered and disorder is explore in Figure 3 for a collection of elevated temperatures typical of synthesis. The effect of higher temperature increases the entropic term and therefore effectively weaken the relative strength of the short range ordering. This leads to a larger region of compositions that exist as a disordered solid solutions. As can be seen in Figure 5 Additionally, as seen in Figure 4, we have the largest uncertainty of whether the transition metals will be ordered or disordered at the high-Ni content region of interest. We therefore repeat our whole procedure of fitting the reduced order model, performing Monte Carlo at 298K to simulate thousands of compositions for the ordered and disordered phase, and feeding the results to a convex hull analysis for a set of 2000 ensembles generated from BEEF-vdW. We then count the number of times each specific composition appears on the convex hull as an ordered phase or as a disordered phase over all of the 2000 convex hulls ultimately generated. We then plot this normalized predominance in Figure 7. This results in a collection of predicted phases near the experimentally used phases around which there is a smearing of uncertainty. We also see that with respect to exchange correlation functional, there is no uncertainty that the majority of the phases are ordered.

Analysis
The model used within this work assumes a fixed layered structure and hence, we limit our discussion to phases below a certain threshold Mn content. Beyond this level, the material will exhibit a layered-to-spinel structural transformation that is seen in the LiMnO 2 end member. Given recent work on the structural phase diagram of the Li-Ni-Mn-Co oxide pseudo-quaternary system, we set this cutoff to be around 40% Mn content at zero Ni content and approximately 50% at zero Co content. 58 We also restrict our analysis to predicted phases that have a low Co content in hops of better satisfying the materials requirement surrounding the scarcity of Co. We show in Table 1 Table   S1 in the Supplementary Materials.
As mentioned before, the 811 phase is a promising candidate for a high Ni-content cathode yet we do not recover this with high prediction confidence (c=0.011). We do however The prediction of the exact compositions of NMC phases will allow more accurate predictions of the properties and degradation of these cathode materials. Further work using the specifically predicted phases here may allow better understanding of the tradeoffs of low Co cathodes and the limitations of the design space. While we do not present a full analysis of these phases, the work here provides a starting point to explore and understand the performance of these new phases and the possibility to achieve high performance ultra low-Co cathodes.
One of the easiest predictions that can be made once the Gibbs Free energy is known is the average voltage of the predicted NMC cathode with respect to Li/Li + given by the Nernst Equation.
V = −∆G f F  The change in Gibbs energy is given by If we assume that the delithiated phase has the same Gibbs energy as the delithiated end members as is discussed in the Supplementary Materials, we can calculate the average voltage A plot of the predicted average voltage over the entire phase can be seen in Figure 8 and a table of the predicted voltage of common NMC phases is presented in Table 2. The prediction is compared to the predictions from AutoLion-1D, 60 a well benchmarked with respect to experiment model which allows for the simulation of an extremely slow discharge to better match the perfectly ideal average voltage predicted by this work. The results of this reproduces the benchmarked data to within 0.1 V for all NMC phases tested.   We see in Table 3 the predicted energies of these defects are rather large and predict extremely small quantities of these defects at both room temperature and even synthesis temperature as compared to some experimentally measured occurrences of this defect. We therefore assume that while this kind of defect is present in experimentally made cathodes, the largest driving force for the creation of the defects is unlikely to be thermodynamic but rather process related and the choice to ignore this effect within the main work is justified.   were used to increase the model accuracy in the high-Ni region of the phase space. One extra calculation of the 532 phase in a 5x2 supercell from Wu el al. 4 was also include to assess the proposed cation ordering of a NMC532 phase. The initial spin state of the DFT training data was randomly chosen to be positive or negative to create a random sampling of spin interaction and in most cases the final energy converged to the same spin state as was initialized. The final spin state, however, was always used in generating the cluster terms. extracted. This process was then repeated by fixing the a and b constants and varying the c constants then fitting to find the ideal constant in the way as before. The resulting structure at the optimized lattice parameters was then relaxed to a maximum force of 0.01eV/Å.

Reduced order Model
It should be noted that the inclusion of all n-body terms in the model would lead to a rank deficient system of equation when trying to fit the interaction coefficients through a system of linear equations. This is because there is an external constraints of all the lattice sites being filled that creates a linear dependence on the number of cations with the number of nearest neighbor interactions, next nearest neighbor interaction, etc. Because of this, we choose to artificially set the all homogenous n-body interactions to be zero. What is left in the model is the relative favorability of a given interaction compared to all other interactions.
The model was fit using multivariate normal regression where the design matrix of number of interactions and the formation energies were normalized to one stoichiometric unit of LiMO 2 .
This prevents overweighting the larger supercells. The magnitude of the spins of each cation is assumed to be constant for all composition of Ni, Mn, and Co with the same level of lithiation. This assumption allows for the reduced order model to be solved more easily since allowing the magnetic moments of each atom to vary continuously and unbounded during Monte Carlo simulations would cause the spin interaction term to be unbounded.
Once the final model was chosen, the convergence of the model as a function of number of training points was tested with leave one out cross validation as seen in Figure S1. The root mean square error levels of at around 100 training points and appears constant at the full 120 point data set.

Monte Carlo Simulations
Once the model is trained, the enthalpy of formation at finite temperature is simulated using a Metropolis Monte Carlo algorithm. For each composition of transition metals, a random initial cation and spin ordering is assumed. Then randomly the spin of a single cation or the position of two cations is swapped. After each swap, the condition for acceptance is given by estimating the probability of the swap happening spontaneously.
Where ∆E is the energy change of the swap. If the energy change is negative then the swap is always kept. If the energy change is positive, it is kept if the probability is less than than a randomly generated number between 0 and 1. This algorithm allows for the thermal population of states and given a sufficiently long simulation, will accurately predict the probability of thermal population.
Each simulation is thermalized via simulated annealing starting at 1500K and then decreasing by steps of 150K to 0K with each step simulated over 100 times the number of lattice sites. The enthalpy is then calculated to be the average enthalpy of all of the sweeps.
This procedure can be seen for an example case in Figure S2. Energy per Formula Unit Figure S2: A plot of the energy per formula unit as a function of step number during the simulated annealing step followed but the 100 steps per number of lattice site thermalization step at 298K. The simulation is for a 9x9 supercell to demonstrate the convergence of a larger supercells.