Thermodynamics of multi-sublattice battery active materials: from an extended regular solution theory to a phase-eld model of LiMnFe<1-y>PO

Phase separation during the lithiation of redox-active materials is a critical factor affecting battery performance, including energy density, charging rates, and cycle life. Accurate physical descriptions of these materials are necessary for understanding underlying lithiation mechanisms, performance limitations, and optimizing energy storage devices. This work presents an extended regular solution model that captures mutual interactions between sublattices of multi-sublattice battery materials, typically synthesized by metal substitution. We apply the model to phospho-olivine materials and demonstrate its quantitative accuracy in predicting the composition-dependent redox shift of the plateaus of LiMn y Fe 1-y PO 4 (LFMP), LiCo y Fe 1-y PO 4 (LFCP), LiCo x Mn y Fe 1-x-y PO 4 (LFMCP), as well as their phase separation behavior. Furthermore, we develop a phase-field model of LFMP that consistently matches experimental data and identifies LiMn 0.4 Fe 0.6 PO 4 as a superior composition that favors a solid solution phase transition, making it ideal for high-power applications.


Introduction
Li-ion batteries are fundamental to the upcoming transition towards sustainable energy production, electric mobility, and energy storage 1 . Although the early storage requirements were satisfied by active materials such as graphite and LiCoO2 2,3 , higher energy densities, sustainability, cheaper elements, and improved safety require developing more sophisticated battery active materials. Li-ion battery electrode materials also have other emerging applications 4 , such as electrochromic displays 5 , ion-tunable electrocatalysis 6 , resistive switching memory [7][8][9] , water desalination and purification 10 , and lithium extraction from brines 11,12 . In all of these applications, the design space for electrode materials with various desired properties has hardly been explored.
Blending or modifying existing electrode materials is a promising method to improve properties, which is gaining attention, albeit with limited theoretical guidance. While the anode materials are moving towards silicon 13 , silicon/graphite composites, or Li-metal 14 , cathode development is running behind, with most advancements focusing on substituting cobalt in layered oxide materials with Ni, Mn, or Al developing LiNixMnyCo1-x-yO2 (NMC) and LiNixAlyCo1-x-yO2 (NCA) cathodes 15,16 . These approaches show the advantages of modifying the composition of an existing cathode with well-established lithiation mechanisms to reduce its cost and environmental impact and to improve energy density and cycle life. Applying the same approach to LiFePO4 (LFP), a phosphor-olivine material introduced by Goodenough and co-workers in 1997, 17 which has advantages over the layered oxides in lower cost and toxicity with greater stability and recyclability 18 , various partial or complete substitutions of Fe with Mn, Co, and Ni have been attempted. Higher redox potential and similar specific capacities can be obtained, improving the overall energy density [19][20][21] while sustaining decent diffusivity and cycle life. Currently, LiMnyFe1-yPO4 (LFMP) [21][22][23] exhibits the most promising characteristics and is rapidly being incorporated into commercial batteries. Therefore, it is crucial to gain a deep understanding of the basic physics of LFMP through modeling. Other materials in the same family, such as LiCoyFe1-yPO4 (LFCP) [24][25][26][27][28][29] and LiCoxMnyFe1-x-yPO4 (LFMCP) 19,20,30 , also display intriguing properties and merit further investigation as well.
First-principles calculations struggle to provide a complete picture of the underlying mechanisms, in part due to the heavy impact of the practical choice of the pseudopotentials on the predicted redox potential 31,32 and partly because the use of Monte Carlo simulations aided by cluster expansion 33 prevent the understanding of the behavior of the material in a realistic battery system at finite temperature.
In order to develop an accurate thermodynamic description, it was essential to model the behavior of individual particles using phase-field methods, which generalize the Cahn-Hilliard formalism [55][56][57][58][59][60][61] for driven electrochemical systems 34,48,52 . This approach has led to realistic models of diffusion and reaction models for materials such as graphite 62,63 , anatase TiO2 45 , LTO 46 , LCO 9 and LFP [64][65][66][67][68] , showing excellent agreement with experiments, guiding researchers to properly understand the reasons for various peculiar behaviors occurring in phase-separating materials and helping companies in the optimization of these kinds of batteries. Recently, phase-field modeling of LFP has succeeded in reproducing a vast dataset of operando x-ray images of nanoparticles cycling at different rates pixel by pixel 69 , while learning the two-phase free-energy landscape, the reaction kinetics of coupled ion-electron transfer 70 , and the heterogeneity of surface reactivity, correlated with variations in carbon coating thickness. An open-source code, MPET (Multiphase Porous Electrode Theory) 48 has been developed to facilitate the implementation of these models for specific cells and control algorithms 71 . With its modular design, users can quickly incorporate phase-field models of the studied material within a porous electrode theory framework, providing insights into both the individual particle and the collective system responses.
In this study, we applied a thermodynamic-based approach to investigate the impact of composition on the performance of phospho-olivine materials. We extended the regular-solution theory 61,72 , originally applied to single-lattice LFP 34,73 , to consider the presence of multiple sublattices for the intercalated species, which have distinct properties and redox potentials. Our theory reveals how interactions between sublattices explain the composition dependence on redox potential and phase transition behavior. By applying the theory to a phasefield simulation of an LFMP half-cell, we gain fundamental insights into the optimal transition metal ratio, which demonstrate the possibility of using mean-field phase-field models to design active materials.

