Investigation of the Impact of High Concentration LiTFSI Electrolytes on Silicon Anodes with Reactive Force Field Simulations

The initial formation cycles are critical to the performance of a lithium-ion battery (LIB), particularly in the case of silicon anodes, where the high surface area and extreme volume expansion during cycling make silicon susceptible to detrimental side reactions with the electrolyte. The solid electrolyte interface (SEI) that is formed during these initial cycles serves to protect the surface of the anode from a continued reaction with the electrolyte, and its composition reflects the composition of the electrolyte. In this work, ReaxFF reactive force field simulations were used to investigate the interactions between ether-based electrolytes with high LiTFSI salt concentrations (up to 4 mol/L) and a silicon oxide surface. The simulation investigations were verified with galvanostatic testing and post-mortem X-ray photoelectron spectroscopy, revealing that highly concentrated electrolytes resulted in the faster formation and SEIs containing more inorganic and silicon species. This study emphasizes the importance of understanding the link between electrolyte composition and SEI formation. This ReaxFF approach demonstrates an accessible way to tune electrolyte compositions for optimized performance without costly, time-consuming experimentation.


Introduction
Due to the increase in lithium-ion battery (LIB) usage and production, it is essential to explore options to produce energy storage systems with high energy densities that are based on cheap, earth-abundant, non-toxic materials. Full-cell LIBs with sulfur cathodes and silicon anodes have generated a lot of attention in recent years due to the high specific capacity of both materials (~1600 mAh g −1 and 4200 mAh g −1 , respectively) [1][2][3]. In order to combine these electrodes into a full cell, it is important to develop electrolytes that are compatible with both of them. The carbonate solvents, commonly used with silicon anodes, undergo detrimental reactions in the presence of lithium polysulfides [4]; therefore, sulfur-compatible electrolytes must be developed for silicon anodes.
With the growing use of LIBs in vehicles or other applications, the safety of the batteries is coming into question. In particular, the electrolyte solution is composed of volatile solvents that are flammable or can break down into toxic by-products. Electrolytes with lithium salts at high concentrations have been shown to have much lower volatility and better thermal stability compared to the conventional, more dilute electrolytes and to be Liquids 2023, 3 133 less likely to produce flammable gases from solvent evaporation or thermal decomposition reactions [5][6][7]. During the initial formation cycles of a battery, the electrolyte undergoes decomposition reactions to form the solid electrolyte interface (SEI) on the surface of the electrode [8,9]. These formation cycles are critical to the long-term function of the cell, particularly in the case of silicon, where the electrode/electrolyte interface is more susceptible to reactions due to the high surface area and expansion of silicon during cycling [10][11][12][13]. The reduction in organic solvent molecules can lead to polymerization reactions, forming a thick, ionically insulating SEI, which increases cell impedance and is detrimental to cell performance [14]. Highly concentrated electrolytes have been shown to result in a more inorganic SEI layer, which would lead to better mechanical and chemical stability of the SEI as well as improved Li + ion conductivity [9,15]. In addition, high salt concentrations have been shown to improve the electrochemical performance of other highenergy-density electrodes. In the case of sulfur cathodes, high lithium salt concentrations have been observed to minimize the polysulfide shuttle effect and to improve capacity retention [16,17]. In the case of lithium metal anodes, high salt concentrations have been observed to result in the non-dendritic reversible plating and stripping of lithium during cycling [18][19][20]. For applications with silicon anodes, it is important to better understand the compatibility of the silicon active material with highly concentrated lithium salt electrolytes in these early formation cycles.
In order to simulate the processes that occur during SEI formation, different theoretical methods can be used and were indeed used in previous studies. The most straightforward way is to simulate a representative cutout of the electrode-electrolyte interface by means of density functional theory (DFT). A hybrid functional, such as B3LYP [21,22] or PBE0 [23], should be able to describe all reactive processes during SEI formation with reasonable accuracy. The problem is, however, that DFT simulations are much too expensive to be performed on a system with several thousands of atoms, which must be included in a realistic SEI cutout over at least the ns time scales needed for the SEI buildup to occur. When SEI formation studies were performed based on DFT, they concentrated on the calculation of intermediates and reaction paths [24][25][26], thus neglecting the dynamical interplay of different reactive processes.
Large length and time scales, on the other hand, can be treated by traditional nonreactive force fields (FF) such as OPLS [27,28]. By avoiding the direct description of electrons and transferring their physics to the initial optimization of force field parameters, the energy and gradient calculations are accelerated by between four and six orders of magnitude in comparison to DFT. The scaling is also heavily reduced, being quadratic to linear with respect to the number of atoms in the large system limit, depending on the treatment of non-covalent interactions. The drawback, however, is that traditional chemistry FFs are incapable of describing bond breakings and formations, i.e., the largest part of chemical reactions. A bond matrix must be defined at the beginning of the calculation and cannot be updated dynamically. This works fine for the time-dependent description of electrolyte solvents since no reactions usually occur in them. When turning to SEI formations, however, reactive processes are crucial, and traditional chemistry FFs are incapable of describing them.
The ReaxFF reactive force field [29] bridges the gap between DFT and nonreactive FFs. By replacing the rigid bond matrix through bond orders dependent on the actual distance of the two atoms, ReaxFF can describe bond fissions and formations. The effort to calculate a ReaxFF gradient is roughly two orders of magnitude higher than with nonreactive FFs, even with quadratic or even sub-quadratic scaling, so that large systems containing tens of thousands of atoms can also be simulated over the range of ns. The central challenge in running ReaxFF simulations is the parametrization of the force field. Its broader versatility in comparison to nonreactive FFs comes hand in hand with a much more sophisticated formulation of the potential energy function, which needs a large number of element-and element-multiplet-specific parameters to be set up. Additionally, parameter interdependencies are more frequent and subtler than in traditional chemistry FFs.
If a new area of application is treated, a considerable fraction of those parameters must be re-optimized, which is far from trivial both in terms of training-set selection and the finding of a global optimum.
Despite those challenges in its setup, ReaxFF has proven to be of high use in many application areas, and hundreds of publications cover both liquid chemistry and material science, i.e., both sides of the SEI [30,31]. When setting the focus on silicon-related studies, the catalysis of the methanol to olefin reaction on a SiO x /AlO x zeolite was simulated with ReaxFF [32], as well as the processes that occurred at the silica-water interface [33], the interaction of SiO 2 with molecular oxygen during surface oxidation [34], the reaction of alkyl radicals on a Si(111) surface [35] and tribochemical reactions of phosphoric acid molecules sheared between two SiO 2 surfaces [36].
LIBs in particular is indeed one of the main application fields of ReaxFF, with studies published covering the structure and dynamics of all LIB parts, namely cathodes (sulfur, sulfur with carbon nanotubes/nanofibers), anodes (carbon, silicon, TiO 2 , Li-metal) and electrolytes (EC, DMC, LiPF 6 , LiTFSI) [37]. In addition, ReaxFF parameters of Li-atoms for usage within organolithium compounds and small Li-clusters with hydrogen and carbon were optimized [38], and the intercalation of lithium into graphene and Si anodes was simulated as well [39][40][41].
Despite numerous applications in the LIB simulation, studies describing the buildup of the SEI with ReaxFF are still scarce. A good example is the study by Kim et al., who simulated the interface between a lithium metal anode and an electrolyte containing EC and DMC [42]. No studies, however, have been published so far that handle the much more complex SEI formed by electrolytes containing ionic liquid counter ions, such as TFSI, which is probably due to the need for ReaxFF to be able to describe eight elements and their mutual interactions at once. This is clearly in the upper range of complexity with respect to the average number of elements simulated in ReaxFF studies.
In addition to ReaxFF, the REACTER algorithm [43] also has the ability to combine the speed of force fields with the quantitative description of chemical reactions occurring using DFT. In REACTER, the system is described as a traditional, nonreactive force field; reaction events are recognized automatically based on the criteria predefined by the user. During a reaction event, the force field description (bond matrix, parameters, and partial charges) switches to the product configuration, whereas several criteria (pre-calculated activation barriers, configuration matching, etc.) can be evaluated to fine-tune the reactivity. Although this approach is also applied to the growth of SEIs in Li-ion batteries [44], it requires a larger amount of user interference in its setup since all possible reactions must be predefined in all relevant details, whereas in ReaxFF, the underlying potential energy description allows for a certain amount of extrapolation in this regard, based on the reference set used. On the other hand, REACTER can be run with essentially the speed of nonreactive force fields, which makes it much faster than ReaxFF. It might therefore be used to overcome the current limitations with respect to the simulation time in our future studies.
The present publication outlines a framework to set up a suitably parametrized ReaxFF force field, running the simulations on an appropriate model and evaluating the results for the extraction of general trends in the case of SEI formation; a framework that can easily be generalized to similar studies in this topic will be considered in future work. To verify the SEI formation mechanism on the silicon active material surface for highly concentrated lithium salt electrolytes, pure porous silicon anodes were used. Using pure silicon anodes ensures that all of the organic components observed with X-ray photoelectron spectroscopy (XPS) come from the SEI. Structured pure silicon anodes were chosen because they have already shown good long-term cycling stability [45][46][47][48][49] and an understanding of SEI composition for these anodes, which can additionally be applied to silicon anodes containing inactive components such as conductive additives and binders. Keeping in mind the desire to use silicon anodes in combination with sulfur cathodes, solvents that are compatible with sulfur cathodes were also chosen, i.e., ether solvents. In particular, the wellknown monoglyme, dimethoxy ethane (DME), together with dioxolane (DOL), were used. Lithium bis(trifluoromethanesulfonyl)imide (LiTFSI) was chosen as the salt since lithium hexafluorophosphide, LiPF 6 , the most popular Li-salt for LIBs [50], reacts detrimentally with ether solvents [51]. The low association between Li + and the [TFSI] − counter ion means that the ionic conductivity would be less impacted at higher concentrations compared to lithium salts, with a stronger association between the Li + ion and anion [52].
For experimental verification, the electrochemical performance of the silicon anode in the early formation cycles in electrolytes with low and high Li-salt concentrations is investigated with galvanostatic measurements. The composition of the SEI after the initial formation cycle is investigated in more detail using X-ray photoelectron spectroscopy (XPS). To gain a better understanding of how electrolytes with high salt concentrations can impact the cycling behavior of pure silicon anodes, the electrochemical performance of the porous silicon anodes during these cycles is linked to the electrolyte-surface interactions determined from the ReaxFF simulations, as well as the resulting composition of the SEI.

