Modeling the Solid Electrolyte Interphase: Machine Learning as a Game Changer?

The solid electrolyte interphase (SEI) is a complex passivation layer that forms in situ on many battery electrodes such as lithium‐intercalated graphite or lithium metal anodes. Its essential function is to prevent the electrolyte from continuous electrochemical degradation, while simultaneously allowing ions to pass through, thus constituting an electronically insulating, but ionically conducting material. Its properties crucially affect the overall performance and aging of a battery cell. Despite decades of intense research, understanding the SEI's precise formation mechanism, structure, composition, and evolution remains a conundrum. State‐of‐the‐art computational modeling techniques are powerful tools to gain additional insights, although confronted with a trade‐off between accuracy and accessible time‐ and length scales. In this review, it is discussed how recent advances in data‐driven models, especially the development of fast and accurate surrogate simulators and deep generative models, can work with physics‐based and physics‐informed approaches to enable the next generation of breakthroughs in this field. Machine learning‐enhanced multiscale models can provide new pathways to inverse the design of interphases with desired properties.


Introduction
The key to designing higher performance and safer lithium ion batteries (LIB) lies in establishing a clear understanding of how the solid electrolyte interphases (SEI) form and evolve over multiple

SEI Complexity and Modelling Challenges
The formation mechanism(s) and structure of the SEI [30] remain poorly understood, due to a) dependence on the specific electrode and electrolyte composition, b) the sheer number of species involved, c) the reaction network between them is vast and intractable. The critical properties and dynamics of the SEI reaction networks emerge on long timescales as crucial reactions are rare events with high kinetic barriers accessible only under applied potential. Additionally, the actual composition and structure depends on the salt and solvent used in the electrolyte and on the surface termination of the electrode. For example, just changing the electrolyte can create a radically different composition and microstructure of the SEI. [7,31] The ongoing reaction between the initial reaction products and the electrolyte produces both organic and inorganic components. The SEI changes upon long-term cycling, resulting in an SEI with a greater content of inorganic species. [32,33] The type of materials and molecules constituting the SEI, as well as the meso/macro scale structure evolve over time. Thus the holistic view of SEI evolution and its impact on battery lifetime, performance and safety needs to resolve a diversity of entangled phenomena spread across length and time scales in both the crystalline and soft phases and interphases; which constitutes a fundamental theoretical challenge.
A dual-layer structure of SEI is typically described with an inner solid inorganic layer and a porous polymeric outer layer. [6,34] The thickness and the formation mechanism of these layers are still debated, for example, it is not clear whether longterm SEI growth proceeds via transportation of negative charge to the SEI/electrolyte interface or by electrolyte diffusing toward the electrode/SEI interface. Alternatively, a mosaic structure [4,35] of the SEI has also been observed, where particles of inorganic crystalline phases are embedded in an amorphous organic matrix. Not only the micro-structure, but also the chemical compositions are also often disputed. [10,36] A core direction in recent research efforts has been aligning experimental and computational studies over many observables across time and length scales to provide the clarity of understanding that has so far eluded the research community. [37][38][39][40][41] From a computational point of view, modelling needs to be performed at multiple scales, ranging from electronic structure or electron transport properties of individual phases or boundaries [43] to mechanical degradation like crack formation [44] (see schematic in Figure 1). Achieving predictive accuracy at each of the length and time scales is a challenge as most purely physics-driven methods suffer from a trade-off between time/ length resolutions, sensitivity, time/length scale of observation, and computational cost. Approximations help us study longer time/length scales, but are often not useful as crucial phenomena happening at a shorter time/length scale are missed. A multi-scale approach [45] has sometimes been utilized by using high accuracy small-time/length scale simulations to parameterize longer scale simulations. Such directions suffer from non-scalability (as for small changes to the system, reparameterization is needed for each scale bridging) and thus can not be used for quickly searching through a large number of structural/dynamical/compositional possibilities. This limits the possibilities of high throughput computation-based design of better SEIs or inverse design of the SEI with desirable properties encompassing the lifetime of batteries.

Quantum Chemistry Methods
Ab initio or quantum chemistry (QC) methods-either based on wavefunction (WF) descriptions or density functional theory (DFT)-provide the full electronic structure of an atomic system within certain approximations. [46] Therefore, these methods are ideally suited to study electrochemical processes, which addressed by modeling techniques discussed in this review: (Å/ps) QM/ DFT methods have been extensively used in order to model SEI formation reactions in terms of the corresponding reaction rates, to study the impact of electrode surfaces on these rates, and to simulate bulk SEI phases. (nm/ns) Reactive MD simulations provide a systematic method to map QM-based reaction networks of solvent reduction and dimerization to the chemical composition and morphology of the resulting SEI structure. (nm/μs) Classical MD simulations can be used to investigate insertion of lithium ions from the electrolyte into the SEI, migration of lithium ions in the bulk SEI, and extraction of lithium ions from the SEI for a given SEI morphology and chemical composition. (m/s) Continuum models adopt the Butler-Volmer kinetics to study the effects of the formation of the SEI on the electrochemical performance of the cell. Image related to "Reactive MD" adapted with permission. [164] Copyright 2021, AIP Publishing. Image related to "Classical MD" adapted with permission. [42] Copyright 2014, American Chemical Society. encompass both electron transfer processes and consecutive (radical) reactions, and have consequently been widely used for modeling the nucleation of the SEI. Due to its favorable computational scaling as compared to WF methods that explicitly include electron correlation, DFT has become the workhorse in most of these applications. Typically, DFT is either employed in ab initio molecular dynamics (AIMD) [47] simulations, in which the system is propagated in time by integrating Newton's equations of motion on the potential energy surface evaluated at the DFT level, or in (static) calculations of molecules, molecule fragments and clusters, or crystal structures, providing access to reaction or transport energies and barriers. In this section, we briefly review these individual aspects in context of SEI formation and discuss their main limitations. For a more detailed review, we refer the reader to the recent review by Wang et al. [5] and an overview article by Li et al. [48] 3. 1