Theory
The mathematical modeling of a closed thermodynamic system starts by defining the Helmholtz free energy , given by = − 61,72 . To determine the properties of a solid solution, we need the temperature , a function for the entropy , linked to all possible configurations the system can have, and the internal energy , representing interactions between the particles of the system. The regular solution model, which is equivalent to a mean-field lattice-gas model with pair interactions, offers an elegant and straightforward way to describe solid solutions, and it was implemented successfully in various battery active materials 34,40,[44][45][46]48,73,74 . However, its application is mainly limited to materials where the intercalated species encounter one lattice type, such as LFP or LTO 46 , or two non-interacting sublattices, such as TiO2 45 . Phase-field models of staging phase transitions in graphite have been developed with multiple, periodic interacting crystal layers, 40 but the parameters must be fitted to experimental data to describe the complex phase diagram of the material 62,63,75 .
To summarize this mean-field theory, we can start considering an active material containing lattice sites that can host intercalating species (e.g., lithium) whose relative concentration is defined as ̃= / . The ε ), which represents the mixing enthalpy, determining whether the material will favor phase separation during (de)intercalation # .
One way to obtain the parameters introduced above is by performing first-principles calculations, attaining the energies of the structure at different fractions of intercalation, and subsequently determining the mixing energy and the redox potential 76 . Moreover, it is possible to measure the chemical potential experimentally during close-to-equilibrium (de)intercalation; fitting the voltage hysteresis gap will then capture the difference between the local minima and maxima of the Ω dependent chemical potential 44,77 .