The ReaxFF Reactive Force Field
The mathematical structure of ReaxFF is fundamentally similar to that of a traditional nonreactive force field, i.e., the energy is calculated solely based on the positions of the atomic cores, and electrons are only treated indirectly by means of damping parameters and correction terms [53].
The force field is adapted to the system of interest by a large set of parameters that affect every atom, bond, angle, and other terms.
The central difference between ReaxFF and nonreactive force fields such as OPLS is that the list of bonds in the system is not fixed but dynamic. Which bonds exist between which atoms are not determined by a fixed, user-given bond matrix but solely by the actual geometry of the atoms that change during the molecular dynamics. Depending on the actual cartesian distance r ij , a bond order (BO) can be calculated between any pair of atoms i and j with: where r σ 0 , r π 0 , and r ππ 0 are the equilibrium distances of sigma, pi, and double pi bonds for the atom pair and p bo1 to p bo6 are atom-pair-specific parameters.
This concept makes it possible to describe chemical reactions, which typically consist of the breakage and formation of chemical bonds. One drawback is that nine parameters are needed for the bond order calculation between a certain combination of elements, in comparison to two or three parameters for the full bond energy in nonreactive FFs. Further, several correction terms taking the total number of bonds to an atom, or the number of lone electron pairs (being just a parameter here) into account, must be included for a reasonable description of the bonding energy. Similar approaches must be undertaken for valence angles and torsions. This results in the total force field (FF) energy taking the form of a sum of several contributions: Besides the usual nonreactive FF terms, including bonding (E bond ), angle (E val ), torsion (E tors ), dispersion (E vdW ), and Coulomb (E Coulomb ), several ReaxFF-specific terms appear as well, such as under-and over-coordination corrections, if the current atom has too many or too few bonding partners (E under , E over ), the terms considering parametric electron pairs (E e−pair ), hydrogen bonds (E H−bond ), conjugated systems (E conj ), and so on.
Coulombic interactions are crucial for the description of long-range interactions and are especially complicated to treat since no fixed charges can be predefined and no simple formula for the change in partial charges during chemical reactions can be given. Therefore, the charges are updated by the iterative electronegativity equilibration method (EEM) [54] in each energy calculation before being plugged into a slightly modified Coulomb term: where q i and q j are the partial charges of the atoms, C is a scaling constant and γ ij is a damping parameter adjusted to the overlap between atoms at close distances.

Evolutionary Algorithms
The optimization of ReaxFF parameters imposes a highly nontrivial task due to their large number (2298 for this study, if all are updated from scratch) and, due to their interdependency, the validity of bond energy parameters depends on atom-specific parameters and those for different correction schemes. Further, most of them have no direct physical meaning (in contrast to nonreactive FFs), which makes it much harder to evaluate if a certain optimized set is reasonable. An elegant and sufficiently powerful way to perform this challenging task is by applying evolutionary algorithms (EA).
In EAs, the nomenclature of the steps within the optimization process is inspired by natural selection and evolution [55,56], although this is a local competition and not global optimization. The practical efficiency of EAs for global optimization, however, derives from their ability to exploit the partial separability of the optimization problem under study. In the beginning, a population of trial force fields, with randomly (within limits) initiated force field parameters, is set up (also called the genetic pool). The fitness (or error sum) of each individual in this pool is defined as: where the sum runs over all n items of the reference data (energies, gradient components, charges, etc.), and the squared differences between the ReaxFF value and a reference (DFT in this study) are added and weighted by a factor σ i , which also takes different units of the reference data into account. The aim is, of course, to decrease the fitness as far as possible, which increases the agreement between ReaxFF data and reference data.
During the optimization process, crossovers between two individuals can occur that exchange blocks of parameters between them, leading to new individuals generated from them. This step exploits partial separability and allows for large jumps in the search space that have a good chance of leading to improvements. Furthermore, mutations can alter the parameter values of single individuals. This corresponds to a test for small-distance improvements in the search space. With a combination of random choice and preference for individuals with a lower error sum and proper balance between the long-distance crossover and short-distance mutation, the EA optimization samples the objective function surface effectively and converges to the global minimum. The convergence towards near-global optima is reached much faster in EAs than with, for example, completely randomized multistart local-search methods, but also, no complete sampling of the search space is needed (as with deterministic methods), which is impossible in the high-dimensional objective function surfaces present in the ReaxFF parameter optimization. With several additional tricks and settings, such as niching (for conserving genetic diversity within the pool) or different ratios of various mutation and crossover operators, the effectiveness of EAs can be raised even more. Further fundamental aspects of EAs are explained elsewhere [57,58].

Force Field Parametrization
Although no published ReaxFF parametrization covers the full system treated in this study, there are at least some publications that treat parts of it. In the originally intended building-block style of ReaxFF parameters, those were used as a starting point for the FF parametrization. The ReaxFF by Liu and Yu [59] covers the full description of LiTFSI and DME/DOL in the interface with lithium metal; the missing interactions are those between the electrolyte and the SiO 2 anode. Some of the required parameters were found in the ReaxFFs by Yun and Pai [60] and Wang and Shi [61]. Even after combining those three force fields, however, there were some missing parameters. These were essentially those for the S-Si bond and various angles containing this bond. To the best of our knowledge, no ReaxFF parameters for these interactions have been published so far; this publication, therefore, is the first one to do so. In addition to these parameters, some other bond terms that were thought to be important for the SEI description but were thus far missing in all three ReaxFF parametrizations were also added.
The coordinates, for which new ReaxFF parameters were optimized, are listed in Table 1.

S-Si
Si The optimization requires a suitable training set of the reference data, produced with a quantum chemical method that can be reproduced as well as possible by the optimized ReaxFF. Considering the reference used for the three publications, from which other parameters for the new ReaxFF were taken [59][60][61], B3LYP/6-311+G** [21,22,62,63] was chosen, as it is the same with which essentially all other reference data were calculated.
The training set was chosen to be a collection of small molecules with up to 21 atoms, where all of them contain at least one of the coordinates listed in Table 1 (Supplementary Material Figure S1). With this collection, two sets of reference data were obtained: first, optimized B3LYP geometries of the molecules, and second, relative B3LYP energies of the bond and angle scans of the newly added coordinates. Due to the high importance of the S-Si bond potential, four bond scans in different molecules were added as a reference for it, in contrast, to 1-2 scans for the angles. The chosen scan coordinates are also shown in the Supplementary Material ( Figure S2).
With the training set at hand, all new force field parameters were optimized simultaneously with the OGOLEM software package, which has the ability to use a wide variety of EA variants and options for the fitting of ReaxFF parameters [56,64]. Only the parameters of the new coordinates were optimized; all the remaining ones were kept at the values optimized by the original authors. This was assumed to be reasonable since the systems treated in the respective publications are similar to those treated here, and a complete re-optimization would have needed a training set spanning all possible intermediates and products within the SEI buildup. With this choice, the search space still had a dimension of 122 parameters. The allowed parameter ranges within the EA were chosen to be in the same order of magnitude as the parameters of similar coordinates, with a variance of ±100%. The global optimization operator was chosen to be a mixture of 10% Arctic and 90% Portugal operators; parent parameter vectors were chosen based on their fitness rank, with a Gaussian width of 0.2 and a mutation ratio of 0.1. The total population of the pool was 400, with 500,000 EA iterations having been conducted. We refer to the online OGOLEM manual [65] for further explanations of these EA setting choices.