.1. Modeling of SEI Formation Reactions
Calculations of Redox Potentials: The SEI formation is initiated by reduction (i.e., electron transfer) processes from the anode to electrolyte molecules in the vicinity. However, typical battery electrolytes are composed of two or more different carbonate solvents, lithium salt(s) as well as additives fulfilling specific tasks. [8] To determine which of these molecular species are most prone to reduction, the estimation of their electrochemical stability-characterized through the decomposition potential in experiments-is an important aspect that has been studied via QC calculations.
Generally speaking, the oxidation or reduction stability expresses the tendency of the molecules and ions in an electrolyte to accept excess electrons (in case of reduction) or to give away electrons (for oxidation). In a most basic approach, the redox stabilities are estimated via the highest occupied molecular orbitals (HOMOs) and lowest unoccupied molecular orbitals (LUMOs) of the molecules, [49,50] or via the closely related ionization potentials (IPs) and electron affinities (EAs), [50] which additionally involve a second DFT calculation of the oxidized/reduced molecule. The latter accounts for the change in molecular relaxation upon oxidation/reduction and is intimately related to the reorganization energy within the framework of Marcus theory. [51,52] Although approximating IPs/EAs by the frontier orbitals (i.e., HOMOs and LUMOs) can yield good results (especially for semi-empirical methods), this is not necessarily always the case such that suitable DFT functionals or WF methods may be required. [50] The situation is further complicated by the intermolecular environment within the liquid electrolyte. That is, especially the presence of lithium ions with their high charge density polarize the surrounding molecules and thus significantly alter the tendency of the latter to accept an electron. Therefore, by including explicit ions or other solvent molecules in the calculations (in addition to an implicit continuum solvent, leading to the so-called cluster-continuum approach), the accuracy of the results can be greatly improved. [52][53][54] Remarkably, for the reduction of EC, the computed redox potential also varies with the position of the coordinating lithium ion. [52] Note that a similar influence of the intermolecular surroundings is generally observed for the oxidation of liquid electrolytes, that is, the presence of an anion substantially lowers the redox potential of a molecule cluster. [52][53][54][55][56] Interestingly, further studies revealed that either a solvent molecule near an anion can be oxidized, resulting in an energy gain due to the Coulombic interaction between the positively charged solvent molecule and the negatively charged anion, or the anion itself is oxidized when the gain in electronic energy is the dominating contribution. [57,58] As a technical remark, it should be mentioned that for DFT calculations of oxidized/reduced clusters, range-separated hybrid functionals [58,59] or functionals with large amounts of (if not full) Hartree-Fock exchange [52,59] tend to yield more accurate results as they mitigate the self-interaction and delocalization errors inherent to DFT. [60] To summarize, mimicking the intermolecular electrolyte environment beyond an implicit solvation model generally results in an improved theoretical description.
So far, our discussion of the calculation of redox potentials has been limited to molecular properties, that is, the frontier orbitals or IPs/EAs of single molecules or molecule clusters. However, it is important to bear in mind that the definition of the redox potential is a thermodynamic property due to its direct connection to the reaction free energy of the electron transfer reaction via the Nernst equation. [61] Therefore, for a quantitative prediction via DFT calculations, it is crucial to also take electrochemical reactions into account. While these reactions are wellknown for textbook examples or other energy storage devices such as redox-flow batteries, [62] the plethora of reactions occurring during the SEI formation poses a formidable challenge in this context, that is, there is no single reaction that is representative for the entire decomposition mechanism of the carbonate electrolyte (even less so for the long-term ageing of the SEI). Nonetheless, quantitative agreement between computed and experimentally measured reduction potentials can be reached when defluorination reactions (i.e., formation of LiF or HF) are taken into account [52][53][54][55]63] (Figure 2a). Likewise, the inclusion of deprotonation reactions, stabilizing the electron hole, in the computation of oxidation potentials usually yields good agreement with experimental values. [52][53][54][55] Reaction Mechanisms: The initial electron transfer processes near the interface trigger a complex cascade of decomposition reactions. Experimentally, it is difficult to separate these individual processes, therefore, QC calculations offer complementary insights as they allow one to determine the energetics of hypothetical decomposition reactions.
Early DFT studies mainly focused on the identification of probable reaction pathways for the degradation of electrolyte components. Pioneering works studied the reduction and initial decomposition mechanisms of ethylene carbonate (EC), [72][73][74] propylene carbonate (PC) [74,75] or additives such as vinylene carbonate (VC). [74][75][76] In these studies, it was found that EC molecules are essentially reduced in two different steps: First, after transfer of an electron, the ring structure of the EC molecule opens via homolytic cleavage of the bond between an ethylene carbon atom and the adjacent alkoxy oxygen atom (one-electron mechanism), which may dimerize to form either lithium ethylene dicarbonate (Li 2 EDC) in combination with ethylene gas or lithium butylene dicarbonate (Li 2 BDC). After a second reduction step of the EC radical anion, it may be further decomposed into ethylene gas and CO 3 2-(two-electron mechanism). Notably, these works also highlight the role of the intermolecular environment within the electrolyte on the reaction pathways, as it was found that both the cyclic as well as the open singly reduced EC radicals are substantially stabilized by solvation effects, as modeled by an implicit solvent. [72,75,76] Moreover, the inclusion of an explicit lithium ion in the DFT calculations significantly favors EC reduction, although this effect is mitigated if further EC molecules are additionally included to mimic the entire lithium ion solvation shell. [72,75,76] Later, oxidative decomposition pathways of carbonates and other solvents were studied in analogous calculations. [59,[77][78][79] Besides DFT calculations on single molecules or clusters, the initial stages of the SEI formation have been modeled in AIMD simulations employing periodic DFT calculations. [51,66,[80][81][82][83][84][85] Again, the one-and two-electron reduction processes mentioned above could be identified. [51,83] However, another mechanism, in which the ring of the reduced EC radical anion opened within the carbonate moiety itself (i.e., bond fission between an alkoxy oxygen and the carbonyl carbon atom) was observed, forming a diolate in combination with carbon monoxide. [51,83] Further examples of AIMD simulations include the study of protective coating layers [69] or the electrolyte degradation at cathodes. [86] Notably, AIMD simulations have the benefit that the real-time evolution of a thermalized system can be directly studied. In contrast, the static DFT calculations discussed above are based on optimized geometries at T = 0 K, from which thermodynamic information can only be derived from approximate molecular partition sums. In this context, Leung and Tenney extracted the free energy of the lithium ion transfer from a model graphite electrode into an EC electrolyte by means of thermodynamic integration. [87] By variation of the number of lithium ions deintercalated from LiC 6 , the surface charge could be adjusted and the respective potentials of lithium deintercalation and EC degradation were estimated. Further studies focused on the role of additives. For instance, Ushirogata et al. [66] studied the decomposition of the SEI-additive VC (see also Figure 2b). By sampling the free energy as a function of reaction coordinates, they could show that the ring opening of reduced VC differs from that of EC. At the same time, VC can scavenge EC radicals, preventing the second reduction step in this way. [66] In contrast, the decomposition of the additive fluoroethylene carbonate (FEC) mainly leads to the formation of LiF, which adhesively interacts with organic SEI components such as Li 2 EDC, likely resulting in an SEI with improved mechanical properties. [84] Note that analogous studies have been carried out for silicon-anode interfaces. [88][89][90][91][92] In summary, the use of periodic boundary conditions offers the advantage that solvation effects in the bulk electrolyte [66,84,85,93] or surface effects arising from the electrode [51,[80][81][82][83][84] are fully captured (the latter will be discussed in more detail in Section 3.1.2).
Fragmentation Schemes and Reaction Networks: To deal with large numbers of interdependent reactions-as occurring during SEI formation-a scheme to identify thermodynamically viable reactions has been introduced by Husch and Korth. [94] In particular, the fragmentation of carbonate molecules and their recombination was automatically carried out based on a few heuristic chemistry-inspired rules. Based on the reaction energies, a ranking of the most probable SEI components was created. Another metric (or descriptor) for the propensity of a given reaction to occur was the similarity between educts and products-containing information of the required atomic rearrangement. Within this context, a very promising approach is the use of chemical reaction networks (CRNs, Figure 2c). Apart from their general application in many fields of (bio) chemistry, [95,96] they have very recently also been applied to the SEI formation. [68] In this approach, the chemical reactions are mapped onto a graph by creating nodes for all potentially occurring chemical species, which are connected by edges representing the individual reactions. Then, each edge is assigned a 'length' (or cost), which in this case is usually either related to the reaction (free) energy or the respective energy barriers. In this way, it becomes possible to utilize pathfinding algorithms. This approach successfully recovered the well-established formation mechanisms of Li 2 EDC discussed above but also demonstrated that their activation by coordinating lithium ions could be different than previously assumed. [68] In our opinion, this is reminiscent of the dependence of the EC reduction potential on the position of the coordinating lithium ion, which could be addressed in future studies. The main challenge for large CRNs is both to identify hypothetical reaction pathways as well as the choice of a proper cost function. While the identification of potential reactions and their reaction free energies can be automated to a large degree, [68,94,96] including kinetic information in the cost function presently remains a challenge.

Impact of the Electrode Surfaces
While the theoretical approaches discussed in the previous section generally provide insights into electrolyte decomposition mechanisms that are in good agreement with experimental data, the interaction with the electrode surface can be crucial. In principle, interfaces can affect molecular degradation reactions in three different ways: First, the interaction of molecules and ions with the interface-either via electrostatic interactions or short-ranged attraction (i.e., physisorption)-may induce local concentration gradients via accumulation or depletion of certain molecules/ions near the interface, second, large electric fields in the vicinity of the electrode polarize nearby molecules and thus may alter chemical reaction pathways, and third, chemical bonding of molecules or generated fragments to the electrode surface (i.e., chemisorption) may stabilize intermediates or destabilize educt molecules, leading to catalytic effects. In this section, we focus on the latter two, whereas accumulation and depletion effects are discussed in Section 3.2.1.
Electron and Ion Transfer across the Interface: As discussed above, many of the reported AIMD simulations include explicit electrode interfaces, [51,[80][81][82][83][84] thereby allowing to directly simulate the electron transfer from the electrode to the electrolyte without a fictitious electron reservoir. [66,85] In this context, it is, however, important to keep in mind that depending on the progress of the SEI formation, two different regimes of electron transfer from the electrode into the electrolyte can be relevant [69] (cf. Figure 2d): At the pristine electrode surface (initial stage), electron transfer is fast, leading to adiabatic behavior, that is, electrons move significantly faster than the nuclei. This is a common assumption for many QC methods, such that the application of standard DFT methods is uncritical. However, after a thin SEI layer has formed, electrons have to tunnel through this layer to further reduce the electrolyte (unless the latter can permeate through the SEI). Consequently, electron transfer becomes non-adiabatic, and more elaborate techniques such as constrained DFT (cDFT) [70] are required. In addition to these electrochemical processes, it was demonstrated that coatings (or adsorbed solvent molecules) significantly reduce the applied voltage. [97] Recently, the computationally cheaper semi-empirical tight-binding DFT (DFTB) approach has been employed to mimic an entire half-cell with a size of several nanometers and composed of a lithium metal anode, a thin SEI based on Li 2 CO 3 and liquid EC [98] (Figure 2e). In this way, the transfer of lithium ions between electrolyte, SEI and electrode was studied in detail, revealing an asymmetric barrier within the SEI region. Here, the energy barrier required to desolvate the lithium ions in the electrolyte contributes about 0.5 − 0.8 eV to the overall barrier in the range of 0.5 − 2.0 eV, depending on the surface charge density.
Apart from AIMD simulations of bulk-like electrolytes at electrode interfaces, periodic DFT calculations of single molecules at the interface allow one to study the energetics of decomposition reactions in more detail. Furthermore, by predefining the reaction products, slower reactions inaccessible to AIMD can be characterized. For instance, Ebadi et al. studied the decomposition of carbonates on a lithium metal surface, for example, revealing that EC fragmentation forming carbon monoxide (CO) is slightly favored as compared to the release of ethylene gas. [99] Note that similar studies have also been carried out for polymer electrolytes at lithium metal. [100] Interface Reactions, Catalytic Effects, and Electric Fields: To investigate the thermodynamics of (electro)chemical reactions at interfaces in general, the well-known computational hydrogen electrode (CHE) approach [101,102] has been introduced. Although grounded originally in the field of fuel cells and electrocatalysis, [103] it has more recently been applied to interfacial reactions occurring in other electrochemical systems like lithium-oxygen [104,105] and lithium-ion batteries as well. [106,107] Essentially, within the CHE approach, the course of, for example, the oxidation/reduction reaction of elemental hydrogen at a certain catalytic interface is tracked via DFT calculations of the intermediates at the uncharged interface. The impact of the applied potential V is accounted for by a Nernstian term μ e = −eV for the electrochemical potential of an electron in a hypothetical reservoir (such as gaseous hydrogen at normal conditions in case of the CHE). Because the choice of the chemical nature of this reservoir is arbitrary (only the corresponding DFT calculations are system-specific), this approach has been generalized to other systems, [102] including battery systems. [106][107][108][109] For instance, the degradation mechanisms of EC on a lithium metal, lithium oxide (Li 2 O) and lithium carbonate (Li 2 CO 3 ) surface have recently been studied using the computational lithium electrode (CLE) setup. [106] By comparing the pristine lithium surface with Li 2 O and Li 2 CO 3 , serving as model surfaces for passivated lithium metal, it was shown that EC most readily reduces to CO 3 2and ethylene gas on the pristine surface, while on Li 2 O and Li 2 CO 3 surfaces the formation of Li 2 EDC and Li 2 BDC were thermodynamically favored. [106] These findings are in line with the large body of experimental data demonstrating that inorganic degradation (or two-electronreduction) products are formed initially (i.e., at the pristine surface), whereas organic products such as Li 2 EDC and Li 2 BDC are rather formed after the inorganic layer has formed. These findings highlight that the direct interaction of electrolytes with the interface can affect the formed product species. In these cases, modeling approaches going beyond those of Section 3.1.1 via the inclusion of explicit interfaces will yield more quantitative results, albeit at a larger computational cost. Furthermore, the CLE approach has been utilized to study the impact of the applied voltage and the lithiation on the surface structure of lithium metal oxide electrodes. [108,109] Despite its success to describe how an external potential affects electrochemical reactions at interfaces, the CHE (and CLE) rests on the assumptions that redox processes are welldefined, that is, electrochemical reactions are clearly discernible from non-electrochemical reactions. Moreover, in the CHE approach it is implicitly assumed that the free energy of the interface exclusively depends on V via μ e , but not on the charges of the surface atoms induced by V. [110,111] However, the electric field close to the interface can be considerable, which has a strong impact on their electronic structure and thus their tendency to undergo electrochemical reactions. [110,111] Insights into the impact of these effects can be obtained by adding/ removing electrons or ions from the periodic simulation cell that is used to model the interface, which is compensated by a homogeneous distribution of the countercharge. [112] In this way, charged electrode surfaces can be mimicked. These calculations revealed that in addition to the Nernstian electron transfer upon oxidation/reduction, a certain amount of electrons is required to maintain the electrode at the thermodynamic equilibrium potential (i.e., the redox potential). [110] Furthermore, the voltageinduced solvent decomposition at magnesium anodes [111] was investigated in this way.
Another caveat of the original CHE (or CLE) model is that it does not account for solvent effects at the interface. Clearly, for high applied voltages as discussed above, the solvent molecules will become oriented and ions dissolved in the electrolyte will form an electric double layer (EDL) due to capacitive charging. Both effects will at least partly screen the electric field from the electrodes, but likely also alter reaction equilibria between reactants close to or adsorbed at the interface. Indeed, it has recently been observed that competing forces acting on the water molecules in an NaCl electrolyte, resulting from chemical interactions with a TiO 2 surface on the one hand and the double layer on the other hand, gives rise to the non-trivial pH dependence of the capacity. [113] Therefore, modeling the EDL via improved DFT techniques [113][114][115] likely will provide important insights for the initial SEI formation as well.

Simulation of Bulk SEI Phases
The inorganic SEI products Li 2 CO 3 , Li 2 O and LiF are believed to form at the early stages of the SEI formation in close proximity to the electrode interface. For a functional SEI, this picture therefore implies that these inorganic compounds should have good ionic conduction properties, while preventing electronic conductance. Here, Chen et al. showed via density-ofstates calculations that Li 2 CO 3 , Li 2 O and LiF exhibit significant band gaps, and thus can be considered as electronic insulators. [116] Creating lithium vacancies in Li 2 CO 3 led to the formation of electronic states close to the Fermi level, which were, however, strongly localized and thus likely do not affect the insulating properties of the material.
On the other hand, the ionic conduction was found to be fast in Li 2 CO 3 and Li 2 O due to the presence of vacancies, whereas LiF is a poor ion conductor. [116] The good ionic conduction of Li 2 CO 3 was also shown by Iddir et al. [117] Later, Shi et al. investigated various lithium diffusion mechanisms in Li 2 CO 3 in more detail, showing that the main mechanism for lithium interstitials is a so-called knock-off mechanism, in which the interstitial lithium atom displaces a neighboring lattice atom into a vacant interstice, [118] rather than a hopping process between distinct interstices (see Figure 2f). Notably, the concentration of these interstitials-formed by excess lithium ions-or other lattice defects, for example, vacancies, depends on the voltage applied to the electrodes coated by the Li 2 CO 3 -containing SEI, such that the ionic conductivities of Li 2 CO 3 at the anode and the cathode may differ. [71] Furthermore, also neutral lithium atoms may diffuse via this interstitial mechanism, giving rise to electron leakage and thus ongoing electrolyte decomposition during the battery cell's lifetime. [71] Finally, if LiF is present in the SEI, leading to LiF/Li 2 CO 3 interfaces, it is thermodynamically favorable that lithium ions are partly transferred from LiF to Li 2 CO 3 , leading to vacancies in LiF and excess interstitials in LiF/Li 2 CO 3 , which significantly enhances the ionic conductivity in the space charge layer region, [43] thus rationalizing the beneficial effect of LiF [119,120] (which itself is a poor ionic conductor [116] ).

Shortcomings
Although DFT calculations can provide valuable insight into decomposition mechanisms leading to the formation of the SEI, essentially two main drawbacks exist: First, DFT calculations scale unfavorably with the system size, making the treatment of larger systems and longer time-scales (such as in AIMD simulations of bulk systems) computationally very demanding. Consequently, with AIMD, only very fast reactions (i.e., occurring on the time scale of up to 100 ps) can typically be accessed, leading to poor statistical sampling. For the same reason, liquid carbonate electrolytes are often modeled with EC as the only solvent, although in practice solvent blends typically contain at least two different carbonates as well as additives. Therefore, the limitation in system size also hinders the exploration of the full concentration range, which likely is an important parameter to understand or even direct the SEI formation. This problem could be solved by ML-force fields, although multi-component systems would still be challenging (see Section 4.1).
Second, the search for potential intermediates or reaction products is essentially based on the chemical intuition of the researcher or inferred from experiments, at least in naive implementations. This substantially reduces the amount of both electrolyte compounds and putative reaction pathways that can be investigated. As a result, it is by no means guaranteed that the true/full reaction mechanism is uncovered. Although heuristic fragmentation-recombination schemes and CRNs [68,94] are very promising to overcome this issue, the inclusion of kinetic information remains challenging, especially when carried out in an automated fashion. Here, kinetic descriptors, [96] for example, based on Hammond's postulate [121] or the Bell-Evans-Polanyi principle-stating that reaction energy and activation energy are correlated [122,123] -may be used in future works. In Section 4.2, we summarize how ML models can work together to address some of these open questions.
In summary, we believe that employing AIMD simulations with an enhanced sampling of reaction events without imposing reaction coordinates that require prior knowledge of the mechanisms-such as nanoreactor schemes [124,125] -will help to overcome both limitations mentioned above and to tackle chemically reactive systems with a large amount of possible reactions. For instance, such a nanoreactor could operate by exerting pressure pulses to a molecular system, which accelerates the crossing of reaction barriers. [124] Additionally, Grimme used the well-known metadynamics approach [126,127] in combination with the fast DFTB method, in which the root mean squared deviation (RMSD) from the initial reactant geometries was used as reaction coordinate. [125] In this way, the reactants become unstable and reactions are enforced, but without necessarily targeting them toward a product state that has to be known beforehand. Notably, both approaches have been successfully applied to reproduce the Urey-Miller experiment-similar in complexity as the SEI formation.

Classical Molecular Dynamics
Classical MD simulations rely on force fields in which the quantum-mechanical description of the electrons is replaced by empirical potentials. This approach therefore is computationally significantly faster, however, due to the neglect of explicit electrons, chemical reactions and electron transfer processes are generally not captured. Consequently, SEI modelling via classical MD simulations has mainly been carried out either (i) prior to the actual SEI formation, that is, to investigate the molecular structure of the liquid electrolyte in the vicinity of charged interfaces (Section 3.2.1), (ii) for the bulk of idealized model SEI structures present after the formation process (Section 3.2.2), or (iii) the interface between such model SEIs and the liquid electrolyte (Section 3.2.3). Furthermore, attempts have been made to introduce chemical reactivity into MD simulations based on empirical force fields-either by explicitly deriving potentials that account for the breaking/formation of bonds or electron transfer (reactive force fields as, e.g., ReaxFF) or by regularly interrupting the MD simulation and converting reactant into product molecules (to be outlined in Section 3.2.4). Note that other reactive force fields than ReaxFF exist, for example, [128,129] which however have not yet been applied to the field of battery science to the best of our knowledge.

Solid/Electrolyte Interface Simulations
The electrolyte structure at the electrode interface is usually studied by simulating an electrolyte slab confined between two charged or uncharged model electrodes (often considered to be planar). With classical force fields, the thickness of the electrolyte slab is typically on the order of 10 nm, whereas the surface area is in the range of 10 nm 2 . [130][131][132][133][134] These dimensions allow for a realistic description of the bulk electrolyte beyond the interface layers and to investigate a possible structural ordering of the molecules and ions that are adsorbed at the interface.
Special care has to be taken to model the external voltage applied to the electrodes. In the most basic implementation, one simply imposes a constant surface charge, that is, each surface atom of the electrodes bears a partial charge (opposing in sign for the two different electrodes) to mimic the voltage of a plate capacitor in vacuo. However, the ion accumulation at the interface generally polarizes the (conductive) electrode and thus affects the distribution of effective surface charges. As a more realistic description, Siepmann and Sprik developed a model [135,136] that accounts for this effect by adjusting the (Gaussian) charges of the electrode atoms self-consistently in such a way that the electrostatic potential at each electrode atom position equals a preset value (i.e., the half-cell voltage). In fact, simulations showed that this constant-voltage approach results in an improved description of the structure of the interface layers as compared to the constant-charge model, [137,138] although these differences are most pronounced at high voltages. [138] Furthermore, the EDL formation dynamics [139] and the associated Joule heating [137] as well as the lithium desolvation process in ether-based electrolytes [140] are affected by the choice of the electrode model. It should be mentioned, though, that due to the iterative determination of the charges of the electrode atoms, the constant-potential method is computationally more expensive than the constant-charge model. Here, Petersen et al. introduced a computationally cheaper polarizable electrode model based on image charges, [141,142] however, this model is limited to planar electrode geometries and can hence not be employed to describe microporous or patterned electrodes. [137,143] Due to the fact that electrochemical processes are generally not captured by classical force fields, classical MD simulations have largely been applied to model materials for non-faradaic energy storages such as supercapacitors. [130,[132][133][134] Nevertheless, in the case of liquid battery electrolytes based on solvent blends (e.g., EC and DMC), it has been demonstrated in MD simulations with polarizable electrodes that EC is preferentially enriched at charged interfaces as compared to DMC due to its larger dipole moment [131] (Figure 3a). For the same reason, the EC molecules become oriented: While for the uncharged interfaces, the carbonyl group of EC tends to be oriented flatly with respect to the surface, applying a voltage forces the carbonyl groups (or the aliphatic carbons) of EC to orient toward positively (negatively) charged electrodes. [142] In total, these orientation effects are somewhat more pronounced at basal graphite surfaces as compared to the edge planes-in particular at high potentials, which may thus impact the precise SEI formation mechanism. [142] It is important to note that in this way, even though classical MD simulations do not explicit account for it, one can nonetheless derive information about a system's reactivity, as due to the accumulation of EC at the interface, it will presumably react most readily. Furthermore, if the structure of the interface layer is known, one may extract representative clusters from MD simulations and subsequently subject them to DFT calculations to compute their redox stability or putative decomposition mechanisms (Section 3.1.1).

Bulk SEI Properties
Another application of classical MD simulations is concerned with the transport properties within the SEI after its formation, which is crucial to understand the interfacial resistance of the SEI layer. Although the main SEI constituents are essentially known (cf. Section 3.1.1, despite ongoing controversy [36] ), its precise structure-for example, the degree of crystallinity, the layering of different compounds or morphological propertiesremains largely unknown, which poses a serious challenge for all modeling approaches. However, classical MD simulations can be employed to compute the ion coordination and lithium transport in various model SEIs or pure SEI components.  [131] investigated the voltage dependence of the density profiles of the different ionic species near the anode. Adapted with permission. [131] Copyright 2011, American Chemical Society. b) Bulk SEI properties. Bedrov et al. investigated the temperature dependence of lithium transport within ordered Li 2 BDC and Li 2 EDC materials and their disordered analogues. [144] Adapted with permission. [144] Copyright 2017, American Chemical Society. c) The SEI/electrolyte Interface. By correlating the free energy landscape of lithium ions with the solvation structure within an SEI consisting of lithium ethylene dicarbonate (Li 2 EDC) and lithium carbonate (Li 2 CO 3 ), Jorn et al. identified two stages for the insertion of lithium ions from the electrolyte to the SEI, namely adsorption and absorption. [145] Such insights are particularly relevant for distinguishing between the entangled effects of lithium exchange with and transport in the SEI on the overall ionic conductivity. Reproduced with permission. [145] Copyright 2020, American Chemical Society. d) Reactive molecular dynamics. Left image: Reactive simulation of EC radicals with explicit electrons within the eReaxFF model. [146] Adapted with permission. [146] Copyright 2016, American Chemical Society. Right image: While the studies sketched in (a)-(c) were focused on extracting lithium ion transport properties from empirical SEI structures, Takenaka et al. have attempted to determine the chemical composition and morphology of the SEI resulting from a QM-based reaction network and starting from a homogeneous salt and solvent distributions. Adapted with permission. [147] Copyright 2014, American Chemical Society.
For instance, Bedrov et al. [144] studied Li 2 EDC and Li 2 BDC in the bulk and compared ordered with amorphous structures (Figure 3b). They could demonstrate that ordered structures show improved transport properties and that these materials are basically single ion conductors at ambient temperatures. [144] Notably, the activation energy for ion diffusion in Li 2 EDC and Li 2 BDC was significantly larger than the estimated barrier associated to the exchange between SEI and liquid electrolyte [42] (see Section 3.2.3), identifying potential bottlenecks for the overall lithium migration in this way. However, earlier work revealed that liquid carbonate electrolytes plasticize the SEI, leading to significantly shorter lithium ion residence times. [148] Such effects may especially prevail in the interface region between SEI and liquid electrolyte (see also Section 3.2.3).
Clearly, the main caveat of these simulations is that the SEI structure has to be entered in an ad hoc manner. Here, reactive MD schemes, to be discussed in Section 3.2.4, seem a promising avenue to overcome current limitations. At the same time, the possibility to focus on well-defined structures also allows one to study limiting cases that may be difficult to realize experimentally.

The SEI/Electrolyte Interface
Finally, classical non-reactive MD simulations have been employed to simulate the interface between the SEI and the liquid carbonate electrolyte. [42,142,145,149] Even though the studied SEI layers were rather thin (up to a few nanometers) due to the limitations in system size, a roughened SEI/electrolyte interface was generally observed for equilibrated structures, which was associated with the intrusion of EC from the liquid electrolyte. [142] Consequently, also the exchange of lithium ions between the electrolyte and the SEI occurs gradually, that is, the lithium ion coordination shell in the electrolyte, consisting of EC and PF 6 − anions, is continuously replaced by EDC 2 − groups from the SEI. [149] This is in stark contrast to the observations from MD simulations of interfaces without SEI, here, lithium ions only come into contact with the electrode surface at rather high applied potentials. [131,142] As a result, the energy barrier for the ion transfer process between SEI and carbonate electrolyte is significantly lower than for the diffusion within the SEI. [42] More recent MD simulations with umbrella sampling [145] revealed that the energy profile for the insertion of Li + into the SEI differs from the opposing release of Li + from the SEI (Figure 3c), mainly due to dissipation related to the reorganization of Li 2 EDC in the SEI. Apart from the SEI/electrolyte interface, Tasaki used classical MD simulations to study the adhesion of Li 2 EDC and lithium 1,2-propylene glycol dicarbonate (Li 2 PDC, formed analogously from PC) at the graphite anode, with the outcome that Li 2 EDC tends to stick slightly stronger to the anode interface. [73] As described in Section 3.2.2, the main drawback of these simulations is that it is usually required to make ad hoc assumptions about the structure and composition of the SEI, which are extremely hard to validate. Nevertheless, in this way well-defined scenarios can be constructed to assess their validity in light of experimental data.

Reactive Molecular Dynamics
Due to the importance of both reactions and diffusion processes, several attempts have been made to model the SEI formation in reactive MD schemes. Essentially, this includes two different approaches: first, the use of reactive force fields such as ReaxFF, [150][151][152][153][154][155] in which bond breaking and formation is-as for non-reactive force fields-described by empirical potentials, and second, hybrid MD/Monte Carlo (MC) schemes, in which the molecular topologies of reactants are changed in frequent intervals. [147,[156][157][158][159][160][161][162][163][164] ReaxFF: In this case, the change in connectivity is achieved by using the bond order of the individual atoms as an additional variable. If any element forms fewer or more bonds than according to its valence number, the bond order potentials lead to an energy penalty. [165,166] In this way, chemical reactions, such as the initial degradation of electrolyte molecules at the anode interface, can be directly studied. For instance, Kim et al. studied the formation of the SEI in carbonate electrolytes with metallic lithium and could confirm its mutlilayer structure with an inner layer composed of inorganic salts and an outer layer containing organic decomposition products such as Li 2 EDC. [150] Similarly, Bertolini and Babluena observed that a porous SEI layer immediately forms in ether-and carbonate-based electrolytes at a lithium metal interface. [155] Interestingly, the pores of this layer were not exclusively filled by decomposition products, but also contained unreacted electrolyte molecules. [155] Moreover, Yun et al. showed that dangling silicon atoms at silicon anode surfaces may catalyze two different ring-opening reactions of EC, leading to the release of carbon monoxide and ethene, respectively. [154] Here, the authors demonstrated that this decomposition reaction may be inhibited on oxidized silicon surfaces, on which Li 2 O is formed with a higher propensity, or via additives. [154] Other ReaxFF studies from the field of battery science focused on the effect of anode coatings, [152] the phase stability of lithiated sulfur cathodes, [153] reactions of EC radicals in bulk carbonate electrolytes [151] or electrolyte degradation at cathodes interfaces. [167] The original ReaxFF implementation [165] suffered from the drawback that the atomic charges were determined by the Electron Equilibration Method (EEM) [168,169] (similar to charge equilibration [170] ), which has several disadvantages, especially for the description of electrochemical reactions. [166,171] In particular, EEM leads to spurious long-range charge transfer which is inconsistent with discrete changes in atomic oxidation states or electron transfer reactions occurring between neighboring atoms only. [171] In other words, due to the fact that a non-integer charge may be distributed across two distant atoms, dielectric materials behave too 'metallic'. More recently, this problem has been circumvented by implementing the ACKS2 scheme, [172,173] which prevents long-range charge transfer by penalizing changes in the atomic electron population between remote atoms, similar to the split charge equilibration method. [174] Another drawback of the original ReaxFF model [165] is that electrons are treated implicitly, prohibiting simulations of electron transfer processes such as oxidation and reduction. To overcome this deficiency, the eReaxFF formalism [175] treats an excess electron or an electron hole in a molecule as a (pseudo) particle. Here, the electron is represented by a Gaussian wave function that interacts with the point-like nuclei. Furthermore, the presence of an electron or hole is coupled to the over-and undercoordination terms in ReaxFF. Indeed, this framework could successfully describe the reduction of EC as well as the cascade of consecutive reactions [146] (left snapshot in Figure 3d), although these decomposition mechanisms have to be carefully parametrized. [176] Notably, the latter work revealed a notable impact of the lithium ion solvation shell on the EC degradation. Furthermore, the authors suggest to use two different parameter sets-one for metallic lithium and one for lithium ions-to incorporate distinct oxidation states. Similar as in hybrid MD/MC schemes (see below), the charge state of the lithium atoms can be altered in MC steps. [176] This is similar in spirit to simulations of a model battery cell without an explicit treatment of electron particles. [177] ML models that can universally track electron density and its transfer and modulation for large scale reactive simulations can usher in a new era (see Section 4.1).
Hybrid MD/MC Approaches: An alternative way of accounting for bond breaking or formation mechanisms consists of the reactive step MD [163,164] or similarly the hybrid MD/MC [147] simulation. These schemes both consist of performing a classic MD simulation that is periodically interrupted in order to insert or remove molecular bonds according to the definition of the involved reactions. The conditions for the insertion or removal of bonds can be subdivided into two types, namely geometric conditions, where the spatial proximity of the reactants is verified, and dynamic conditions, corresponding to the probability of a reaction to occur per unit time given that the geometric conditions are met. The geometric conditions are directly derived from classic force field parameters, since these mainly consist of inter-molecular distances. On the other hand, the dynamic conditions are derived from quantum-mechanical calculations of the energy landscape associated with the corresponding reactions. Within the hybrid MD/MC scheme the reaction probability is defined based on the energy difference between the initial and final states of the reactants. In contrast, in the reactive step MD scheme the reaction probability is estimated from the energy difference between the initial state and the transition state (or, alternatively, by directly determining the reaction rates from AIMD simulations if the reactions are sufficiently fast [164] ). While the hybrid MD/MC scheme appears to be computationally effective, a severe shortcoming of this method is the possible overestimation of reaction rates involving highly energetic transition states. The reactive step MD scheme is therefore expected to yield more accurate insights into dynamic mechanisms related with the SEI formation. Despite the conceptual simplicity of the reactive step MD scheme, it should be noted that the operation of bond rearrangement may involve non-trivial difficulties related with the smoothness of the resulting force field. The suggested approach to avoid such non-smooth variations mainly consists of performing local relaxations of the reacted molecules after each reaction. [164] In fact, the physical consistency of this method with regard to structural and dynamic system properties has been verified by means of a comparison with ab-inito MD computations performed on a relatively short time scale, [164] thus providing a theoretical basis for related future computations.
Coming to a more applied perspective, a number of investigations of the initial stage of SEI formation based on the hybrid MD/MC has been recently reported. [147,158,178] These studies have mainly attempted to determine the spatial distributions of the various reaction products formed due to solvent reduction near the anode and subsequent dimerization of the reduced solvent molecules. By performing this procedure for a wide range of salt concentrations, Takenaka et al. [158] attempted to correlate the SEI thickness with the salt concentration, suggesting that an increase of salt concentration may be related to a decrease of the SEI thickness due to slower diffusion of salt driven reaction products. In a closely related work, [147] the hybrid MD/MC method has been used to investigate the dependence of the morphological properties of the SEI on the structural properties of the solvent molecules or the role of additives. [156,157] A more recent study [178] was focused on the formation of the multi-layer structure of the SEI related with the decomposition of lithium hexafluorophosphate salts and ethylene carbonate solvents in contact with an amorphous silicon anode.

Shortcomings
The inherent limitation of both non-reactive and reactive force fields is twofold: First, the force field parametrization can be laborious, making the inclusion of novel compounds difficult. Generally speaking, the challenge of this task increases with the complexity of the simulation model. For instance, although polarizable force fields are considered to be more accurate (especially for charged systems), [179][180][181][182][183] additional parameters like atomic polarizabilities enter the model, which need to be determined. This is even more true for ReaxFF as its parametrization is not based on atom types (e.g., an ether oxygen atom or a carbonyl oxygen atom of EC), but rather contains several potential and interpolation functions. This leads to a substantially large number of required force field parameters per element, which requires special fitting techniques that can make the parametrization elaborate. [165,166] Besides the force field itself, hybrid MD/MC approaches require additional parametrization for the individual reactions. As discussed in Section 3.1.4, this is particularly problematic when the potentially occurring reactions are unknown. Here, ReaxFF should in principle be able to uncover such reactions, while hybrid MD/MC approaches are computationally cheaper.
A second inherent limitation of conventional force fields is that the functional form of the interaction potentials is predefined, such that even the most thorough calibration or the inclusion of additional interaction terms not necessarily capture all relevant effects. One particular benefit of ML-based force fields [184][185][186] is that they do not suffer from this drawback due to the fact that any form of the interaction potentials can be captured by the ML-based interpolation schemes. Recent developments in ML models even let us track electron transfer during simulations. These aspects will be further discussed in Section 4.1.

The Continuum Scale
Physical models use a group of partial differential equations to describe the physical processes that occur at the interfaces and through the porous structure of the various compositions of batteries. However, physical models can be computationally complex and are of limited fidelity due to the lack of directly measured physical parameters. [187][188][189][190] Computationally feasible approaches are adopted to reduce the complexity of such models (e.g., ignoring certain processes related to the degradation mechanism). There are different types of physical models for simulating battery interfaces at different length scales that offer varying levels of system description details and computational complexity. For instance, the use of molecular dynamics (MD) [82,93,191] and density functional theory (DFT) [192][193][194] to study the effect of the SEI layer on the intercalation of Li ions (Section 3.1). For high fidelity continuum scale simulations one can use phase field models to simulate the formation [195] and the micro-structure morphology evolution [196] of the SEI in lithium ion batteries. The investigation of the formation and growth of the SEI film can be done also with meso-scale interfacial modeling [197] and multi-scale modeling approaches [11] at a higher length/time scale.
While all the above-listed approaches are well exploited, continuum models are mostly used to simulate the electrochemical performance of batteries due to their lower computational cost. Continuum models are usually coupled with physical degradation models such as the formation and growth of the SEI during the charging/discharging process of batteries. However, the use of continuum models requires the assumption of a homogeneous behavior of the cell, the formation and structure of the SEI and reaction occurring at the SEI/ electrolyte interfaces. However, the formation of the SEI is a complex inhomogeneous process thus making the predictions made by continuum models phenomenological with a significant level of uncertainty. Continuum models come in different complexities depending on the assumptions and length scales (a) microscale: particles of the electrode active materials (b) mesoscale: porous electrodes, and (c) macroscale: cell design. The various assumptions, continuum models and how they are coupled with the physical degradation model are described in the sections below in the order of increasing complexity.

Single Particle (SP) Models
The single-particle model (SPM) developed by Zhang et al. [198] incorporates the effect of transport phenomena within electrodes in a simplified manner (Figure 4). In SPM, diffusion and intercalation occur within a single electrode particle. The SPM can be applied to study the electrochemical performance of cells at the sandwich level by considering the anode and cathode each as a single particle with the same surface area as the electrode. [203] The major drawback of the SPM is that it does not consider the effects of concentration and potential in the solution phase between the particles. [190,203,204] This makes it inappropriate for simulating the electrochemical performance of cells at high current rates where the effect of the potential in the solution phase is significant. [204] These limitations can be alleviated by using high order polynomial functions derived in an ad-hoc manner [205][206][207] and in a systematic asymptotic manner [208] for the solution phase concentration and potential within the electrodes, and parabolic function for those within the separator. The SPM is simple and does not require much computational effort which makes it suitable for lifetime modeling of Li ion batteries by incorporating the effects of the degradation mechanisms such as the increase in the film resistance and the loss of active Li ions during the formation and growth of the SEI [209][210][211] but at limited conditions of low current rates and thin electrodes. [203] The formation of the SEI in SPM is assumed to be a singlestep rate reaction involving the reduction of ethylene carbonate into one of the main products such as Li 2 CO 3 and the growth is assumed to be homogeneous in the axial direction. These assumptions make it easy and computationally cheap to predict the effect of the SEI on the cell performance but deviate from the complex structure and heterogeneous growth of the SEI observed experimentally. [4,33,212] Colclasure et al. [213] on the other hand incorporated detailed chemical kinetics and multicomponent species transport into the SEI film growth model coupled with a single graphite particle model surrounded with liquid electrolyte. Even though several SEI kinetic reactions were considered in this model, the growth of the SEI was studied with a 1D model and based on the homogeneous growth assumption.

Pseudo-2D (P2D) Model
The pseudo-2D (P2D) model ( Figure 4) is an extension of the ohmic porous-electrode models, which accounts for the solid and solution phase potentials and current, but neglects spatial . Computational demands and system complexity and resolution of continuum models. The single particle (SP) model [190,198] is suitable for simulations at low current rates if insignificant changes in the concentration of the electrolyte are expected. The pseudo-2D (P2D) [190,199] model describes the electrochemical reactions by assuming a Butler-Volmer kinetics and includes mass transfer in the solution and solid phases. The pseudo-3D (P3D) multiscale model [200] provides a better understanding of the dynamic operation of Li ion batteries for large-scale application where the generation of heat is critical. The SP, P2D and P3D models focus on the homogeneous formation and growth of the SEI. Using phase-field models [201] the heterogenous formation of the SEI can be modeled using different phase diffuse parameters to represent the various components to the SEI. Image of Phase-field model: Reproduced with permission. [201] Copy right 2018, American Chemical Society. 4D-resolved model can be used to form heterogenous SEI layers and predict their spatial distribution. [202] Image of 4D-resolved model: Reproduced under the terms of the CC-BY 4.0 license. [202] Copyright 2021, The Authors, published by Wiley-VCH.
variation in the concentrations. The ohmic porous-electrode model assumes either a linear, Tafel, or exponential kinetics to describe the electrochemical reaction and incorporate additional phenomena such as the dependency of conductivity as a function of porosity. [214,215] Unlike the ohmic porous-electrode model, the P2D model describes the electrochemical reactions by assuming Butler-Volmer kinetics and includes diffusion in the solution and solid phases. The P2D model was developed by Doyle et al. [199] to describe the galvanostatic charge/discharge of a lithium anode/solid polymer separator/insertion cathode cell using the concentrated solution theory. The generic nature of the P2D model permitted the incorporation of a variety of battery compositions which lead to the development of a number of similar advanced models. [203,[216][217][218][219] The predictive capability of the P2D model is better than the SPM due to the addition of many internal variables. However, this comes at an extra computational cost than the SPM because of the non-linear Partial Differential Equations (PDEs) employed to describe the transport phenomena of various species, electrochemistry and thermodynamics in batteries. The use of reduced-order models (ROM) where these non-linear PDEs are reduced to simpler Differential Algebraic Equations (DAEs) [190] has been an effective means of solving the issues related to the computational cost and complexity of P2D models.
Christensen and Newman [220] adopted the P2D model to develop a continuum scale mathematical model to simulate the growth of the SEI and transport of Li ions and electrons through the SEI film on the negative electrode. The reaction mechanism used for the formation of the SEI was kept simple by adopting just a portion of the mechanism proposed by Aurbach [212] for the formation of Li 2 CO 3 while neglecting the formation of (CH 2 OCO 2 Li) 2 . Even though the P2D model is based on porous composite electrodes, the SEI models developed are assumed to be a 1D film growing on a planar electrode surface. This modeling approach has been used by many researchers even in this current era to estimate the rate of the growth of the SEI, the film resistance and the irreversible capacity loss due to the formation of the SEI. [221][222][223][224] However, the major question remaining is whether Li 2 CO 3 is the only major constituent of a real SEI in Li ion battery and whether the growth of the SEI is homogeneous in one direction.

Pseudo-3D (P3D) Multi-Scale Model
The P3D multi-scale model is a 1D + 1D + 1D electrochemical and thermal model developed to simulate the transport phenomena of species, multicomponet mass and charge transfer, and variations in the temperature gradient in Li ion batteries (Figure 4). Each dimension represents a specific spatial scale in a hierarchical manner. The first dimension represents the nanoscopic level where the transport of the electroactive species in the particles of the active material is governed by diffusion. The second dimension is the microscopic level where the description of multicomponent mass and charge transport as well as heat production occurs in a single repeating unit comprising the anode, separator, cathode, and current collectors. The last dimension is the macro-scale level which describes the heat transport across the cell. The three scales are coupled through respective boundary conditions. [225] The P3D model was developed mainly to provide a better understanding of the dynamic operation of Li ion batteries for large-scale application where the generation of heat is critical. However, this model is computationally expensive and involves the consideration of several assumptions and approximations leading to a variety of shortcomings such as the inability to keep track of the thermal effects of electrochemical parameters [226] and the need for empirical input from experiments or use of volume-averaged equations for solid-phase intercalation. [227] A detailed P3D model that considers the layered structure of cell stacks, the case of a battery pack and the gap between both battery components was developed by Chen et al. [227] to analyze the thermal behavior of Li ion batteries. Herein, a simplified thermal model was assumed to speed up the calculation. A more reliable approach to simulate the heat production is to express the molar enthalpies and entropies as a function of the state of charge [200] but it comes at a computational cost. The P3D model was extended to include a detailed chemical reaction and a multi-phase chemistry to predict the thermal and ageing effects of Li ion batteries. [223] Owing to the complexity of the extended model, a time-up-scaling methodology was adopted to accelerate the simulations. The degradation mechanism considered in this model was a single-step charge transfer reaction for the formation of SEI and the growth was assumed to be limited by diffusion and homogeneous. However, the main SEI component assumed in this model was (CH 2 OCO 2 Li) 2 and not Li 2 CO 3 as assumed in some of the P2D models and the SP model. This raises a possible concern about the component of SEI that should be considered in modeling the degradation mechanism in continuum models and the assumptions concerning the homogeneous growth of the SEI.

4D-Resolved Model
4D-resolved models are 3D in space models coupled with time to simulate the performance of electrochemical systems ( Figure 4). The 3D model part involves the generation of 3D spatially resolved microstructures. These 3D microstructures are usually computationally generated via coarse-grained molecular dynamics (CGMD) using LAMMPS software [228] or random location algorithms such as GeoDict software. One major limitation of existing 3D resolved models is the use of effective parameters for porosity and tortuosity to implicitly represent the carbon-binder domain (CBD) [229,230] or merging the CBD with the active material as a single solid phase. [231] This inhibits the capturing of the effect of the loss of the actual active surface area and the electron heterogeneities due to the percolation network of the CBD since the initiation of the SEI formation is mostly at the graphite/electrolyte interface.
Recently, a computationally efficient and versatile algorithm capable of meshing multiple phases via voxelization was developed by A. Franco et al. [232] This algorithm permits the generation of multiphase geometries with distinct CBD and active material phases filled with electrolyte. The coupling of the multiphase 3D-microstructures with time (4D-resolved model) has been implemented to assess the impact of CBD spatial location and its transport properties toward Li ions on the electrochemical response, and the identification of spatial operation heterogeneities inside the electrode. [233] 4D-resolved models have also been applied to assess the heterogeneous growth of the SEI and its impact on the electrochemical performance of the graphite porous electrode. [202] Even though the 4D-resolved model was capable of forming heterogeneous SEI layers and predicting their spatial distribution, there is still more room for improvement for the description of the SEI formation step in the model since it does not consider the various components of the SEI.

Phase Field Models
Phase-field models are computational models which describe the micro-structure evolution of interfaces of material systems as a function of space and time (Figure 4). Phase-field models have the diffuse interface feature, which introduces one or more phase-field variables to account for the state of the phase of the system. The evolution of the phase-field variable is directly related to the phase transition process, hence there is no need for phase-boundary tracking. [234] This reduces the computational time required for a given problem, but is comparatively higher than that of the P3D and the subsequent lower-dimensional models. However, more phase-filed variables are required for systems that involve concentration, temperature and applied electric field. The derivation of the evolution equations of phase-field variables is based on the principles of local equilibrium [235] and free energy minimization, [236] which are normally non-linear partial differential equations. Current phase-field models describe the kinetic phenomena in batteries where the electrolyte is usually treated as an ideal solution. These kinetic phenomena only require a reasonable (but not realistic) thermodynamic data-based free energy functional of a practical electrochemical system. This inhibits the currently available phase-field models from reflecting on the individuality of the investigated system or interphase in the SEI.
The formation of the SEI during the charging process of batteries usually involves the electrode and the electrolyte (two phases). This can be successfully handled by current phase-field models developed for Li dendrite deposition and growth [201,[237][238][239] if the diffusion of lithium ions through the SEI layer during (de) intercalation is neglected. However, most of these models are based on linear kinetics and are not applicable to problems far from equilibrium. This can be resolved by developing nonlinear phase-field models which incorporate the Butler-Volmer reaction kinetics to predict the growth of the lithium deposit, [240] and invariable thermodynamic phase-field models to analyze dendritic patterns. [241] But, the effect of the structure and composition of the SEI on the dendritic growth was not considered in these models. This can be rectified by extending the existing phase-field models for Li deposition by introducing a noise field to the solid Li/ electrolyte interface to account for the effect of the formation of the SEI [242] and using numerical reconstruction method to model the geometric micro-structure of the SEI. [243] A major limitation of this extended model is that the introduced noise field depends on some artificial assumptions for the diffusion across the SEI which raises questions about the uncertainty regarding the composition of the SEI and the type of electrolyte.

Shortcomings
The main limiting factors hindering the use of continuum models to describe the formation and growth of the SEI are the assumptions made-SEI structure and growth being homo geneous and the main product of the SEI considered. These limitations can be overcome by considering the heterogeneity of the SEI using phase-field models to describe the multilayer mosaic structure and growth of the SEI and using different phase-field variables to describe the phase changes of the different SEI products. This will require the modification of existing phase field models to a highly meticulous and concordant phase-field model for the multi-component multiphase SEI where the effects of (de) intercalation reaction during the operation of batteries can be considered. Obviously, more input parameters will be needed for such models and critical attention should be given to the way they are obtained. These input parameters can be obtained via experiments, MD simulations, DFT calculations and can be treated as fitting parameters. However, the impact of each of these input parameters on the output para meters (for instance capacity loss and thickness of SEI) will vary in both space and time. Such phase-field models are expected to be computationally extremely expensive. Using ML surrogate models might help in reducing the computational load (see Section 4.3). Additionally ML models can be used to conduct uncertainty analysis to ascertain the impact of the input para meters on the output parameters (see Section 5.1).

A Way Forward with ML
Two distinct kinds of machine learning models hold particular promise in the battery interface simulation domain. Supervised models can be used for regression tasks, where traditional physics-driven models are computationally very expensive or need sequential numerical solutions (e.g., molecular dynamics). Moreover, models can learn to predict properties derived from time dependent simulations directly or to expedite those via reliable approximations. We will discuss some of the most promising future directions in the application of ML in SEI simulations.

Machine Learning Potentials
By construction of ML force fields, the overarching goal is to achieve the accuracy of DFT calculations with the computational efficiency of classical force fields, [13,186,244] allowing an appropriate size of the simulation systems and the time scales needed for proper statistical sampling of many interface properties and quantities of interest. These models should learn the correlation between structure and energy/forces, that is, the relevant interactions purely from data obtained by running actual simulations are in principle only limited by the quality and quantity of the reference data used to train them. This also captures the current limitation of ML potentials in interface simulation. Another benefit of ML-based force fields is that they do not depend on any predefined functional form.
Therefore, also complex correlations can be learned. Although it is in principle possible to augment classical force fields with additional terms to capture more complex interactions (such as polarization [183] ), this usually leads to an increased computational cost. In contrast, the computational cost of ML-based force fields only depends on the choice of the model, but not the complexity of the interactions. The state of the art models capable of handling complex chemistry can be classified into three main categories ( Figure 5 )-a) Gaussian process-based learning [245] with efficient structural fingerprints like SOAP (Smooth Overlap of Atomic Positions). [246] These are very data efficient, but difficult to generalize when the system contains many types of atoms and the dataset size required to sample the phase space is very large. b) Neural network potentials [247] operating on predetermined structure fingerprints of the local environment. [248] These models can handle large chemical complexity and large data sets help but the locality of structural fingerprinting often limits how long-range interaction can be learnt. c) Graph neural networks [249,250] that learns both the fingerprints as well the mapping from that to the energy and forces. These are the most recent developments and particularly well suited for SEI modelling, due to their high flexibility, the possibility of learning long-range interaction via message passing, [251] automated learning of the most efficient fingerprints and the graph data structure being native to the atomic structure ( Figure 5). However, it is important to keep in mind that even highly accurate ML force fields will require enhanced sampling techniques such as metadynamics [126,127] (see Section 3.1.4) to tackle slow processes occurring on scales of microseconds or longer.
Building ML potentials for solid-state electrode material simulations [20,24,252] is easier than for the SEI, as key properties can be obtained with an ML potential even if it is only accurate in the vicinity of the equilibrium structure and only captures short-range interactions. On the contrary, for effective SEI simulations, ML potentials need to handle long-range anisotropic electrostatic interactions and polarization, as well as be accurate in far-from-equilibrium structures that constitute reaction pathways. Thus the real application of ML potentials in SEI models has thus far been essentially non-existent with only a few examples of ML potentials for the solid-water interface. [253][254][255] The large variations in inter-and intra-molecular forces only adds to the complexity of the challenge, which can be addressed going forward with a two-pronged approach. The data generation for model training needs to go beyond the often used method of using optimization and (or) MD trajectories, but actively searching and selectively including useful reactive structures intelligently with methods like, BUSS (bootstrapping uncertainty structure selection) [256] or meta-dynamics [257] without brute force sampling all possible higher energy parts of the potential energy surface, which is an impossible task for complex systems. Intelligent sampling of data has helped to build reactive and accurate ML models for simple reactions. [258] On the other hand, models need to be designed to be expressive toward both long-and short-range interaction, covalent and non-covalent interactions and capable of modelling field effects. The three classes of ML models are currently under development to provide such capabilities. For example, modifications that allow non-local charge transfer in fingerprint-based models [259] or graph neural network models [260] are important for accurate electrostatics.

Electron Density and Other Properties
Access to electron density can help us in creating a detailed understanding of SEI formation and dynamics-from bonding and interaction type [90,261] to assigning redox reactions on atoms via changes in oxidation states. [262,263] Charge density can be used as a descriptor to map Li ion migration paths in solid-state phases [264,265] or intercalation locations. [266] While ML potentials let us break the cost-accuracy trade-off for total energy and force estimations (and derived properties from those), doing the same for volumetric data like density distribution accurately is much more of a challenge. All the major groups of ML models popular in the ML potentials community have recently also been tried out to predict density, that is, Gaussian process, [267] neural networks with fingerprints [268] and graph neural networks ( Figure 5). [269] Predictions have been tried out with basis functions [270] but also on a grid. [271] Although most of the existing effort has been for simpler systems, our work has shown high predictive accuracy for disordered and partially lithiated NMC cathode materials and organic liquid electrolytes by using a probe-based message passing graph neural network. [271] This approach might be most suitable for complex systems like the SEI with many phases, chemistries and element types involved. Similarly other important observables at the electronic length scale, like the density of states, can also be evaluated with graph networks at very little cost. [272]

Reaction Network Graphs
The ability to predict the chemical reaction mechanisms and products of the SEI, based only on information about the molecules and materials that are formed during SEI reactions [68] is of importance. However, the immense complexity of the SEI reaction network has prevented such an effort till recently. Knowledge of the structures and free energies of all the energetically accessible reaction intermediates, transition states and products provides sufficient thermodynamic and kinetic information about the reaction network to make quantitative predictions. An automated workflow based on a systematic generation of relevant species using general bond fragmentation/recombination operations can be used to explore SEI reaction networks. [68,[273][274][275] However, the crucial obstacle is that vast parts of the potential intermediate/product phase space are typically unknown. It is challenging to obtain reliable energetics of all potential reaction intermediates and transition states with high accuracy DFT simulations. Recently, structural and thermodynamic properties have been calculated with DFT for over 17,000 unique electrolyte and interphase species. [276] While empirical correlations exist between thermodynamics and kinetics between competing reactions (Bell-Evans-Polanyi principle [123] ), for practically infinite reaction networks with many almost degenerate branches, like that related to SEI formation, can not be explored solely on thermodynamic numbers. Thus, thermodynamic properties of possible phases are only the tip of the iceberg, information about reaction path determining barriers is missing. The data volume originating from (reaction) network connections within this library of species is enormous. To map out large reaction networks, we need to build autonomous reaction explorer models, [277] but it requires access to simulators that provide accurate energetics of reaction intermediates and transition states. ML potentials can be used for this purpose as discussed in the previous subsection. However, such surrogates are expected to have varying levels of uncertainty in the predicted energies [278,279] requiring reaction network analysis with uncertain data [280] and ML models can also help in this aspect both for molecular [281] and solid-state reactions. [282]

Surrogate Models for Phase-Field Simulations
Beyond atomic-scale, access to high accuracy/resolution, 3D phase-field simulations considering many different SEI phases together is the long term goal toward a clear understanding of SEI formation and micro-structure evolution. At a higher scale, outcomes from such models will also help us go away from nonphysical homogeneous models of SEI that are used in stateof-the-art simulations (see Section 3.3.5). Due to the sequential nature and grid-based approach to solving concerned equations, high-resolution 3D phase field with many phases can not be performed in a viable time frame. [283] Additionally obtaining reliable parameters (like interface energies, diffusivity etc.) from high fidelity QM simulations can also become a bottleneck when many phases are involved. ML potentials can help obtain high-fidelity parameters without expensive QM simulations as described in Section 4.1, but it can play a pivotal role in making high-resolution 3D meso-scale phase-field simulations viable. ML can be used to either build surrogate models that forego the sequential simulation and predict long term equilibrium structures directly or ML models can replace the partial differential equation-based sequential solution by a ML-based one to accelerate the simulations. In the direct prediction regime, final outcome structures or feature based representations of them are generated based on the initial seed and input parameters. Only recently, recurrent neural networks like LSTM, [284,285] Bayesian optimization [286] and auto-regressive Gaussian processes [287] have been tried for this purpose. Sampling the phase space of the complex SEI micro-structure would need large datasets, which would make Bayesian models less cost-effective. We also envision image-to-image translation models like convolutional U-Net [288] or latent space-based generative models like variational auto encoders [289] that can be used for this purpose. Also convolution-based models are expected to be more efficient due to the ease of learning spatial correlations that are important in micro-structure evolution. [290] As autoregressive models are restricted to learn spatio-temporal correlations based on the datasets, it limits the predictive ability on a long time horizon. Thus, compared to end-to-end prediction, a more efficient strategy can be to use ML to work with Figure 5. Machine learning models for accelerating atomic-scale simulations need to represent atomic neighbourhoods as fingerprints. Atom-centered symmetry functions [184,248] and smooth overlap of atomic positions (SOAP) [246] are very popular fingerprints. Instead of designed fingerprints, graph neural networks can learn representations. [249,250] For fitting atomic structure representations to target properties like energy and forces [13] or electron density [271] or density of states, [272] one can use regular neural networks [186] or Gaussian processes [245] or graph neural networks, [249] among many others. Images related to "Atom centered symmetry functions": Adapted with permission. [346] Copyright 2018, AIP Publishing. Images related to "Smooth overlap of atomic positions (SOAP)": Adapted with permission. [347] Copyright 2020, the Royal Society of Chemistry. Images related to "Graph representation with node vectors" and "Graph neural networks": Reproduced with permission. [250] Copyright 2018, American Physical Society.
coarse time/length resolutions while maintaining accuracy. [290] Using ML to accelerate time evolution in phase-field models should be especially useful if the ML models work directly with micro-structure images, which are high dimensional and thus learning very long time scale correlations might not always be feasible.
Generation of final micro-structure directly is challenging as the inital seed is stochastic in phase-field simulations and thus a well-defined distribution of the output images can not be learnt as a function of it. However the structural characteristics of the output solution, for example, domain size and distribution, are fundamentally linked to the simulation parameters and can be learnt in a data-driven manner from a large number of phase-field simulations. Thus a low-dimensional representation of the micro-structure can be very effective in building lightweight ML models that forecast micro-structure evolution accurately in the latent space, [285] which later can be used for the reconstruction of the micro-structure. [291] Long time scale phenomena like SEI crack formation [44] typically are simulated with phase-field models considering the SEI to be a homogeneous phase. [292] Such simulations will be far more reliable if detailed information on the inhomogeneous structure of the SEI is available through high-resolution 3D models performed through ML surrogates or accelerators.

Scale Bridging with ML
The traditional multi-scale modelling approach [45,293] works by sequentially connecting smaller-scale models (that treat the system at finer details with corresponding physics equations) to macro-scale models (working at coarser spatial and temporal resolution) by means of parameter passing (results from finerscale models become the input for higher-scale models). This is a bottom-up approach. Scale bridging in the opposite direction is known as top-down approach. Due to the approximations involved in the models used at different scale levels, straightforward parameter passing propagates errors that tend to increase along the simulation chain. This drawback can be mitigated partly by opting for iterative coupling between scales, [294] but it involves significantly higher computational complexity. The need for expert intervention in performing such models within model uncertainty bounds makes automated workflow-based simulations difficult and limits its applicability in SEI simulations and predictive interphase design. For a specific set of interlinked models, error calibration can be done. Nevertheless, multi-scale models being fundamentally distinct from the typical single-scale models, the calibration of a multi-scale model is intricate and specific to a micro-model connection. [295] SEI modelling requires not just multi-scale but also multiparadigm simulations. [293] Thus necessary multi-method simulations and uncertainty handling need to be flexible toward the set of methods being bridged. Machine learning and data workflow automation can significantly improve the fidelity of multiscale modelling in multiple ways. [296] Development of Systematic Multi-Scaling Techniques by Optimizing and Automating Model Coupling: a small change in electrolyte composition requires re-parametrization of the multi-scale model of the SEI. Thus, the evaluation of new electrode/electrolyte combinations for a better SEI is a long, slow, and costly process. The need for automated scale bridging is a major challenge in SEI modelling. [297] Making the data transfer between different scales flexible through scale-bridging APIs and data ontologies like the battery interface ontology Bat-tINFO (https://github.com/BIG-MAP/BattINFO) developed in the BIG-MAP project, the scale bridging process can become autonomous and move away from the construction procedure by hand.
Reduction in Computational Cost and Uncertainty Propagation with ML Surrogates: we can adaptively replace specific models within a multi-scale modeling setup with ML-based surrogate models to reduce the computational cost. [298,299] Such a setup allows us to probe contributions of individual scale bridging parameters to error propagation in the model hierarchy through a local functional analysis (e.g., Taylor series). [300] For example in a bottom-up approach, global sensitivity analysis carried out for a phase-field model for uncertainty in output relative to uncertainty in the input parameter space (those obtained from, e.g., DFT simulations) indicates parameters (e.g., specific interface energy) that need higher-fidelity computation to help reduce the uncertainty of the phase-field model output. Furthermore, this scheme allows multi-scaling with multi-fidelity models [301] based on uncertainty bounds.
Automatic Coarse Graining with ML: Coarse graining enables simulations of longer time and length scales. Coarse-grained models use similar constitutive equations but finer details of the system description are overlooked. However, the coarse graining model must not change the important part of the simulation outcomes and be consistent with the conclusions drawn from a model at a finer level of detail. Empirical coarse-graining depends on an effective model definition approach to achieve this. [302] The assessment of parameters for coarse-grained simulations can be defined as a machine learning problem. It involves two coupled learning problems: defining the mapping from a finer representation to a reduced representation, and parameterizing a model defined for coarse-grained representation. A supervised ML [303] or generative modeling framework (variational auto-encoder) [304] can be used to unify the tasks of learning discrete coarse-grained variables, decoding back to finer details, and parameterizing coarse-grained models. Benefits of ML-based coarse graining at the atomic scale have been demonstrated but it holds good promise for the continuum scale as well, [305,306] the fundamental approach being similar.

Deep Descriptors from ML
The black-box nature of complex machine learning models stops us from obtaining scientific insights from them. To make ML models physics-driven, we also need to ensure the predictions from deep models come from the combinations of data features that are appropriate in the context of the real underlying physical phenomena. Methods for interpreting deep learning models are generally known as 'explainable AI (XAI)'. [307] Using XAI to understand how complex ML models work can provide new scientific insight and design principles. Machine learning models can also autonomously generate scientific hypotheses without human intervention. [308] We can use deep learning for understanding some of the most complex systems like battery interphases where descriptors might be in high dimensional space beyond the capacity of manual construction. With XAI, we can discover complex descriptors [309] going beyond simple linear descriptors. [102] Descriptors can help us to efficiently search for high performance materials by locating regions of interest in phase space with limited experimental/computational cost. [22,[310][311][312] For example, the charge density can be a descriptor to ionic migration paths in solid-state electrolytes, [265] or the electronegativity of atoms are descriptors to intercalation potentials. [313] Approximate property estimation through descriptors is very fast in screening target systems. [314] Data-driven descriptor discovery can be done with statistical/symbolic models or ML models. [315][316][317] Both of them have been successfully used in the materials science community and batteries as well. [22,318,319] In ML model development, we compromise between accuracy and explainability. Most accurate models (e.g., deep neural networks) are the least explainable, and the most explainable models (e.g., linear regression) have low accuracy. Structure-property relations portrayed with machine-learned descriptors can be verified with high-fidelity computational or experimental data in a similar way human scientific hypotheses are tested. For the battery interfaces, comprehensible models with reasonable accuracy have been adopted. [320,321] However, there is a clear trend in deep learning models quickly gaining popularity in this application space opening up the possibility of applying XAI methods for making those models interpretable.

Inverse Generative Design
Deep generative models like variational autoencoders (VAE) [322] or generative adversarial networks (GAN) [323] can learn chemical constraints in electrode/electrolyte composition and structures at the atomic scale [327] as well as the continuum scale [328] from large datasets (Figure 6). These models can come up with new electrode/electrolyte combinations that are novel but chemically viable. Beyond the chemical rules learnt from data, prior expert knowledge can be utilized for the generation of practical battery systems. [329,330] Although most of the current application of generative models at the atomistic scale focused on small molecules, few practical cases of solid state materials have been demonstrated. [331][332][333] At the microstructure scale, generative models for solid materials have been very successful even for complex heterogeneous systems and limited datasets. [325,[328][329][330][334][335][336][337] Generative models can be used for inverting the design process of better battery interphases by going from properties to structure-and composition-conditional generative models that learn correlations in structural space conditioned to the properties. [15] Thus, favorable SEI forming electrolyte compositions can be probabilistically generated. With appropriate training structure-property datasets, conditional VAE (cVAE) [326,338,339] and conditional GAN (cGAN) [324,340] can be trained. Successes in conditional generative models for molecules and bulk solids have been reported recently. Batteries are interfacial systems incorporating numerous condensed-matter phases. Figure 6. Generative models like variational autoencoders [322] and generative advarserial networks [323] can be used for learning patterns in data like atomic structure graphs and 3D micro-structure. Trained models are used to generate new physics/chemically valid structures that were never part of the existing datasets. When such generative models are trained to generate structures conditioned to properties, [324][325][326] we can perform inverse design of interfaces that possess properties aimed at high-performance safer batteries. The SEI being heterogeneous in nature, generative models dealing with the SEI must be able to work with multiple data types-from small graphs of electrolyte molecules, periodic crystal graphs of solid phases and 2D/3D images of meso/ micro-structures.
High performance and safe battery interphases also require optimization for multiple target properties. Thus, such conditional composition generation needs to be done for multiple properties at once which is very challenging. The observed properties are controlled by the structure at multiple length scales, [341,342] and their evolution over time as well. It is very data-intensive and requires datasets across multiple length and time scales (BIG-MAP data paper). We will need new types of generative models that can work with solid/liquid/ organic/inorganic phases and different forms of data representation like graphs and voxels. For doing so, with a multiresolution and multimodal generative model, we can extract representations for each data stream and fuse features using crossmodal multihead attention layers and train in a multitask setting. [343] To represent data at multiple scales, a hierarchical VAE [344,345] that works with multiple latent spaces for multiple length scales that are conditioned to each other [37] can be used ( Figure 6). As it is trained for structural correlations across scales it can be used to perform scale bridging through oneshot generation of multiscale structures without the need of step by step parametric multi-scaling.

Outlook
While classical, electronic structure and continuum level modelling have unquestionably provided valuable insights into the nucleation and compositional nature of the SEI over the last two-to-three decades, [5] it is equally evident that none of these traditional modelling scales or multi-scale approaches bridging them have been able to capture the fundamental physico-chemical complexity of the spatio-temporal evolution of the SEI. The advent of machine and deep learning techniques in the areas of batteries and electrochemical interfaces have added a new dimension to the quest for predictive SEI models, as outlined above. It is, however, important to stress that such data-driven approaches must be both physicsinformed and uncertainty-aware [279] in order to ensure predictive accuracy for spatio-temporal SEI models. Advanced multi-scaling approaches in latent space, for example, using hierarchical latent space models, also provide new possibilities for parameter-free multi scaling and to identify deep descriptors for SEI evolution, as does explainable AI. Given the complexity of the SEI, the need for access to suitable training data is enormous and autonomous workflows and automated approaches to ingest multi-sourced and multifidelity pose a critical bottleneck, which must also be solved in order to harvest the full potential of such hybrid physicsinformed and uncertainty-aware data-driven SEI models. Youssef Mabrouk pursued his undergraduate and graduate studies at the Department of Physics of the University of Munich, where he focused on condensed matter physics and statistical physics. During his graduate studies he was working as an assistant at the Institute of Materials Physics in Space at the DLR, where he wrote his thesis about the data-driven discovery of governing equations. Mr. Mabrouk is now pursuing his doctoral studies at the Helmholtz-Institute Münster, where he is focusing on atomic-scale modelling of battery electrolytes and interfaces. The use of data-driven techniques to facilitate collaboration between theory and experiment is also a focus of his project. His group builds and applies methods for data driven models for multiscale simulations of energy materials. He is co-PI in several collaborative research projects focused on machine learning-based multi time/length scale inverse design and WP leader of the "Battery Interface Genome" WP of the BIG-MAP project.