Multi-sublattice model
Aiming to describe a general multi-sublattice material in which the lattice sites are mixed uniformly, the usual regular solution model must be modified to account for the various interactions that the intercalated species can have in the material.
We consider a lattice composed of intercalating sites, divided into sublattices, of which = belong to the sublattice . Assuming that the lattice sites of different types are fixed in space after the synthesis of the entropy for the sublattice , having an occupation ̃= / , is calculated using the solid solution approach, resulting in The internal energy of the general system with different sublattices is computed using a mean-field pair interaction approach so that the effect of the sublattices interactions arises naturally. Using the same logic behind the regular solution model, we can start considering a scheme represented in Figure  1, in which, for clarity, only two types of sites are considered. It shows how an intercalated particle in the occupied lattice site will interact with other occupied sites of type , the vacancies of the sites , but also with the vacancies and the occupied sites of type , leading to 9 different interaction energies for a 2-sublattice material.
The evaluation of the interaction energies between different sites ε , in such complex systems, is challenging and would involve atomistic quantum mechanical computations. As a first approximation, we calculate them as an average between the same-site interactions ε = (ε + ε )/2.
Assuming that the site-site interactions ε are not influenced by the ratio between the compounds , the total internal energy will so be obtained considering that an atom or a vacancy in the sublattice , having close neighbors, will interact with = neighbors with an energy , or depending on their occupation state, but also with = neighbors with an energy , or . Applying this concept and considering all the possible combinations, the internal energy is Where Ω can be calculated as and it will coincide with Ω = (ε − 1 2 ε − 1 2 ε ), the value for the single lattice structure where only one kind of interaction is present.
Equation 2 links the properties of the original materials, which are summarized in Ω , Ω , … , Ω , to the properties of the mixed compound. We can thus conclude that this one equation describes the system in all its possible compositions, becoming a powerful tool for alloy engineering.
Distinguishing now between the absolute energy dependence on concentration (first two terms is eq. 2) and the enthalpy of mixing (last term in eq. 2), the free energy of the system can be rewritten as Where Θ is the standard chemical potential for the sublattice , which is correlated with the standard half-cell potential of the reacting sublattice Θ so that in the case of Li-ion batteries Θ = − Θ / .
Once the regular solution theory has been reformulated for a multi-sublattice active material, it is crucial to analyze the analytical solution to gain insight into the impact of various alloying elements on the material's behavior before obtaining the free energy functional needed to construct a comprehensive phase-field model.
To predict the behavior of the system is necessary to build an dimensional energy space and follow the concentration path that minimizes the energy extracted to transform the system from a completely deintercalated state, where ̃1,̃2, …̃= 0, to a fully occupied system where ̃1,̃2, …̃= 1. In this way, it is possible to numerically obtain a solution for (), and from it, a voltage curve for a homogenous single particle system can be obtained ().
If the difference in standard chemical potential between the lattice sites is significant compared to Ω we will observe a series of redox plateaus in the voltage curve. Conversely, a more complicated energy path will be followed (see supplementary information, section 2). Limiting ourselves to the first case, we can analytically calculate the chemical potential, and so the voltage curve of the various plateaus as This simple formulation allows us to analytically capture how the system behaves depending on the compositions , , … , and the known factors Ω , Ω , … , Ω .
In fact, considering an intercalation process, we can deduce that in case 1 Θ < 2,3,… Θ the intercalated species will initially sit in the lattice sites here defined as "1" keeping ̃2 ,3,… = 0, and then, once ̃1 ~ 1 the second lattice sites will react, and so on. Therefore, the effective standard chemical potentials measured will be From equation 7, we can conclude that, due to the enthalpic contributions of the surrounding sublattices, the effective chemical potential, and so the measured redox potential Θ , will be decreased by a value corresponding to the sum of the pair mixing energy coefficients, weighted by their corresponding stoichiometry. For the second plateau, we instead expect that since 1 ~ 1 the redox potential will increase by a factor 1 Ω 12 and be reduced by a factor ∑ Ω 2

=3
. Moreover, it is worth noticing how the effective enthalpic interaction Ω = Ω is now a function of the stoichiometry of the compound. In the case of phase-separating materials, this will impact the voltage hysteresis gap and the overall phase separation behavior with respect to the pure original lattice. Knowing that if Ω < 2 no phase separation will occur during the plateau of the specie " ", we can directly calculate the compositions that assure a solid-solution transition < 2 /Ω depending directly on temperature and the mixing energy of the original compound.