Molecular Dynamics Simulations
In order to study the SEI buildup, the interface between the electrolyte and anode was simulated with reactive molecular dynamics (MD) using the newly parametrized ReaxFF. For all calculations, the LAMMPS software package was used [66][67][68][69]. A Nose-Hoover thermostat was used to simulate a canonical ensemble at the desired temperatures in all cases [70]. An amorphous SiO 2 surface in contact with two different electrolyte solutions, 0.7 M TFSI and 4 M TFSI, within a 2:1 (v/v) mixture of DME and DOL, should be simulated. In the first place, both the anode and the electrolyte part were generated separately. For the electrolyte solutions, the following procedure was applied: single boxes with pure DOL and DME, each consisting of 512 molecules, were set up and equilibrated with an NpT ensemble using the Nose-Hoover barostat at 1 bar and 298 K. The averaged relative volume of the boxes resulted in a translation of the 1:2 (v/v) ratio into a 1:139544 (mol/mol) ratio of DOL and DME. Based on this, a translation of the macroscopic LiTFSI concentration to molecule numbers was fitted. Details concerning this procedure are explained in the Supplementary Material (see Table S1 and Equation (1)).
After determining the correct concentrations, boxes elongated into z-direction were set up, which contained 384 DOL, 536 DME, and 60 LiTFSI molecules for the 0.7 M LiTFSI and 228 DOL, 324 DME and 428 LiTFSI molecules for the 4 M LiTFSI system.
The amorphous SiO 2 bulk was set up by placing Si and O atoms randomly onto a silicon diamond cubic crystal structure, where the first benchmark that ran suggested that the exact distribution of atoms within the bulk had no significant influence on the initial reaction dynamics. Both the electrolyte boxes and the SiO 2 bulk were equilibrated by running 1 ns NpT dynamics at 1 bar at 298 K, with a time step of 1 fs. Finally, both parts were brought together, where the x and y dimensions of the electrolyte part were held fixed during equilibration to the final dimensions of the anode part such that both sides of the interface had the same area. After removing the z-periodicity of both parts, electrolyte molecules with dangling bonds were plugged together on one side, whereas dangling bonds in the SiO 2 part were kept in place in order to accelerate the SEI growth. This modification serves to indirectly include the influence of an applied electric current on the SEI building process, which occurs much slower (on timescales too long to converge within reasonable ReaxFF simulations, as the benchmarking indicated) if a saturated anode surface is set up.
Starting from this structure, 1.4 ns of reactive MD was propagated for the 0.7 M and 4 M systems at 298 K and 700 K in an NpT ensemble at 1 bar, respectively, with a time step of 0.5 fs. The simulations at different temperatures were conducted in order to check if the SEI buildup could be accelerated without changing its mechanism significantly. Since such changes were indeed observed, all evaluations within the main manuscript are restricted to the 0.7 M and 4 M systems at 298 K, whereas the 700 K results are only shown in the Supplemental Information as a further benchmark regarding the limitation of the simulation method used.

Preparation of Porous Silicon Anode
In order to investigate the SEI composition on the silicon active material surface, pure silicon anodes were chosen. The anodes were structured to enable more stable cycling. The preparation of mesoporous silicon anodes, using a highly p-doped (80 mΩcm) silicon (100) substrate material, was performed with a large area electrochemical etching process, utilizing an "in-line" etching tool. Electrochemical etching was performed using a topdown approach, with an applied current of 13 mA cm −2 and an etching time of 14 min. The resulting mesoporous silicon layer possessed porosities of 43% and a thickness of 6 µm.

Electrochemical Characterization
Two different electrolyte solutions were prepared in a glovebox under an Ar atmosphere with both a low and a high-concentration electrolyte. Both electrolyte solutions were prepared with monoglyme (DME) and DOL in a 2:1 (v/v) ratio of DME:DOL. LiTFSI was used as the lithium salt in the electrolyte, and the concentrations were varied to make 0.7 M and 4 M electrolytes. All chemicals were obtained from Sigma Aldrich. Prior to their characterization, liquid electrolytes were sealed in polyethylene pouch cells in the glove box under an Ar atmosphere. Raman spectroscopy analysis was performed using an alpha 300 RA set-up. An excitation wavelength of 532 nm was used.
Electrochemical characterization of the porous silicon anodes was performed using a half-cell set-up with metallic lithium as the counter and reference electrode. The tests were performed in home-built inert electrode housings using Celgard separators. The electrolyte solutions with high and low lithium salt concentrations, as described previously, were used. A constant current-constant voltage (CC-CV) cycling procedure was carried out on an Astrol battery cycler (Astrol Electronic AG) with a capacity limitation of 75% of the theoretical capacity of silicon: 3150 mAh g −1 . A voltage range of 0.1 V-0.8 V vs. Li/Li + at 0.1 C was used. As this particular study focuses on the SEI composition from the initial formation cycles, only ten cycles were used for the electrochemical verification. This number was chosen to ensure that complete formation was achieved and that a stable coulombic efficiency was reached. Additionally, cyclic voltammetry was performed to investigate the degradation of the electrolyte on the surface of the porous silicon anodes in the initial formation cycle. The voltage was first decreased from the open circuit potential to 1.9 V, with a rate of 0.28 mV s −1 , after which the voltage was further decreased to 0.1 V, with a rate of 0.0194 mV s −1 (corresponding to a rate of 0.1 C).

Postmortem Characterization
The ex-situ postmortem characterization of the chemical composition of the SEI, formed in the first cycle, was performed using X-ray photoelectron spectroscopy (XPS). Prior to characterization, the anodes were rinsed in an Ar-filled glove box with a solution of DME/DOL (2:1 (v/v)) to ensure the removal of all soluble SEI components. The samples were then vacuum dried for 5 min. XPS measurements were conducted using an instrument from Omicron Nano-Technology, GmbH, Taunusstein, Germany, with an Al-anode source (240 W). The spectral analysis was performed using CasaXPS software, and the spectra were referenced to aliphatic carbon at 285 eV.

Results
The ReaxFF simulations were successfully implemented to investigate the interaction between high-concentration, ether-based electrolytes and a silicon oxide surface (representative of the silicon anode surface chemistry). Experimental verification was later performed by investigating the SEI composition after the first formation cycle using X-ray photoelectron spectroscopy (XPS). The SEI composition after ten cycles was also investigated to determine how the SEI was developed after the initial formation cycles. The results are displayed and discussed below.