Application to Phospho-olivine cathodes
The available literature and the need for a reliable electrochemical model make phospho-olivine cathodes an excellent candidate for a critical verification of the developed theory by testing its implications on the redox potential shift and the order of the phase transition. We focus on the most promising iron substitutions in LiFePO4 (LFP) by Mn and Co. Cyclic voltammetry and voltage curve data from the literature can provide the necessary information to implement in the model. For Ω and Θ , we can rely on previous regular solution models for LFP 41 were Ω = 4.63 and Θ = −3.422 . However, for LiMnPO4 (LMP) and LiCoPO4 (LCP), no previously developed models are available, so we must obtain Ω = 7. 44 and Θ = −4.09 from fitting the voltage curve from Tasaduk et al. 20 and the cyclic voltammetry of Kobayashi et al. 78 , respectively. Obtaining values for LiCoPO4 presents particular challenges due to its instability with the current electrolyte 24 , limiting the available data. Nevertheless, we note that the peak separation of the cyclic voltammetry in the work of Jalkanen et al. 79 is close to that of LFP. Since the peak separation is proportional to the voltage gap due to phase separation, we can assume that Ω ≈ Ω and that Θ = −4. 78 can be determined from the peak midpoint. Fig. 2 Comparison between the calculated (lines) and the measured (empty dots) redox potential at various Mn and Co substitutions ( ). The green line is the redox potential shift of the Co plateau in LiCoyFe1-yPO4, while the dark orange is its counterpart in the Fe plateau. The blue line is the redox potential shift in LiMnyFe1-yPO4, and the light orange one is the redox shift in the corresponding Fe plateau.
Starting by analyzing the shift in redox potential, we can compare the midpoints of the cyclic voltammetry results at different Fe substitutions from references [57] and [59] with the redox shift obtained from the theory. The expected dependence of the redox potential on the Mn content for the plateaus in LiMnyFe1-yPO4 will so be Making the same calculations for the LFCP, we can see in figure 2 how the theory predicts the experimental values without any fitting parameters. A strong indication of the validity of the theory is in the quantitative prediction of the shift in the Fe plateau Δ Θ = ( − 1)Ω which clearly does not depend on the original redox potential but only on the value of the average enthalpic interaction Ω , which is different for the case of Co and Mn substitution. For the LiCoPO4 case, we could have been more precise by considering the double redox plateau present, seemingly linked to a staging behavior 28,29 . At present, we chose to neglect this effect, assuming it would pose minimal effect on the pair interaction energies, therefore, the chemical potential of LCP can be approximately modeled in the same fashion as that of LFP and LMP.
These results establish a strong foundation for the mean-field theory and offer a clear explanation for how the redox shift can be attributed to the interaction between surrounding lattice sites, in agreement with the Monte Carlo calculations of Malik et al. 33 in which the shift in redox potential was attributed to the effect of the pair interaction energies. It is remarkable that the straightforward assumption of averaging the two interaction energies, as demonstrated in eq. 3, continues to hold true, even when the phenomenon is rooted in atomistic behavior. The predicted equilibrium phase transition behavior aligns with the available data in the literature. The experiments of Jalkanen 79 and Kobayashi 78 show the CV peak separation of the corresponding plateau to be linearly dependent on the composition of the olivine material. Further, the works of Ravnsbaek 80-82 and Strobridge 29 analyze the operando XRD profiles of different compositions of LFMP and LFCP, exposing the absence of phase separation for the corresponding sublattice if the composition is below the one calculated from the model. Therefore, the mathematical theory and the associated phase diagram can become tools for the practical engineering of these alloys. For example, they enable the selection of a composition range in which the material behaves as a complete solid solution (red triangle), at least from an equilibrium thermodynamics perspective. We expect a wider zone (light red) in which the real system can behave as a solid solution due to the stabilizing effects of coherency strain 41 and auto-inhibitory intercalation reactions 52 and the relation between the particle dimensions and the phase separation front 29,83 . Within the solid solution region, it is also possible to select the composition that minimizes the Co content (orange curve for LFCMP in figure 3). The complete solid-solution behavior is confirmed by experiments on both LiMn1/3Co1/3Fe1/3PO4 30 , and LiMn0.3Co0.2Fe0.5PO4 84 in which the systems show a monotonically decreasing voltage curve.
These conclusions come directly from the analytic application of the extended regular solution theory and the consistent calculations of Ω , Ω and Ω without the need for ab-initio simulations. The mean-field model also helps in the phenomenological description of the system. The dilution of a sublattice, and its consequent reduction in first neighbors, weakens the attractive interaction between the intercalated atoms allowing the entropic contribution to take over, leading to a solid solution. Since this effect is subtle, a temperature change will also lead to different behavior (see supplementary information, section 2). This conclusion differs from the one reported by Malik et al. 33 in which the disappearance of the phase separation for certain compositions was attributed to the reduction in the Li composition difference between the initial and the final state.
Finally, the combination of the predictions for the redox potential shift and the phase transition naturally leads to the possibility of calculating the voltage curves for every composition, including LFMP, LFCP, and LFMCP, in good accordance with the previously cited experimental results.

Phase field modeling
Having demonstrated a correspondence between experimental data and the analytical solution of the model, it is now interesting to investigate other factors that may come into play when the system is out of equilibrium. This can be accomplished by employing a complete phase-field simulation, which takes into account factors such as coherency strain and gradient penalties. Beginning with a simulation of a single particle, the behavior will be observed during both charging and discharging, with the aim of gaining insight into the collective dynamics that arise at various compositions.
Given the wealth of available experimental data in the literature and the potential for commercial applications 85 , we have narrowed our focus to LFMP simulations. We implement our model in the open-source code MPET 48 , freely available in its GitHub repository. The complete set of equations and the parameters of the phase-field simulations can be found in section 1 of the supplementary information, alongside a comparison of the simulated and measured voltage curves at equilibrium.