Force Field Parametrization
The results of the ReaxFF parameter optimization conducted for the description of the SEI buildup are presented first. In Figure 1, the total error sum is plotted against the EA iteration number. The convergence was reached quite fast, at a fitness of approx. 200,000, and the final value of the optimization is 185,288.6. The fast convergence and small deviations between the results of different EA runs (not shown here) indicate an easy EA problem, which is The convergence was reached quite fast, at a fitness of approx. 200,000, and the final value of the optimization is 185,288.6. The fast convergence and small deviations between the results of different EA runs (not shown here) indicate an easy EA problem, which is not surprising since a number of mostly independent parameter sets belonging to different coordinates were optimized. The search space is then effectively decomposed into several subspaces, with each dimension from 7 to 15 (depending on the coordinate type). Such low-dimensional problems are trivial to solve for OGOLEM (but are still hard/impossible to handle for other probabilistic or deterministic methods). More significant is how well the optimized ReaxFF reproduces the data collected in the training set. This is shown for the geometry part in Figure 2. Most bond lengths are reproduced well (less than 0.1 Å deviation), building four clusters at 0.95, 1.1, 1.5, and 1.7 Å, while a small fraction showed larger deviations, up to 0.9 Å. This indicates, however, that all the structures were stable, and no dissociations occurred. For bond angles, the average deviations were between 10 and 15 degrees, but no outliers are present here, in contrast to the bond lengths. Therefore, angles are also reproduced reasonably. The dihedral angles are reproduced with good quality, showing somewhat smaller deviations than bond angles, around 10 degrees. In summary, the reproduction of the geometries is reasonable, albeit not of outstanding quality. The convergence was reached quite fast, at a fitness of approx. 200,000, and the final value of the optimization is 185,288.6. The fast convergence and small deviations between the results of different EA runs (not shown here) indicate an easy EA problem, which is not surprising since a number of mostly independent parameter sets belonging to different coordinates were optimized. The search space is then effectively decomposed into several subspaces, with each dimension from 7 to 15 (depending on the coordinate type). Such low-dimensional problems are trivial to solve for OGOLEM (but are still hard/impossible to handle for other probabilistic or deterministic methods). More significant is how well the optimized ReaxFF reproduces the data collected in the training set. This is shown for the geometry part in Figure 2. Most bond lengths are reproduced well (less than 0.1 Å deviation), building four clusters at 0.95, 1.1, 1.5, and 1.7 Å, while a small fraction showed larger deviations, up to 0.9 Å. This indicates, however, that all the structures were stable, and no dissociations occurred. For bond angles, the average deviations were between 10 and 15 degrees, but no outliers are present here, in contrast to the bond lengths. Therefore, angles are also reproduced reasonably. The dihedral angles are reproduced with good quality, showing somewhat smaller deviations than bond angles, around 10 degrees. In summary, the reproduction of the geometries is reasonable, albeit not of outstanding quality.   Figures S4-S6). The S-Si bond length profile (a) shows a qualitative agreement, although with different shapes. The ReaxFF curves have a much narrower minimum and reach their dissociation energy at approx. 3.5 Å, whereas the DFT curves converge at a much slower rate. The characteristic parameters, however, are reproduced satisfactorily. The ReaxFF equilibrium length is somewhat larger than the DFT references, but the absolute deviation is only around 0.1 Å. The dissociation energies deviate between 50 and 60 kJ mol −1 . The angle scans show quite a variety of behaviors. ReaxFF shows very satisfying qualities for some of them (b), whereas other profiles show large deviations in the energetics (c). Most profiles, however, reproduce a correct equilibrium angle and do not diverge.
The reproduction of reference data seems not to be of outstanding quality. This, however, could be expected in advance: (1) ReaxFF, especially when parameterized for such a large collection of elements as conducted here, usually does not have the ability to reproduce whole potential energy surfaces with chemical accuracy. The aim is rather to obtain the correct relative activation barriers and structural parameters in order to make fast and qualitative simulations of complex reaction networks and alternative reaction paths accessible.
In fact, the level of agreement achieved here is on par with typical ReaxFF fittings published in the literature. (2) Only a small subset of the important parameters was provided for optimization.
All element-specific parameters and all bond parameters besides those for the S-Si bond were kept fixed. The reason for this is that those global parameters have a strong influence on any chemical situation that might occur during the simulation.
Much better results regarding the present training set could be achieved if more parameters are released for optimization. Then, however, the large landscape of chemistry considering the reactions of DOL with TFSI, and products containing F, S, and N atoms, etc., would have been altered as well. A huge training set covering the information and references to all these possible reactions must be provided if this should be conducted.
Liquids 2023, 2, FOR PEER REVIEW 10 The S-Si bond length profile (a) shows a qualitative agreement, although with different shapes. The ReaxFF curves have a much narrower minimum and reach their dissociation energy at approx. 3.5 Å, whereas the DFT curves converge at a much slower rate. The characteristic parameters, however, are reproduced satisfactorily. The ReaxFF equilibrium length is somewhat larger than the DFT references, but the absolute deviation is only around 0.1 Å. The dissociation energies deviate between 50 and 60 kJ mol −1 . The angle scans show quite a variety of behaviors. ReaxFF shows very satisfying qualities for some of them (b), whereas other profiles show large deviations in the energetics (c). Most profiles, however, reproduce a correct equilibrium angle and do not diverge. The reproduction of reference data seems not to be of outstanding quality. This, however, could be expected in advance: 1) ReaxFF, especially when parameterized for such a large collection of elements as conducted here, usually does not have the ability to reproduce whole potential energy surfaces with chemical accuracy. The aim is rather to obtain the correct relative activation barriers and structural parameters in order to make fast and qualitative simulations of complex reaction networks and alternative reaction paths accessible. In fact, the level of agreement achieved here is on par with typical ReaxFF fittings published in the literature. 2) Only a small subset of the important parameters was provided for optimization. All element-specific parameters and all bond parameters besides those for the S-Si bond were kept fixed. The reason for this is that those global parameters have a strong influence on any chemical situation that might occur during the simulation. Much better results regarding the present training set could be achieved if more parameters are released for optimization. Then, however, the large landscape of chemistry considering the reactions of DOL with TFSI, and products containing F, S, and N atoms, etc., would have been altered as well. A huge training set covering the information and references to all these possible reactions must be provided if this should be conducted.
Since most parameters were taken from three papers covering interface reactions with TFSI and DOL/DME, all this chemistry is already included, and there is no need and no rational cause to overwrite all this. Such a global re-optimization could, however, be performed if electrolytes/electrodes with a totally new composition could be simulated.
Nevertheless, the optimized force field can be expected to lead to correct simulation results regarding the SEI composition since this largely depends on the relative size of the reaction barriers. Since most parameters were taken from three papers covering interface reactions with TFSI and DOL/DME, all this chemistry is already included, and there is no need and no rational cause to overwrite all this. Such a global re-optimization could, however, be performed if electrolytes/electrodes with a totally new composition could be simulated.
Nevertheless, the optimized force field can be expected to lead to correct simulation results regarding the SEI composition since this largely depends on the relative size of the reaction barriers.

Element Distributions
In the next sections, the results obtained from the ReaxFF sampling trajectories will be shown. In Figure 4, two screenshots of the model battery systems with 0.7 and 4 M TFSI within the electrolyte are shown during SEI formation. Further screenshots can be seen in the Supplementary Material ( Figures S7-S9). Obviously, the simulated system is too large and too complex to draw conclusions from such visual inspections alone.
The simplest way to analyze the results of the reactive dynamics is to look at how the eight different elements are distributed in space. Since we are interested mainly in the direction perpendicular to the SEI interface, i.e., the z-axis, the evaluation program packed all atoms of the respective elements into the system with 200 bins along the z-axis and separately for each written time step (one per ps, i.e., 1400 frames in total). The results are shown in Figure 5 for selected elements (O, F, and Li). Complete results, as well as simulations at different pressures and temperatures, are shown in the Supplementary Material (Figures S10-S14). Since the system is periodic in the z-direction, essentially, two SEIs can be seen at z ≈ −18 Å and z ≈ 50 Å.

Element Distributions
In the next sections, the results obtained from the ReaxFF sampling trajectories will be shown. In Figure 4, two screenshots of the model battery systems with 0.7 and 4 M TFSI within the electrolyte are shown during SEI formation. Further screenshots can be seen in the Supplementary Material (Figures S7-S9). Obviously, the simulated system is too large and too complex to draw conclusions from such visual inspections alone. The simplest way to analyze the results of the reactive dynamics is to look at how the eight different elements are distributed in space. Since we are interested mainly in the direction perpendicular to the SEI interface, i.e., the z-axis, the evaluation program packed all atoms of the respective elements into the system with 200 bins along the z-axis and separately for each written time step (one per ps, i.e., 1400 frames in total). The results are shown in Figure 5 for selected elements (O, F, and Li). Complete results, as well as simulations at different pressures and temperatures, are shown in the Supplementary Material ( Figures S10-S14). Since the system is periodic in the z-direction, essentially, two SEIs can be seen at z ≈ −18 Å and z ≈ 50 Å.
Here, and in the following, the SEI region is defined as the interval in the -direction, where a clear mixture between both initially separated phases (anode and electrolyte) can be located. This means that elements that existing initially only in the anode (Si) or in the electrolyte (F, N, …) can now be found in the same spatial region, indicating that an interface including components of both phases is formed. With respect to the element concentrations, those regions that experience significant changes (a depletion in the enrichment of elements) during the dynamics were allocated as the SEI region. In the following sections, it turns out that this definition is in good agreement with that of regions where new kinds of chemical bonds exclusive to the interface are formedT).   The simplest way to analyze the results of the reactive dynamics is to look at how the eight different elements are distributed in space. Since we are interested mainly in the direction perpendicular to the SEI interface, i.e., the z-axis, the evaluation program packed all atoms of the respective elements into the system with 200 bins along the z-axis and separately for each written time step (one per ps, i.e., 1400 frames in total). The results are shown in Figure 5 for selected elements (O, F, and Li). Complete results, as well as simulations at different pressures and temperatures, are shown in the Supplementary Material ( Figures S10-S14). Since the system is periodic in the z-direction, essentially, two SEIs can be seen at z ≈ −18 Å and z ≈ 50 Å.
Here, and in the following, the SEI region is defined as the interval in the -direction, where a clear mixture between both initially separated phases (anode and electrolyte) can be located. This means that elements that existing initially only in the anode (Si) or in the electrolyte (F, N, …) can now be found in the same spatial region, indicating that an interface including components of both phases is formed. With respect to the element concentrations, those regions that experience significant changes (a depletion in the enrichment of elements) during the dynamics were allocated as the SEI region. In the following sections, it turns out that this definition is in good agreement with that of regions where new kinds of chemical bonds exclusive to the interface are formedT). As expected, two parts of the system can be easily identified: The SiO2 anode from z ≈ −20 to z ≈ 50, where the oxygen concentration is high, and the electrolyte from z ≈ 50 Å to z ≈ 140 Å, where F and Li concentrations are high, and the oxygen concentration is lower. The time-dependent abundance of Li shows approximately horizontal filaments since only one ion is added per LiTFSI unit, building a population of just 60 in total for the case of the 0.7 M system. Lithium ions tend to penetrate the anode part for both the 0.7 and 4 M systems, where this behavior is largely enhanced at higher concentrations. The interface remains clearly visible at 0.7 M for O and F concentrations during the whole simulation time, the SEI grows to a thickness of less than 10 Å, whereas a much broader blurred region between the phases develops at 4 M. There, the SEI region with intermediate O and F concentrations is around 30 Å thick (z ≈ 30 to z ≈ 60 Å) at the end of the simulation time. Especially from the O concentration at 4 M, it appears that the SEI is formed mainly by the penetration of electrolyte atoms into the anode material, which might play the role of a sponge that sucks significant portions of the electrolyte (or its decomposition products) into it. It should be noted that the flow of Li ions from the anode to the electrolyte usually contributes to SEI formation in the formation cycles of a battery; however, this was not considered explicitly in this study. Thus, all Li ions originate from the initially placed LiTFSI units. Since our simulation time is very short in comparison to the time needed for several formation cycles, it was assumed that this contribution could Here, and in the following, the SEI region is defined as the interval in the -direction, where a clear mixture between both initially separated phases (anode and electrolyte) can be located. This means that elements that existing initially only in the anode (Si) or in the electrolyte (F, N, . . . ) can now be found in the same spatial region, indicating that an interface including components of both phases is formed. With respect to the element concentrations, those regions that experience significant changes (a depletion in the enrichment of elements) during the dynamics were allocated as the SEI region. In the following sections, it turns out that this definition is in good agreement with that of regions where new kinds of chemical bonds exclusive to the interface are formedT).
As expected, two parts of the system can be easily identified: The SiO 2 anode from z ≈ −20 to z ≈ 50, where the oxygen concentration is high, and the electrolyte from z ≈ 50 Å to z ≈ 140 Å, where F and Li concentrations are high, and the oxygen concentration is lower. The time-dependent abundance of Li shows approximately horizontal filaments since only one ion is added per LiTFSI unit, building a population of just 60 in total for the case of the 0.7 M system. Lithium ions tend to penetrate the anode part for both the 0.7 and 4 M systems, where this behavior is largely enhanced at higher concentrations. The interface remains clearly visible at 0.7 M for O and F concentrations during the whole simulation time, the SEI grows to a thickness of less than 10 Å, whereas a much broader blurred region between the phases develops at 4 M. There, the SEI region with intermediate O and F concentrations is around 30 Å thick (z ≈ 30 to z ≈ 60 Å) at the end of the simulation time. Especially from the O concentration at 4 M, it appears that the SEI is formed mainly by the penetration of electrolyte atoms into the anode material, which might play the role of a sponge that sucks significant portions of the electrolyte (or its decomposition products) into it. It should be noted that the flow of Li ions from the anode to the electrolyte usually contributes to SEI formation in the formation cycles of a battery; however, this was not considered explicitly in this study. Thus, all Li ions originate from the initially placed LiTFSI units. Since our simulation time is very short in comparison to the time needed for several formation cycles, it was assumed that this contribution could be neglected during the first merge of both phases described here (which corresponds to the very initial buildup of the battery from its parts).
The z-dependent concentrations of all nine elements are plotted for the first-and last-time steps in Figure 6 As expected, two parts of the system can be easily identified: The SiO2 anode from z ≈ −20 to z ≈ 50, where the oxygen concentration is high, and the electrolyte from z ≈ 50 Å to z ≈ 140 Å, where F and Li concentrations are high, and the oxygen concentration is lower. The time-dependent abundance of Li shows approximately horizontal filaments since only one ion is added per LiTFSI unit, building a population of just 60 in total for the case of the 0.7 M system. Lithium ions tend to penetrate the anode part for both the 0.7 and 4 M systems, where this behavior is largely enhanced at higher concentrations. The interface remains clearly visible at 0.7 M for O and F concentrations during the whole simulation time, the SEI grows to a thickness of less than 10 Å, whereas a much broader blurred region between the phases develops at 4 M. There, the SEI region with intermediate O and F concentrations is around 30 Å thick (z ≈ 30 to z ≈ 60 Å) at the end of the simulation time. Especially from the O concentration at 4 M, it appears that the SEI is formed mainly by the penetration of electrolyte atoms into the anode material, which might play the role of a sponge that sucks significant portions of the electrolyte (or its decomposition products) into it. It should be noted that the flow of Li ions from the anode to the electrolyte usually contributes to SEI formation in the formation cycles of a battery; however, this was not considered explicitly in this study. Thus, all Li ions originate from the initially placed LiTFSI units. Since our simulation time is very short in comparison to the time needed for several formation cycles, it was assumed that this contribution could be neglected during the first merge of both phases described here (which corresponds to the very initial buildup of the battery from its parts).
The z-dependent concentrations of all nine elements are plotted for the first-and lasttime steps in Figure 6