Simulation results
Single particle simulations Fig. 4 Evolution of the normalized Li concentration inside a 150 nm particle upon lithiation. The blue lines correspond to the Mn plateau, the orange lines to the Fe plateau. The labels show the average composition of the particle.
The single particle simulations in equilibrium (C/1000) (fig 4) show how modifying the enthalpic contribution at different compositions affects the system. Given the fraction of Mn in LiMnyFe1-yPO4, for the cases of = 0.6 and = 0.8, we obtain Ω < 2 at ambient temperature so that the particle, as expected from the analytical calculations, behaves as a solid solution during the (de)intercalation of the Fe plateau and phase separate when in the Mn plateau. In contrast, for intermediate values of within the range of 0.2 to 0.4, although the effective interaction energies Ω and Ω exceed the critical threshold of 2 for phase separation, the coherency strain provides a stabilizing effect that leads to the transformation of the particle into a solid solution. The insertion direction in our one-dimensional model is the one in which the coherency strain is minimum, which coincides with the preferential direction for phase separation 41 . This implies that the observed solid-solution behavior will remain consistent when considering a three-dimensional particle.
This claim requires further experimental verification, keeping in mind that the composition where we experience suppression of the phase separation may slightly differ from the one observed in simulations due to the documented sensitivity of the calculated coherency strain values on ab-initio simulation parameters 86,87 . However, the single-particle simulations rationalize how the composition LiMn0.4Fe0.6PO4 with its probable solidsolution behavior can limit the problems due to the measured low Li diffusivity 88 .

Porous electrode simulations
We applied our single-particle model in a porous electrode system to explore the collective many-particle behavior in a real cathode. In particular, it is interesting to simulate the effect of the composition on the possible suppression of phase separation, already known for LFP , and on the lithium distribution along the depth.
We distinguish the phenomena by performing simulations at two different Mn content: the first consists of simulating the charge and discharge process of a thin electrode at C/10 to assess the collective dynamic, and the second involves a 1C discharge of a thicker electrode to focus on the transport limitations. All the simulations are done on an ensemble of 400 particles with lognormal size distribution.
The computed evolution of the concentration during the C/10 cycle was collected in a probability distribution and converted to normalized volume expansion (see supplementary information, section 1) to better compare the result with the work of Ravnsbaek 81 et al. The results, shown in figure 5, not only strongly agree with the work of Ravnsbaek et al. 81 , but they also offer a thermodynamically consistent explanation of them. Both in experiments and simulations, a bimodal volume distribution is present in the plateaus where the theory predicts phase separation. Moreover, it is essential to state that, considering the cases of LiMn0.2Fe0.8PO4 and LiMn0.4Fe0.6PO4 where, as discussed in the previous paragraph, we predict a solid solution transition for the single particle, the collective dynamics are dominated by the non-monotonic shape of the chemical potential and its subsequent concentration-dependent exchange current density 52 . Since we are close to equilibrium conditions (C/10), we can conclude that the origin of this bimodal distribution is established by the inter-particle separation (mosaic lithiation) in which the smaller particles are more lithiated than, the bigger ones, as also observed for LFP 97 .
Focusing on the asymmetry between charge and discharge, for LiMn0.4Fe0.6PO4, both in simulations and experiments, the phase separation of the Fe plateau is present only during charging. From this observation, Ravnsbaek et al. 80 suggested that the intrinsic order of the phase transition in LiMn0.4Fe0.6PO4 depends on the direction of the transition and attributed the reason to coherency strain effects. Our consistent thermodynamic model enables us to reinterpret these conclusions offering a physical explanation for the observed experimental behavior. Our simulations indicate the significance of the collective auto-inhibitory and auto-catalytic behavior upon lithiation and delithiation, respectively 52 . Due to the asymmetric concentration dependence of the exchange current density upon delithiation, an enhanced particle-by-particle reaction is observed, while during lithiation, the inter-particle separation is suppressed, as previously observed in LFP 97 and NMC 54 porous electrodes. While this phenomenon is only observed when cycling LFP at high rates, it is instead already present at C/10 in the Fe plateau of LFMP due to the low Ω . The only observed mismatch in figure 5 occurs in discharge case of LiMn0.2Fe0.8PO4 and can be attributed to various phenomena not included in the model such as a possible metastable phase for the Fe plateau 98,99 , effect of particle size 83,100 or non-linear dependence of the volume on the Li concentration.
Finally, it is expected that already at 1C, in the case of a thin cathode, in the range 0.2 < < 0.4, the interparticle separation is completely suppressed during lithiation (see supplementary information, section 2), even if low intra-particle diffusivity may lead to different experimental results. To complete the picture, we studied the effect of the composition on the lithium depth profile during a 1C discharge of a commercial-like cathode. The simulated half-cell has a cathode thickness of 80 a porosity of 30%, the transport limitations of Li in the electrolyte are therefore not negligible anymore. Due to the absence of the nucleation barrier, solid solution materials tend to lithiate relatively uniformly along the depth, even in the case of strong transport limitation. On the other hand, phase-separating materials show steeper gradients in the Li concentration 101,102 . The mixed phase separationsolid solution behavior of LFMP makes it an exciting candidate to evaluate this effect since the different compositions affect the phase transition. As shown in figure  7, if LiMn0.8Fe0.2PO4 is used, the Mn plateau will show an inhomogeneous redox activity, typical of phaseseparating materials, while the last 20% of the discharge, corresponding to the Fe plateau, will have a uniform reaction along the depth. Similar considerations can be made on LiMn0.4Fe0.6PO4 in which a small but visible peak in redox activity can be seen in both plateaus due to the non-monotonicity of the chemical potential, which is still present even if the single particle is not phase separating.
We can so conclude that the Li distribution along the depth of the electrode at different states of charge is severely affected by the concentration of Mn. Taking LFP as a reference, it is clear that the change in the degree of phase separation helps guarantee uniform lithiation along the depth, and the case of LiMn0.4Fe0.6PO4 is the one that most favors homogeneity. This result significantly impacts the cycle life, thanks to the possibility of reducing current hot spots, and offers a new route for a composition-based optimization of a commercial electrode.