Bond Abundancies
As a first step in studying the chemical reactions happening during the SEI formation, it is instructive to look at the spatial distribution of chemical bonds. The location and species of formed or broken bonds are shown in Figure 7.
For the 0.7 M LiTFSI system, the most abundant bonds during the whole dynamics are O-Si, C-H, and C-O, where the first type is present in the anode bulk, and the other

Bond Abundancies
As a first step in studying the chemical reactions happening during the SEI formation, it is instructive to look at the spatial distribution of chemical bonds. The location and species of formed or broken bonds are shown in Figure 7.
Liquids 2023, 2, FOR PEER REVIEW 14 the anode material, and the fluoride seems to play an important role since its bonding partners vary significantly during the dynamics.

Bond Distributions
As in the case of atoms, the spatial distribution of chemical bonds also provides a good indication of the shape of the SEI. The SEI can essentially be located where new bond species are formed, connecting atoms of the electrolyte and the anode material. Due to combinatorics, many different bond species can be built during the SEI buildup; three representative examples are shown in Figure 8. More distributions are plotted in the Supplementary Material (Figures S15 and S16).
The O-Si bond shown in Figure 8a,d behaves similar to the O atoms in Figure 5. During SEI buildup, the anode consisting mainly of those bonds shrinks in size, and this behavior is more significant in the 4 M system. C-Si and S-Si bonds, on the other hand, are absent at the beginning of the dynamics since they can only be built between both electrolyte and anode phases. Their location, therefore, is a good marker for the current extent of the SEI at different times during the dynamics. As already seen before, the SEI becomes broader in the 4 M system, and it mainly grows into the former anode material (see Figure 8e,f). The 4 M LiTFSI system again shows much more dynamics and reactivity. O-Si bonds are most abundant; their number rises and reaches a maximum at approximately 200 ps, probably due to two competing processes where bulk rearrangement leads to more crystalline SiO 2 and then consumption due to SEI formation. C-H (falling) and C-O (slightly rising) bond abundancies exhibit similar time dependencies as in the 0.7 M system. F-C bonds now appear frequently as well (due to the much higher number of TFSI ions), and their number falls during the SEI formation process. The number of F-Si bonds rises from zero, in the beginning, to almost 2000 at the end. Si-Si bonds, on the other hand (zoomed plot), are consumed almost entirely, as in the 0.7 M system. C-Si interconnections rise from zero to more than 100 (in comparison to 70 in the 0.7 M system), confirming the existence of a broader SEI region. O-H and O-F bonds gain importance as well during the dynamics. Sulfur is also involved in the SEI formation, S-Si bonds (whose parameters were newly introduced into ReaxFF in line with this study) grow from zero to almost two hundred. With the larger diversity and an absolute number of different bonds, it is much more difficult to draw a deduction concerning the actual reaction mechanisms when knowing only the absolute bond numbers. The main points of attack are, again, the defects within the anode material, and the fluoride seems to play an important role since its bonding partners vary significantly during the dynamics.

Bond Distributions
As in the case of atoms, the spatial distribution of chemical bonds also provides a good indication of the shape of the SEI. The SEI can essentially be located where new bond species are formed, connecting atoms of the electrolyte and the anode material. Due to combinatorics, many different bond species can be built during the SEI buildup; three representative examples are shown in Figure 8. More distributions are plotted in the Supplementary Material (Figures S15 and S16). A broader view is given in Figure 9, where the spatial densities of the eleven most abundant bonds after 1.4 ns of dynamics are plotted in the SEI region. Similar plots covering the whole z-coordinate range can be found in the Supplemental Material ( Figure  S17). Whereas the main components of both former phases (anode: O-Si, electrolyte: C-H) mark the intersection of both phases, a larger number of less abundant species is distributed quite broadly along the SEI. The C-Si bond distribution again marks the extent of the SEI, which is roughly 2 nm thick in the 0.7 M system and 4 nm thick in the 4 M system.  Figure 5. During SEI buildup, the anode consisting mainly of those bonds shrinks in size, and this behavior is more significant in the 4 M system. C-Si and S-Si bonds, on the other hand, are absent at the beginning of the dynamics since they can only be built between both electrolyte and anode phases. Their location, therefore, is a good marker for the current extent of the SEI at different times during the dynamics. As already seen before, the SEI becomes broader in the 4 M system, and it mainly grows into the former anode material (see Figure 8e,f).
A broader view is given in Figure 9, where the spatial densities of the eleven most abundant bonds after 1.4 ns of dynamics are plotted in the SEI region. Similar plots covering the whole z-coordinate range can be found in the Supplemental Material ( Figure S17). Whereas the main components of both former phases (anode: O-Si, electrolyte: C-H) mark the intersection of both phases, a larger number of less abundant species is distributed quite broadly along the SEI. The C-Si bond distribution again marks the extent of the SEI, which is roughly 2 nm thick in the 0.7 M system and 4 nm thick in the 4 M system. A broader view is given in Figure 9, where the spatial densities of the eleven most abundant bonds after 1.4 ns of dynamics are plotted in the SEI region. Similar plots covering the whole z-coordinate range can be found in the Supplemental Material ( Figure  S17). Whereas the main components of both former phases (anode: O-Si, electrolyte: C-H) mark the intersection of both phases, a larger number of less abundant species is distributed quite broadly along the SEI. The C-Si bond distribution again marks the extent of the SEI, which is roughly 2 nm thick in the 0.7 M system and 4 nm thick in the 4 M system.

Molecule Abundancies and Distributions
Even more interesting than the number and location of chemical bonds in the systems is how many molecules of individual species are present during SEI buildup at different parts of the interface. In this study, this analysis concentrated on small gas phase mole-

Molecule Abundancies and Distributions
Even more interesting than the number and location of chemical bonds in the systems is how many molecules of individual species are present during SEI buildup at different parts of the interface. In this study, this analysis concentrated on small gas phase molecules with up to seven atoms in them. How the molecules were determined from the location of the chemical bonds, by using openbabel [71] for the generation and comparison of SMILES strings [72] for extracted local geometries, is explained in the Supplementary Material (Algorithm S1). For both 0.7 and 4 M systems, the number of those small molecules was counted for each MD screenshot (every 0.5 ps). The most abundant molecules, summed up over all the time frames, are shown in Figure 10, with their numbers in the final time frame indicated as well (see also Supplementary Material, Figures S18 and S19).
Ethylene is the most abundant small molecule in the 0.7 M system, followed by methanol, formaldehyde, hydrogen, carbon dioxide, and hydrogen fluoride. If we compare this with (b), where the same list is given for the 4 M LiTFSI system, we see many similarities but also some obvious differences. First, the total number of small molecules is much larger; the most abundant species are counted nearly 50 times in the last frame, compared to 19 and 12 in the 0.7 M system. The ordering is also somewhat different: the abundance of ethylene is now only the third highest, with carbon dioxide and formaldehyde taking the top positions. Hydrogen fluoride, methyl hypofluoride, and water follow them. Since, as another study suggested, ethylene is one of the main degradation products of DOL [73], this finding coincides well with expectations. As ethylene is a volatile gaseous by-product, it can be expected that, over longer time scales, it will escape the SEI. cules with up to seven atoms in them. How the molecules were determined from the location of the chemical bonds, by using openbabel [71] for the generation and comparison of SMILES strings [72] for extracted local geometries, is explained in the Supplementary Material (Algorithm S1). For both 0.7 and 4 M systems, the number of those small molecules was counted for each MD screenshot (every 0.5 ps). The most abundant molecules, summed up over all the time frames, are shown in Figure 10, with their numbers in the final time frame indicated as well (see also Supplementary Material, Figures S18 and S19).  Ethylene is the most abundant small molecule in the 0.7 M system, followed by methanol, formaldehyde, hydrogen, carbon dioxide, and hydrogen fluoride. If we compare this with (b), where the same list is given for the 4 M LiTFSI system, we see many similarities but also some obvious differences. First, the total number of small molecules is much larger; the most abundant species are counted nearly 50 times in the last frame, compared to 19 and 12 in the 0.7 M system. The ordering is also somewhat different: the abundance of ethylene is now only the third highest, with carbon dioxide and formaldehyde taking the top positions. Hydrogen fluoride, methyl hypofluoride, and water follow them. Since, as another study suggested, ethylene is one of the main degradation products of DOL [73], this finding coincides well with expectations. As ethylene is a volatile gaseous by-product, it can be expected that, over longer time scales, it will escape the SEI.
In order to study the reaction dynamics, which leads to the formation of the different molecular species, the time-dependent numbers of the most abundant ones are plotted in Figure 11. In the 0.7 M system, we see that the formation of ethylene is mainly finished In order to study the reaction dynamics, which leads to the formation of the different molecular species, the time-dependent numbers of the most abundant ones are plotted in Figure 11. In the 0.7 M system, we see that the formation of ethylene is mainly finished after 0.6 ns, and other less abundant species then also have essentially constant numbers; only the number of formaldehyde and oxygen molecules slightly decreases after 0.7 ns. This indicates that the SEI formation is essentially finished after 0.6-0.7 ns, with only subordinated reactions consuming some formaldehyde and oxygen still taking place. (Of course, rare reaction events would not be visible on this time scale but could lead to appreciable changes on a millisecond time scale). The picture is again much more interesting for the 4 M system. Here, no real equilibrium is reached in the simulated 1.4 ns; only the number of ethylene molecules reaches a plateau after 0.4 ns. Formaldehyde and carbon dioxide, which are less abundant between 0.4 ns and 0.6 ns, are still being produced, leading them to outplay ethylene. Fluoride species, such as methyl hypofluoride, hydrogen fluoride, and hyperfluorous acid, that are produced from the decomposition of TFSI anions, are important as well; hydrogen fluoride is produced first and reaches a plateau after 0.8 ns, whereas the number of hyperfluorous acid molecules rises until the end of the dynamic.
The rather short simulation time was sufficient to reach the first kind of chemical equilibrium in the SEI region for the 0.7 M system. For the 4 M system, on the other hand, many processes are still not complete, and even longer simulations are needed. Since, however, the general trends concerning the elemental compositions of the SEI are already in qualitative agreement with the experimental findings, and a complete simulation of the buildup seems to be out of reach even for computationally cheaper methods such as RE-ACTER [44], it follows that our approach was justified for giving reasonable mechanistic insights into the process.  The rather short simulation time was sufficient to reach the first kind of chemical equilibrium in the SEI region for the 0.7 M system. For the 4 M system, on the other hand, many processes are still not complete, and even longer simulations are needed. Since, however, the general trends concerning the elemental compositions of the SEI are already in qualitative agreement with the experimental findings, and a complete simulation of the buildup seems to be out of reach even for computationally cheaper methods such as REACTER [44], it follows that our approach was justified for giving reasonable mechanistic insights into the process.
Additional information can be gained from the spatial distribution of the molecules. Since the number of molecular species present at one point in time was quite small, the total trajectory was split into five parts, for which molecular spatial abundances were summed up, so that those time-resolved and also sufficiently space-resolved distributions were gained.
These distributions are plotted in Figure 12a,b for the 0.7 M system (see also Supplementary Material, Figures S20 and S21). As could be expected from the previous findings, there is no significant variation between (a) and (b) since the SEI formation is essentially finished after ca. 0.4 ns. All molecules besides oxygen can only be found in the SEI region, and the oxygen originates presumably from local defects in the amorphous SiO 2 structure, where excess oxygen atoms dimerize. The SEI region itself is somewhat shifted towards larger z values when comparing it to the bond concentrations shown in Figure 8. There we see that the SEI spreads from z = 40 to z = 60. Small molecules are thus preferentially formed at the electrolyte side of the SEI. The anode side presumably consists mainly of larger networks of molecule parts attached to the dangling SiO 2 chains (see below). Within the SEI region consisting of small molecules, a fine structure can be found with the methanol molecules nearer to the anode side (z = 50 to z = 57) and the formaldehyde molecules nearer to the electrolyte side. The most abundant ethylene can be most probably found in the center of the molecular SEI.
the SEI region consisting of small molecules, a fine structure can be found with the methanol molecules nearer to the anode side (z = 50 to z = 57) and the formaldehyde molecules nearer to the electrolyte side. The most abundant ethylene can be most probably found in the center of the molecular SEI.
If those findings are compared with the time-dependent distributions for the 4 M system in Figure 12c,d the SEI there is much broader. In the first fifth of the trajectory, the SEI covers essentially the same region as for the 0.7 M system, but it grows further, and finally, small molecules can be found in almost the entire electrolyte region. Concerning the different species, carbon dioxide is preferably located near the anode, whereas formaldehyde, ethylene, and methanolate are spread far into the electrolyte region. We can thus follow that the SEI has a layered structure: small molecules are preferentially located at the electrolyte side, with different molecules existing at different layers of the molecular SEI, in contrast to the rather solid-state SEI at the anode side.

Reaction Mechanisms
A better impression of the reaction mechanisms can be gained if trajectories are visually inspected. Since the systems contain more than 10,000 atoms, no systematic evaluation is possible in that way. Instead, the decomposition of the three electrolyte molecules, DOL, DME, and (Li)TFSI, was tracked for some examples. One example of a DOL molecule in the 0.7 M TFSI system is shown in Figure 13. Initially, the molecule is located some Å away from the electrolyte-anode interface. After 90 ps, it has reached the already growing SEI (beige Si atoms in the background). A proton ejected from some other SEI reaction attacks one of the ether bridges at 91 ps. The DOL molecule then decomposes into one ethylene: formaldehyde and a FOH molecule some ps later. The FOH merges with a methyl group at 98 ps and becomes a methanol molecule after the F atom has left. If this is compared with the molecule abundances in Figure 14, ethylene, methanol, and formaldehyde are indeed the most common small molecules, indicating that the mechanism shown could be an important one. If those findings are compared with the time-dependent distributions for the 4 M system in Figure 12c,d the SEI there is much broader. In the first fifth of the trajectory, the SEI covers essentially the same region as for the 0.7 M system, but it grows further, and finally, small molecules can be found in almost the entire electrolyte region. Concerning the different species, carbon dioxide is preferably located near the anode, whereas formaldehyde, ethylene, and methanolate are spread far into the electrolyte region. We can thus follow that the SEI has a layered structure: small molecules are preferentially located at the electrolyte side, with different molecules existing at different layers of the molecular SEI, in contrast to the rather solid-state SEI at the anode side.