Conclusions
In this study, we expanded the regular solution theory to explain and predict the behavior of phosphor-olivine cathodes. The inclusion of multiple sublattices and their interactions provided an elegant explanation for the shift in redox potential and phase separation behavior. The mean-field theory formalization offers an intuitive understanding of these phenomena, which can enable research on new active materials. This approach can serve as a valuable alternative to computationally extensive ab-initio calculations, delivering clear insights starting from simple concepts instead. Our phenomenological description of the mathematical derivation demonstrates that the redox shift is due to the interactions of the non-reacting sublattice. Furthermore, we found that a redox plateau that previously showed phase separation can transform into a solid solution. This transformation occurs due to the reduced number of closest neighbors within the same sublattice, which lowers the effective interactions of the intercalated species.
The application of our model to well-studied materials such as LFMP, LFCP, and LFMCP and their possible compositions shows how quantitative and accurate this theory is, even if the examined system is considered complex. The subsequent application in a phase-field framework was able to reproduce and explain various experimental results whose interpretations were incomplete and lacked mathematical support. The firm conclusion about the absence of phase separation in low Mn content LFMP is still to be confirmed experimentally. However, the proposed mechanism to explain the operando XRD peak shift sheds light on the importance of considering multi particles behavior when experimenting with a collective system such as a halfcell. Finally, this model strongly indicates the optimal composition for a high-power cathode, showing how LiMn0.4Fe0.6PO4 may be an excellent candidate thanks to its solid solution behavior and low transport-induced inhomogeneity. To verify this claim, further experiments are necessary, and a 2-dimensional model, able to capture the known transport limitation in the particle, is also advised. The interplay between the concentrations of the two sublattices may play an essential role in explaining the out-of-equilibrium behavior opening the route for optimization of the (dis)charging procedure to exploit these effects 103 . We expect that the new theory may be applied to other popular active materials such as LiMn1.5Ni0.5O4 (LMNO) or the various compositions of NMC, explaining the effect of metal ratios on the performances. To do so, it will be necessary to consider the structural modifications occurring in the spinel or the layered structure that, at the moment, are not taken into account in the theory. It is finally hoped that overcoming these limitations may expand the domain of this theory in such a way that, as particle dimension, porosity, and thickness, also the composition of the materials can be included in the parameters to optimize when a battery is designed, improving cycle life and energy efficiency.

Code and Data availability
MPET is available as a Bitbucket repository: https://bitbucket.org/bazantgroup/mpet/src/master/ but at the moment, it does not contain the LFMP model in the master branch. The reader is thus invited to visit the GitHub repository of the code to explore the different branches: https://github.com/TRI-AMDD/mpet.git. The opensource nature of the code makes also possible to reproduce the described simulations using the parameters described in table 1 of the method section of the supplementary information.