Reaction Mechanisms
A better impression of the reaction mechanisms can be gained if trajectories are visually inspected. Since the systems contain more than 10,000 atoms, no systematic evaluation is possible in that way. Instead, the decomposition of the three electrolyte molecules, DOL, DME, and (Li)TFSI, was tracked for some examples. One example of a DOL molecule in the 0.7 M TFSI system is shown in Figure 13. Initially, the molecule is located some Å away from the electrolyte-anode interface. After 90 ps, it has reached the already growing SEI (beige Si atoms in the background). A proton ejected from some other SEI reaction attacks one of the ether bridges at 91 ps. The DOL molecule then decomposes into one ethylene: formaldehyde and a FOH molecule some ps later. The FOH merges with a methyl group at 98 ps and becomes a methanol molecule after the F atom has left. If this is compared with the molecule abundances in Figure 14, ethylene, methanol, and formaldehyde are indeed the most common small molecules, indicating that the mechanism shown could be an important one.

Reaction Mechanisms
A better impression of the reaction mechanisms can be gained if trajectories are visually inspected. Since the systems contain more than 10,000 atoms, no systematic evaluation is possible in that way. Instead, the decomposition of the three electrolyte molecules, DOL, DME, and (Li)TFSI, was tracked for some examples. One example of a DOL molecule in the 0.7 M TFSI system is shown in Figure 13. Initially, the molecule is located some Å away from the electrolyte-anode interface. After 90 ps, it has reached the already growing SEI (beige Si atoms in the background). A proton ejected from some other SEI reaction attacks one of the ether bridges at 91 ps. The DOL molecule then decomposes into one ethylene: formaldehyde and a FOH molecule some ps later. The FOH merges with a methyl group at 98 ps and becomes a methanol molecule after the F atom has left. If this is compared with the molecule abundances in Figure 14, ethylene, methanol, and formaldehyde are indeed the most common small molecules, indicating that the mechanism shown could be an important one. The example mechanism for a DME molecule is given in Figure 14. The molecule is located near the interface at the beginning; after 10 ps, one ether oxygen is attacked by a Liquids 2023, 2, FOR PEER REVIEW 20 silicon atom, probably with a dangling bond. Shortly after, the bridge is broken, a methanol substituent is bonded to the silicon network, and a methoxyethane radical is built. A total of 4 ps later, this radical decomposes into an ethylene and a methanol molecule, where probably, an external proton is involved in breaking the ether bond since the products include one additional hydrogen. Again, all products belong to the most abundant molecules, making this mechanism likely to be the main one. Finally, an example TFSI decomposition mechanism is shown in Figure 15. In the beginning, it is already placed in the immediate vicinity of the anode surface. Then, at 1 ps later, a silicon atom attacks one sulfur atom, and shortly after, it is bound to two sulfur atoms in a tetrahedral arrangement. Two ps later, both CF3 groups disaggregated. The central N(SO2)2 group essentially remains intact, successively absorbed deeper into the amorphous SiO2 network, noticeable by the larger and larger number of silicon atoms that are bound to the group from different sides. One carbon and one oxygen atom become part of a carbon dioxide molecule (visible at 9 ps). The fluorine atoms, on the other hand, form hydrogen fluoride or react with oxygen atoms from the surrounding (76 ps). The example mechanism for a DME molecule is given in Figure 14. The molecule is located near the interface at the beginning; after 10 ps, one ether oxygen is attacked by a silicon atom, probably with a dangling bond. Shortly after, the bridge is broken, a methanol substituent is bonded to the silicon network, and a methoxyethane radical is built. A total of 4 ps later, this radical decomposes into an ethylene and a methanol molecule, where probably, an external proton is involved in breaking the ether bond since the products include one additional hydrogen. Again, all products belong to the most abundant molecules, making this mechanism likely to be the main one.
Finally, an example TFSI decomposition mechanism is shown in Figure 15. In the beginning, it is already placed in the immediate vicinity of the anode surface. Then, at 1 ps later, a silicon atom attacks one sulfur atom, and shortly after, it is bound to two sulfur atoms in a tetrahedral arrangement. Two ps later, both CF 3 groups disaggregated. The central N(SO 2 ) 2 group essentially remains intact, successively absorbed deeper into the amorphous SiO 2 network, noticeable by the larger and larger number of silicon atoms that are bound to the group from different sides. One carbon and one oxygen atom become part of a carbon dioxide molecule (visible at 9 ps). The fluorine atoms, on the other hand, form hydrogen fluoride or react with oxygen atoms from the surrounding (76 ps). Finally, an example TFSI decomposition mechanism is shown in Figure 15. In the beginning, it is already placed in the immediate vicinity of the anode surface. Then, at 1 ps later, a silicon atom attacks one sulfur atom, and shortly after, it is bound to two sulfur atoms in a tetrahedral arrangement. Two ps later, both CF3 groups disaggregated. The central N(SO2)2 group essentially remains intact, successively absorbed deeper into the amorphous SiO2 network, noticeable by the larger and larger number of silicon atoms that are bound to the group from different sides. One carbon and one oxygen atom become part of a carbon dioxide molecule (visible at 9 ps). The fluorine atoms, on the other hand, form hydrogen fluoride or react with oxygen atoms from the surrounding (76 ps).

Electrochemical Characterization
In order to investigate the SEI formation on the silicon active material surface, pure porous silicon anodes were cycled in the presence of varying lithium salt concentrations. The anodes were cycled in half-cells with low current densities (0.1 C) against a Li-metal counter and reference electrode and subjected to galvanostatic charge/discharge tests. The results are shown below in Figure 16.

Electrochemical Characterization
In order to investigate the SEI formation on the silicon active material surface, pure porous silicon anodes were cycled in the presence of varying lithium salt concentrations. The anodes were cycled in half-cells with low current densities (0.1 C) against a Li-metal counter and reference electrode and subjected to galvanostatic charge/discharge tests. The results are shown below in Figure 16. The initial formation cycle is shown in Figure 16a, where the samples are cycled with a CC-CV mode until a cut-off potential of 3150 mAh g −1 . For the anode cycled in 0.7 M LiTFSI, there is a small plateau observed during charging around 2.0 V, and a second, much larger plateau beginning at 1.3 V was much higher than the potential for the lithiation of silicon (which normally begins at 0.8 V [4,74]). The first plateau around 2.0 V is The initial formation cycle is shown in Figure 16a, where the samples are cycled with a CC-CV mode until a cut-off potential of 3150 mAh g −1 . For the anode cycled in 0.7 M LiTFSI, there is a small plateau observed during charging around 2.0 V, and a second, much larger plateau beginning at 1.3 V was much higher than the potential for the lithiation of silicon (which normally begins at 0.8 V [4,74]). The first plateau around 2.0 V is caused by electrolyte components beginning to deposit on the surface of the anode [75] and likely corresponds to the peak observed for both electrolytes in the cyclic voltammetry measurements shown in the Supplemental Material ( Figure S22). This plateau may also correspond to a reduction in trace H 2 O and other trace oxygen-containing impurities [76], although these would have a negligible contribution to the overall SEI composition. Below 1.5 V, the plateau is likely caused by electrolyte decomposition and SEI formation reactions. This plateau begins at lower voltages for the 4 M electrolyte (1.2 V), which is characteristic of highly concentrated electrolytes [77].
In the initial formation cycle (Figure 16a), the reversible capacity of the silicon electrode in 4 M LiTFSI is much higher (1226.4 mAh/g) than in the 0.7 M electrolyte (476.1 mAh/g). This is because irreversible electrolyte degradation, resulting in SEI formation, contributes to a more significant amount of the charge capacity measured for the 0.7 M electrode than the 4 M electrode in the first cycle. Therefore, there is less lithium present in this electrode for the delithiation process.
Though the discharge capacities after the formation cycles of the anodes cycled in both electrolytes are both relatively high after ten cycles, the electrolyte containing 4 M LiTFSI obtained the highest discharge capacity of 2830.2 mAh g −1 Si (Figure 16b), compared to the discharge capacity of 2463.3 mAh g −1 Si obtained by the 0.7 M electrolyte. In addition, in the first cycle, a higher coulombic efficiency was obtained with the 4 M electrolyte (38.9%) compared to the 0.7 M electrolytes (15.1%) (Figure 16c) (related to the higher reversible capacities observed in Figure 16a). Even after reaching a stable point, the coulombic efficiency for both electrodes was still below 100%, implying that there would still be irreversible losses for both electrodes during extended cycling. However, the coulombic efficiency for the electrode in the 4 M electrolyte (~90%) is higher than for the electrode in the 0.7 M electrolyte (~82%). This means that lower irreversible losses and more stable cycling can be expected for this electrolyte.

Post-Mortem Investigation
In order to further investigate the chemical composition of the SEI layer, porous silicon anodes, which cycled in the electrolytes with both LiTFSI concentrations (4 M and 0.7 M), were investigated ex-situ with X-ray photoelectron spectroscopy after one cycle. The results are shown below in Figure 17.
The element percentages, calculated from the resolved XPS spectra, are shown in Figure 17a. From this, it is clear that the anode cycled with 4 M LiTFSI possesses much higher percentages of fluorine, sulfur, and nitrogen, all elements unique to the salt anion, and lower percentages of carbon and oxygen. What is particularly interesting is the higher silicon percentage. As the XPS has an attenuation depth of around 10 nm, this seems to indicate that the insoluble SEI for the 4 M sample is thinner than that of the 0.7 M sample. This apparent contradiction to the ReaxFF results is addressed in the following section. The resolved C 1s spectra (Figure 17c) show carbon species originating from the decomposition of the solvent molecules in the binding energy range of 285-290 eV [78][79][80]. The peak at 293 eV, originating from the -CF 3 group of the TFSI − anion, is seen to increase in intensity for the silicon anode cycled in the 4 M LiTFSI electrolyte. This indicates the presence of more species in the SEI originating from the TFSI − anion. This is confirmed with the F 1s spectra (Figure 17b). The broad amorphous peak appearing in the Si 2p resolved spectra (Figure 17d) appears at higher binding energies for the silicon anode cycled with 4 M LiTFSI (101.8 eV compared to 100.6 eV). There may be more lithium silicates, such as Li x SiO y , or organosilicates, in the SEI of the silicon anode cycled with 0.7 M LiTFSI [79], while more inorganic species are formed with 4 M LiTFSI. This, combined with the high Si atomic percentage, would indicate that more silicon species are incorporated into the SEI for the 4 M electrolyte than for the 0.7 M electrolyte. Details of the relative peak percentages of different species observed in the resolved XPS spectra are shown in Table S2 and the  elemental percentages are shown in Table S3 in the Supplementary Material. lombic efficiency for the electrode in the 4 M electrolyte (~90%) is higher than for the electrode in the 0.7 M electrolyte (~82%). This means that lower irreversible losses and more stable cycling can be expected for this electrolyte.

Post-Mortem Investigation
In order to further investigate the chemical composition of the SEI layer, porous silicon anodes, which cycled in the electrolytes with both LiTFSI concentrations (4 M and 0.7 M), were investigated ex-situ with X-ray photoelectron spectroscopy after one cycle. The results are shown below in Figure 17. The element percentages, calculated from the resolved XPS spectra, are shown in Figure 17a. From this, it is clear that the anode cycled with 4 M LiTFSI possesses much higher percentages of fluorine, sulfur, and nitrogen, all elements unique to the salt anion, and lower percentages of carbon and oxygen. What is particularly interesting is the higher silicon percentage. As the XPS has an attenuation depth of around 10 nm, this seems to indicate that the insoluble SEI for the 4 M sample is thinner than that of the 0.7 M sample. This apparent contradiction to the ReaxFF results is addressed in the following section. The resolved C 1s spectra (Figure 17c) show carbon species originating from the decomposition of the solvent molecules in the binding energy range of 285-290 eV [78][79][80]. The peak at 293 eV, originating from the -CF3 group of the TFSI − anion, is seen to increase in intensity for the silicon anode cycled in the 4 M LiTFSI electrolyte. This indicates the presence of more species in the SEI originating from the TFSI − anion. This is confirmed with the F 1s spectra (Figure 17b). The broad amorphous peak appearing in the Si 2p resolved spectra (Figure 17d) appears at higher binding energies for the silicon anode cycled with 4 M LiTFSI (101.8 eV compared to 100.6 eV). There may be more lithium silicates, such as LixSiOy, or organosilicates, in the SEI of the silicon anode cycled with 0.7 M LiTFSI [79],

Discussion
During typical SEI formation processes, the products of salt anion reduction reactions tend to be inorganic compounds (e.g., LiF or Li 2 O), while the reduction in solvent molecules results in the formation of insoluble SEI components (e.g., Li 2 CO 3 ) or partially soluble semi-carbonates and polymers [81]. The SEI structure tends to have inorganic, insoluble components accumulating at the electrode surface, with the components becoming more polymeric towards the electrolyte [82]. This expected SEI structure was seen both in the ReaxFF simulations, as well as the post-mortem XPS investigations. What is interesting is how the SEI structure differs for high LiTFSI salt concentrations, as this directly influences the cycling performance of the silicon anodes.
From the galvanostatic cycling results, it was seen that the formation was completed much faster for the 4 M LiTFSI electrolyte since the reversible capacities are much higher ( Figure 16a) and a stable coulombic efficiency is reached sooner than for the 0.7 M electrolyte (Figure 16c). In the ReaxFF simulations, Li penetration into the anode was enhanced at higher concentrations, and the SiO x interface appeared blurred. For low salt concentrations, the SiO x interface remained clearly visible for much longer (Figure 5), suggesting slower SEI formation and lower surface interactions, corresponding to the results observed from the experimental cycling. The broad interface and high penetration of electrolyte atoms into the anode material observed from the ReaxFF simulations in Figures 5 and 6 support the theory that more silicon species are incorporated into the SEI for 4 M LiTFSI (Figure 17). At high salt concentrations, the O-Si bonds are more abundant and increase over time during their formation (Figure 7). In both the ReaxFF simulations and the XPS results, F plays a significant role in the reactions for the 4 M electrolyte, and the species in the SEI are observed to be more inorganic, with a higher number of smaller molecules throughout the SEI ( Figure 11). As the organic C species observed in the resolved XPS significantly decreased for the 4 M electrolyte (Figure 17), it could be concluded that there is also a lower polymerization or other reactions associated with a reduction in the solvent molecules.
In the highly concentrated electrolyte, the SEI appears to grow into the anode surface (Figure 8), coinciding with the observation that more silicon species are incorporated into the SEI. From the cycling results, the passivation of the porous silicon anode occurs more quickly with high Li-salt concentrations than with lower concentrations. This is evident from the fact that the CE for the 4 M electrolyte was higher than the 0.7 M electrolyte in the initial cycle (38.9% and 15.1%, respectively) since SEI formation reactions contribute more to the initial charge capacity for the 0.7 M electrolyte. For both electrolytes, the native oxide on the surface is directly involved in SEI formation reactions (shown by the consumption of Si-Si bonds in Figure 7). Additionally, the SEI for both electrolytes appears to be constructed from inorganic decomposition products directly at the surface of the anode, where the first layer mainly involves SiO x species and organic decomposition products at the electrolyte side. A similar SEI construction was observed in previous studies [83,84]. After formation for the 4 M electrolyte, the organic decomposition products at the electrolyte side were minimal, corresponding to the highly inorganic SEI observed with XPS. In the ReaxFF simulations, the SEI was thinner in the case of 0.7 M TFSI, in contrast to the XPS data, which suggests that the 4 M SEI is thinner. This is probably due to the fact that the time scale of the ReaxFF simulations (1.4 ns) is still too short for a complete simulation of the converging growth process. In addition to this, comparatively slow transport processes might, e.g., move gaseous species away from the SEI, which still constitutes a significant part of the SEI in the 4 M model. Nevertheless, the ReaxFF simulations were able to reproduce many of the trends with respect to the elemental composition of the interface.
The SEI formation mechanism is likely also linked to the different solvation structures in the low and highly-concentrated electrolytes. The solvation structure of the electrolyte is made up of solvent-separated ion pairs (SSIPs), contact ion pairs (CIPs), and aggregates (AGG), and the latter dominate at high salt concentrations. In the literature, CIPs (when ions are in contact with each other) have been shown to be preferentially reduced prior to SSIPs (where ions are separated by solvent molecules) [85]. As seen in the work by Yamada et al. [86], with higher salt concentrations, there is a change in the lowest unoccupied molecular orbital (LUMO) of the electrolyte from the solvent molecule (for low salt concentrations) to the TFSI − anion. More CIPs were observed in the 4 M electrolyte, which can be observed with Raman spectroscopy (see Figures S23 and S24 in the Supplementary Material). Since the CIPs are preferentially reduced at the anode surface, the early passivation occurring in the first cycles for the 4 M electrolyte is due to this reduction in the CIPs on the surface. Since more CIPs are present at high concentrations, this reduction can occur more quickly, and the passivation of the surface is reached earlier.

Conclusions
The ReaxFF simulations were successfully completed, and, combined with the results from the experimental verification, it can be concluded that penetration of lithium into the silicon surface occurs more easily for highly concentrated electrolytes. At high LiTFSI concentrations, more SiO x species are incorporated into the SEI, and more TFSI − anions take part in SEI formation reactions, allowing for more efficient SEI formation, as well as a more inorganic SEI layer. This leads to higher coulombic efficiencies being reached in the initial formation cycle (38.9%, compared to 15.1% for the 0.7 M electrolyte) and higher discharge capacities after formation (2830.2 mAh g −1 Si and 2463.3 mAh g −1 Si for 4 M and 0.7 M LiTFSI, respectively). From these results, high LiTFSI salt concentrations appear to enhance the performance of silicon anodes; however, further investigations are still required to determine the implications for long-term cycling and for full-cell applications with sulfur cathodes.
This ReaxFF approach was shown to correspond well with the experimental results and offers the possibility of gaining a deeper understanding of SEI formation mechanisms and predicting the performance of silicon anodes with different electrolyte chemistries without undergoing time-consuming experimental testing. Arbitrary spatial and temporal resolutions can be gained during the simulations, allowing for the characterization of different sublayers in the SEI in the vicinity of the anode or the electrolyte. Further, molecular reaction products can be identified and classified automatically, elucidating the complex chemical mechanisms during the formation. The general setup of the ReaxFF optimization, as well as its application during reactive molecular dynamics, together with the evaluations conducted by utility scripts developed with the focus on reactive interfaces, facilitate future simulations of SEI formations in different LIBs, with only moderate translation efforts required.