The solid-electrolyte interphase (SEI) layer is a critical constituent of battery technology, which incorporates the use of lithium metals. Since the formation of the SEI is difficult to avoid, the engineering and harnessing of the SEI are absolutely critical to advancing energy storage. One problem is that much fundamental information about SEI properties is lacking due to the difficulty in probing a chemically complex interfacial system. One such property that is currently unknown is the dissolution of the SEI. This process can have significant effects on the stability of the SEI, which is critical to battery performance but is difficult to probe experimentally. Here, we report the use of ab initio computational chemistry simulations to probe the solution state properties of SEI components LiF, Li2O, LiOH, and Li2CO3 in order to study their dissolution and other solution-based characteristics. Ab initio molecular dynamics was used to study the solvation structures of the SEI with a combination of radial distribution functions, discrete solvation structure maps, and vibrational density of states, which allows for the determination of free energies. From the change in free energy of dissolution, we determined that LiOH is the most likely component to dissolve in the electrolyte followed by LiF, Li2CO3, and Li2O although none were favored thermodynamically. This indicates that dissolution is not probable, but Li2O would make the most stable SEI with regard to dissolution in the electrolyte.

Lithium-based energy storage technology has been and continues to be massively utilized. From personal electronics to larger scale uses such as off-grid energy storage, there is a wide demand for the convenience and storage capabilities of high-performance battery technology.1 The most common type of lithium battery technology is the lithium ion battery. This battery system has energy densities of ∼200–500 Wh/L, which makes the system quite advantageous for many applications.2,3 However, there is a large demand for greater performance from energy storage technologies in the form of more storage capacities and greater cycle life.3,4 Given the general requirements, there are many possible battery candidates that may be able to meet these goals. On the other hand, developing novel energy storage systems is greatly challenging. One type of next generation energy storage technology that meets higher performance and does not involve completely reinventing the wheel is to use systems such as lithium–sulfur or lithium–air batteries, which still use lithium technology but not the lithium ion system specifically.5 These systems utilize other cathodes with much greater capacity (the theoretical capacity of 3860 mAhg−1), but in order to access that higher capacity, a lithium metal anode must be used.6,7

The lithium metal is highly reactive, which makes attempts to use it in a stable battery challenging. As the electrolyte comes into contact with the electrode, lithium begins to react with the electrolyte components, decomposing them into a mixture of different chemical species, which coat the lithium metal forming the solid–electrolyte interphase (SEI) layer.8–10 The SEI is a difficult to avoid reality of utilizing the lithium metal in a battery technology,11 which has strong impacts on the performance of the battery in different ways. The SEI can passivate the electrode, but certain degrees of passivation may impede the necessary electrochemical reactions that make a battery useful. Electronic insulation and/or ionic insulation is not a desired quality for an electrode to have. The mechanical stability of the SEI is also closely related to the growth/suppression of lithium dendrites.12,13 Research has suggested that SEI cracking is one of the nucleation points for the growth of dendrites, which introduces both performance and safety challenges with lithium metal battery operation.14,15 Moreover, the chemical stability of the SEI, i.e., its resistance to dissolution or chemical decomposition, is key to the SEI passivation role.16,17

In order to engineer and harness the SEI to create a practical battery, the fundamental properties and processes of the layer must be characterized. There are many different aspects of the SEI that must be examined to achieve a complete picture such as the mechanical properties to understand cracking.18,19 There are both experimental and computational approaches to characterizing the properties of the SEI, which then further guide researchers into developing a SEI, which does not affect battery performance.18,20–22 One fundamental process is dissolution of the SEI into the electrolyte.23 Anytime there is a solid phase in contact with a liquid phase, there is a possibility that the solid phase may dissolve into the liquid phase. If an SEI component is prone to dissolution, then it would not provide the desired stability for a lithium metal battery with long cycle life.

Native layers of inorganic compounds (Li2CO3, LiOH, and Li2O) are always present in the lithium foils due to environmental contamination24 and play an important role in the initial SEI layers formed by oxidation of the Li foils in contact with the electrolyte, and LiF is a very common component of the SEI that appears even earlier than the organic SEI products because of the higher reduction potential of the F-containing anions. Unfortunately, it can be difficult to experimentally observe and investigate solution state properties as fine as solvation structure and dissolution dynamics. This makes for a good use of computational chemistry methods, which have been long developed for studying liquid systems in order to close the gap of understanding.25 Solubility and dissolution of SEI components have been investigated before mainly with classical molecular dynamics and in systems focusing mainly on lithium ion battery electrolytes.16,17 These studies provide useful information but with forcefield approximations and not as applicable to the lithium–sulfur battery system. However, with ab initio computational chemistry methods, more accurate predictions and observations can be made of the solution side of the SEI interface, which can provide critical information regarding the stability of the SEI, dissolution, and nucleation. In this report, we investigate the solution state properties of LiF, LiOH, Li2O, and Li2CO3 by using ab initio molecular dynamics (AIMD). The solubilities are obtained from free energy differences between the solid and liquid states. The free energies in the system are calculated from the vibrational density of states (VDOS) using a technique based on the velocity autocorrelation function from ab initio molecular dynamics simulations.

The progressions of the systems over time were achieved with ab initio molecular dynamics (AIMD) simulations as implemented by the Vienna ab initio simulation package (VASP). VASP is a first-principles simulation program that incorporates periodic boundary conditions (PBCs) and with a high degree of performance. The periodic boundary conditions in this case are used to represent the bulk electrolyte by repeating the cells infinitely in the x, y, and z directions. Only one cell needs to be simulated in order to obtain the behavior of the entire bulk. The electron–nucleus energy terms are incorporated with the projector augmented wave pseudopotentials, while electronic correlation uses the Perdew–Burke–Ernzerhof (PBE) functional. This functional utilizes the generalized gradient approximation (GGA) and has been well studied and referenced in the literature. Density functional theory (DFT) was used to determine the electronic structure and the subsequent forces for the AIMD. The canonical (NVT) ensemble was used at 300 K with a 0.5 damping parameter Nose thermostat used for temperature control. An electronic convergence parameter was set at 1 × 10−6 eV, and a 400 eV cutoff was selected. Surface Brillouin zone integration was executed with a (2 × 2 × 2) Monkhorst–Pack k-point grid. For all the simulations in this work, a time step of 1 fs was used and 5 ps of each simulation were reserved for equilibration time, while the remaining time of the simulations was used for production time.

The systems were constructed by placing an SEI molecule in a 10 × 10 × 10 Å cell. The cell was then amorphously packed with the electrolyte to complete the representation of the solution phase. The electrolyte was represented by using pure 1,3-dimethoxyethane (DME) with the studied species in question: LiF, LiOH, Li2CO3, and LiF. A pure DME electrolyte was chosen to limit the size of the PBC box size to a computationally feasible level. In a conventional electrolyte, lithium salts, additives, and a mixture of solvents would be used, but based on previous studies of solvation in similar electrolytes, one would expect the majority of the ionic species’ primary solvation shell to be DME. Thus, any SEI fragment that dissolves is expected to be solvated by the solvent and less by the other electrolyte components. Besides the simulations of one SEI molecule (ion pair) in the solvent (∼1.5M corresponding to one molecule in a 10 × 10 × 10 Å box), we evaluated the structure in solution of dissociated molecules and small clusters. The dissociated and cluster simulations used the same box size and cell parameters but changed the initial configuration of the simulation. Finally, solid-state systems were also characterized. The AIMD simulations of solid-state systems were built using bulk crystal structures from online databases.

For the estimation of the free energies of the AIMD simulation, a previously established technique based on the vibrational density of states (VDOS) or vibrational spectra was used.26–30 The VDOS is equivalent to a density of states (DOS), but instead of the broadening of electronic levels due to small perturbations, the vibrational modes are broadened. The VDOS for each system is determined from the AIMD trajectories by calculating the velocity autocorrelation function for each atom in the system and subsequently applying a Fourier transform. This determines the atom spectral density or power spectrum, which is then mass-averaged across the system. The resulting function represents the population of the system at different vibrational frequencies, which can then be used with traditional statistical mechanics formulas to estimate the energy due to translation and vibration, which allow the free energy to be determined from the ground state energy.31 Scheme 1 displays a flowchart of the type of calculations done for obtaining the free energy and the corresponding equations.

SCHEME 1.

(a) Flowchart for obtaining free energy from AIMD simulation. (b) Equations used in the process. For more details in the methodology, we refer the readers to the original reference for this method (Refs. 29–31).

SCHEME 1.

(a) Flowchart for obtaining free energy from AIMD simulation. (b) Equations used in the process. For more details in the methodology, we refer the readers to the original reference for this method (Refs. 29–31).

Close modal

In order to characterize the stability of the passivation layer of lithium in terms of solubility in the solvent, an accurate picture of the SEI components in the solid and liquid phases must be created. The thermodynamic favorability of the SEI components in the liquid phase vs the solid phase is described by the differences in the configurations and energetics between the states. Therefore, the better the representation of the liquid and solid states is, the more accurate the description of the dissolution becomes. In the case of a bulk solid structure, relatively few assumptions need to be made about the structure as opposed to the liquid state. The solvation structure, ionic association status, and energy of interactions between the solvate and the solvent are critical knowledge required to accurately describe the solution state of ionic dissolved SEI components. We first aim to create an accurate representation and understanding of the structure of dissolved SEI components in a DME solvent to characterize the possible dissolution of the SEI into the electrolyte.

1. Contact ion pairs

The first systems investigated were LiF, LiOH, Li2CO3, and Li2O molecules (contact ion pairs) in DME, as shown in Fig. 1. Each system was run for ∼15 ps in order to achieve equilibration. The time evolution of the energy and temperature of each system (see the supplementary material, system variables, Figs. S-46–S-62) were used to confirm that systems had approached an equilibrated state.

FIG. 1.

Initial contact ion pairs for (a) LiF, (b) LiOH, (c) Li2CO3, and (d) Li2O with DME.

FIG. 1.

Initial contact ion pairs for (a) LiF, (b) LiOH, (c) Li2CO3, and (d) Li2O with DME.

Close modal

First, we explored the structure of these liquid systems with radial distribution functions (RDFs). The Li–X RDFs for all tested species are displayed in Fig. 2. For the LiF system, there is a strong peak less than 2 Å that is from the fluorine atom. This is indicative of the Li–F bond of the SEI molecule. From observing the distance directly (Fig. S-1), the Li–F bond oscillates slightly, but the average is near 1.7 Å. In the LiOH system, a similar behavior is encountered. The first and strongest peak is a double peak between the Li and oxygen atoms in the system. The first peak is centered around ∼1.5 Å, while the second peak is centered in ∼2 Å. The first distance is confirmed as the Li–O bond from the LiOH molecule in Fig. S-2. Thus, the second peak must be caused by the DME–oxygen atom coordination. In all the tested systems, there is a sharp Li–O peak in the ∼2 Å range, which is the characteristic of lithium oxygen bond distances for DME–lithium interactions. Li2CO3 does not exhibit the double Li–O oxygen peak. The Li–O bonds from the carbonate are merged into one peak with the DME–O solvation atoms. The Li–O CO3 bond distances, in Fig. S-3, are longer than those observed for LiOH, so the double peak has become merged into one. The final SEI molecule that was investigated was Li2O. The Li–O RDF exhibits a double peak similar to LiOH unlike Li2CO3. In Fig. S-4, the Li–O bond distance of the Li2O molecule is ∼1.7 Å, which is more similar to that of LiOH than Li2CO3, which explains the double peak over one broader peak. As with all the other SEI components, the radial distribution functions demonstrate that the primary interactions of the SEI species are solvation of Li+ and O from the DME interaction. There are other well-defined peaks in Fig. 2, such as Li–C and Li–H, but given the greater distance from the SEI species in question, these peaks are secondary effects of the DME–oxygen coordination with the lithium ions.

FIG. 2.

Radial distribution functions for Li–X pairs in DME: (a) LiF, (b) LiOH, (c) Li2CO3, and (d) Li2O. For X, turquoise represents F, red represents O, black represents C, purple represents Li, and green represents H.

FIG. 2.

Radial distribution functions for Li–X pairs in DME: (a) LiF, (b) LiOH, (c) Li2CO3, and (d) Li2O. For X, turquoise represents F, red represents O, black represents C, purple represents Li, and green represents H.

Close modal

In order to examine all possible interactions in the system with radial distribution functions, one must investigate the other species (O–X, H–X, and C–X). For instance, the F–X RDF of the LiF system (Fig. S-5 B) shows only one clearly defined peak, which is due to the interactions between fluorine and lithium. This is the Li–F bond from the fluorine ion’s perspective and is located at the same position as the lithium’s Li–F RDF. In the Li–F RDF [Fig. 2(a)], there is a sharp peak for the Li–O interactions, which shows the presence of a DME solvation shell, but for the F–X RDF, the Li–F curve is the only sharp peak in the system. The corresponding O–X RDF (Fig. S-5 A) has strong peaks for carbon, hydrogen, and lithium less than 3 Å. Based on the location of the carbon and hydrogen peaks along with the geometry of the DME molecule, these peaks are simply the molecular structure of DME. The O–Li peak is due to DME solvating lithium with its oxygen atoms. There is a weak O–F peak present, which is not seen on the F–O RDF and is located at a distance of greater than 3 Å. This implies that the O–F RDF peak is due to fact that O is closely coordinating with Li, which is subsequently closely coordinating with fluorine. This demonstrates that Li–O solvation is the dominating type of solvation, and there is little evidence of solvation of the fluorine ion. The C–X and H–X RDFs show the structure of DME along with indirect coordination due to O–Li interactions. In the LiOH, Li2CO3, and Li2O systems (Figs. S-6/7/8), the RDFs show no clear solvation of their respective counterions much like in the LiF system. This further demonstrates that the dominant solvation is Li+ with the O atom from DME.

One major assumption that has been made thus far in this work is that the SEI molecules are always in contact ion pairs and closely coordinating, which will change the solution structure. However, the SEI may have other ionic association such as dissociation or solvent-separated ion pairs (SSIPs).

2. Solvent-separated SEI

We first investigated the solvation structure of SEI components in SSIP configurations by running AIMD simulations of each SEI molecule in a SSIP (separated by ∼7 Å) in the initial configuration of the simulation. The idea is to observe if the separation is stable or not and to characterize the behavior of the system. The initial configurations of each system can be found in Fig. S-9. For the systems with two lithium ions, both of the ions were separated from the anion, resulting in 2 Li+ and O2− and CO32−.

In the LiF simulation, we can observe that there is a sharp drop in the energy during the first picosecond of the simulation, which demonstrates that rearrangements are occurring in the system, lowering the overall energy. The geometry of the system was analyzed, as displayed in Fig. 3. The Li–F distance dropped from ∼8 Å into the range of 5–6 Å although the distance even dipped as low as 4 Å during the simulation. These distances are too large for the molecule to be considered recombined, so this simulation is stable during the simulation time as a solvent-separated lithium and fluoride ion. This subsequently changes both the Li and F radial distribution functions. The Li–X RDF no longer has a close strong fluorine peak, but instead the DME Li–O solvation peak and the Li–F peak are much further away corresponding to that distance. The F–H RDF shows new behavior, which is that there is a defined peak that is closer than the others, thus proving that the fluoride ion interacts with DME via hydrogen atoms. The lithium ion is solvated by Li–O interactions, and fluorine atoms are solvated via F–H interactions with DME.

FIG. 3.

LiF separated in DME. (a) Li–F distance as a function of time. Radial distribution functions for (b) Li–X and (c) F–X. For X, turquoise represents F, red represents O, black represents C, green represents H, and purple represents Li.

FIG. 3.

LiF separated in DME. (a) Li–F distance as a function of time. Radial distribution functions for (b) Li–X and (c) F–X. For X, turquoise represents F, red represents O, black represents C, green represents H, and purple represents Li.

Close modal

Unlike LiF, the LiOH separated system quickly recombined in less than 0.5 ps, as seen in Fig. S-10. This implies that the SSIP configuration of LiOH is not favorable but that would require further investigation and evaluation of the energies of the respective states to make confident conclusions. For the Li2CO3 separated system, we first examined the Li–C distances (Fig. S-11). This metric makes it easy to observe the distance between the two ions. The first lithium ion rapidly approaches the cluster and is recombined within ∼1 ps. The second lithium ion exhibits a different behavior where the recombination is almost instantaneous (<250 fs), but then, the distance increases before recombining at ∼1 ps. This dynamic behavior is clarified in Fig. S-12, which shows the lithium–oxygen distances. The first lithium ion very quickly begins coordinating with the second oxygen atom in the first moments of the simulation and begins to coordinate with the first oxygen atom before bouncing away and recombining with the third oxygen atom. For the second lithium ion, the lithium ion quickly coordinates with the first lithium ion and then bounces away from the third oxygen before settling with the second oxygen. The Li2O separated system also completely recombines like Li2CO3 and LiOH but on a longer timescale. Figure S-13 displays the Li–O distances of the system. The first lithium ion of the system quickly recombines in less than 1 ps, while the second lithium stays in the 4–5 Å range until 6 ps when it recombines in less than 250 more femtoseconds. There is a LiO anion that lives for ∼6 ps before Li2O is completely reformed. The second lithium ion does not experience a large drop in coordination, which demonstrates that the system was able to combine the ions and maintain solvation structure.

3. Ionic SEI

Another possible ionic association configuration is for the ions to be completely dissociated. The structure of the dissociated SEI in the pure solvent can be probed by running the dissociated SEI ions in their own simulation cells to represent a complete separation. Instead of simulating Li+ with F separated by the maximum distance within a cell, there would be two cells: one consisting of Li+ and another with F. This would be equivalent to infinite separation, while the separated system is ∼8 Å in the initial configuration. The infinite separation simulation is practical and sensical from a theoretical point of view but does introduce some problems from the computational side. In these simulations, the cells will have an overall net charge due to the ions, which introduces problems regarding the bleeding of charge from one box to another across the periodic boundary conditions. This condition forces the inclusion of computational artifacts to balance the charge in the cell, which are known to reduce the accuracy of results.32–34 

The following systems were run following the previously described methodology: Li+, F, OH, LiCO3, and CO32− (Fig. 4). We attempted to simulate O2− and LiO to describe the dissociation of the Li2O system, but the anions were highly unstable and would cause unrealistic behavior such as abstraction of hydrogen atoms from DME. These reactions were not observed for the solvent-separated system discussed before, so perhaps the background charge might have something to do with the instability. The simulation of Li+ in DME is not a novel simulation but is important to obtain the trajectories for comparison with other work in this report. The RDF for the lithium ion system is as expected with strong Li–O coordination and the structure of DME from the other RDFs (Fig. S-14).

FIG. 4.

Initial configurations of ions in solution: (a) Li+, (b) F, (c) OH, (d) LiCO3, and (e) CO32− in DME.

FIG. 4.

Initial configurations of ions in solution: (a) Li+, (b) F, (c) OH, (d) LiCO3, and (e) CO32− in DME.

Close modal

The F system provided further confirmation for solvation behavior that was previously seen in the LiF separated system. From the F–H RDF curve in Fig. S-15, there is a sharp hydrogen peak, which is closer than any other interactions in the system, thus providing further proof of the solvation of F with the hydrogen atoms of the DME. Since there are two different modes of interaction between the F and Li with DME, it is possible for one of the DME molecules to solvate both F and Li+ simultaneously. In OH with the DME system (Fig. S-16), their RDFs are very similar to the RDFs of the pure DME systems. There are peaks in the O–H and H–O RDFs of ∼1 Å, which correspond to the OH anion structure, but the remaining peaks are not significantly changed in behavior, which leads to the conclusion that OH does not interact with the DME in any significant way unlike F. The completely dissociated case, CO32−, also indicates no specific coordination with DME similar to OH. These two anions show no change in the structure of the solution to create a specific interaction compared to pure DME (S-18).

The LiCO3 system’s radial distribution functions (Fig. S-17) are representative of strong Li–O coordination due to the combined effects of the CO3 anion DME. There is no presence of carbonate solvation interactions in the other curves based on referencing C–X, O–X, and H–X to the pure DME RDFs. This provides further evidence that the primary solvation interaction between DME and the SEI species is with the oxygen atom in the DME molecules, but if there is no lithium, then exceptions can occur such as in the LiF system.

4. Clusters

The original model for the contact ion pair SEI in the liquid state was based on a simple assumption that an SEI molecule would be isolated from other molecules and with a full primary solvation shell. However, research into liquid structures of other ionic systems has demonstrated that there are other possibilities for liquid structures involving the presence of solvated clusters.35,36 This implies that the most favorable structure for some solvates is not isolated but instead interacting closely with other molecules.

We investigated this effect by creating systems with two of each contact ion pair SEI molecule in DME and simulated the systems using the previous methodology. Nothing has fundamentally changed about the systems other than the concentration of the SEI and the initial configuration. Each of the SEI molecules were placed close together to start the system in a cluster formation so that the system will be representative of an SEI cluster in DME instead of just two separate molecules (Fig. S-19).

Figure 5 displays the Li–X radial distribution functions for the SEI clusters in DME systems. The curves are quite similar to the previous cases with the exception of Li–Li interactions and some different Li–O peak behavior. For the LiF system, the Li–F RDF shows a large and narrow peak for Li–F, which indicates that the Li–F bonds of the Li–F molecules are the most common type of interaction. The presence of a clear Li–O peak demonstrates that Li–O coordination is still present in the cluster systems. Unlike in the single LiF case, there is a Li–Li peak, which is located at the distance of ∼3–4 Å, which shows the structure of the LiF cluster. This is the indirect coordination of one lithium with another via fluorine atoms. The F–X RDFs confirm these observations with the strong and close (<2 Å) F–Li peak and a F–F peak at ∼3–4 Å. There is no change in the type of observed solvation interactions between LiF and DME. Figures S-20 and S-21 show the analysis of the geometry of the LiF cluster in the simulation. There is a significant amount of rearrangement in the Li–Li and F–F distances in the first picosecond of the simulation, which shows that the initial configuration was not favorable. From the F–Li distances, the first lithium ion is equally shared between two fluorine atoms, whereas the second lithium is only close to one of the fluorine atoms. The shape of the cluster is linear, which corroborates the distance over time graphs. For the LiOH cluster system, the Li–O peak shows a broadness, which was a deviation from the lone LiOH RDF, which had a double peak attributed to DME solvation and Li–OH interactions. The LiOH cluster Li–O peak is more akin to the Li2CO3 in DME simulation, which is due to the Li–CO3 oxygen interactions being at the same length scale as the Li–DME interactions. For the LiOH cluster, the interactions are at the same lengths as in the lone simulation, but Li interacting with other OH in the system adds interactions at a longer length scale that broaden the peak. The geometry of the system can be explored in Figs. S-22/23. From the snapshot of the system, we can observe that the LiOH cluster is intact and is in a linear configuration much like the LiF system. One of the two Li ions is shared closely between the OH anions, and the other lithium ion is less bound.

FIG. 5.

Li–X radial distribution functions for clusters of 2 in DME: (a) LiF, (b) LiOH, (c) Li2CO3, and (d) Li2O. For X, turquoise represents F, red represents O, black represents C, purple represents Li, and green represents H.

FIG. 5.

Li–X radial distribution functions for clusters of 2 in DME: (a) LiF, (b) LiOH, (c) Li2CO3, and (d) Li2O. For X, turquoise represents F, red represents O, black represents C, purple represents Li, and green represents H.

Close modal

The radial distribution function for the Li2CO3 system demonstrates many similarities to the lone Li2CO3 system. The largest deviation in previous behavior arises from the Li–Li curve because there is an additional peak in the cluster simulation. The additional peak is located at a distance of ∼5 Å, and based on distance evaluations, it is found to be caused by the indirect coordination of lithium ions from one carbonate ion to another. The distance evaluations for Li–C can be found in Figs. S-24/25. We can observe that the first carbon atom is about equally close to all the lithium ions in the system in a range of ∼2.5–3.5 Å although these geometries are not as constant as seen in previous systems. The second carbon atom is also approximately equally close to all the lithium ions in the system, which implies that this cluster is not linear as the LiOH and LiF systems were determined to be. The snapshot of the system at equilibrium (Fig. S-26) shows a different outcome, which is a linear configuration for the Li2CO3 molecules, but they are interacting across the periodic boundary conditions, thus giving the equally shared lithium ion conclusion from the distances. The interaction of the Li2CO3 molecules across the boundary in this way is not desirable and would be more accurate to simulate the cluster in a larger box although it would be much more computationally challenging.

Since Li2O is the other larger SEI, we might expect to run into similar problems with the periodic boundary conditions and cell size, but from Fig. S-27, we can observe that the cluster is not in a linear configuration interacting across the cell boundaries but instead in a staggered configuration that provides sharing of two of the lithium ions between the oxygens from cluster to cluster. This is further corroborated by the distance plots in Figs. S-28 and S-29, which show that the first oxygen closely coordinates with the first three lithium ions, whereas the second oxygen atom only closely interacts with the last two. The third lithium ion is shared equally among the Li2O molecules. The Li–Li radial distribution function also shows this behavior with two defined peaks: one for lithium that is sharing the same oxygen and another for indirect coordination between the two Li2O species.

In general, the shapes of the Li–X RDFs are the same, which proves that the solvation structure of the SEI species did not fundamentally change. From the remaining RDFs (Figs. S-30/31/32/33), the O-RDF shows conventional Li–O solvation behavior, while the C–X and H–X RDFs primarily show the molecular structure of DME.

Although radial distribution functions give an informative view of the structure of a liquid, they are not capable of providing energetics and dynamic or time-dependent information. However, the solvation dynamics and energetics of the system are also required to obtain a full picture of the SEI dissolved into the electrolyte. In previous work, we created solvation structure maps to help analyze this information.37 Since the primary interaction between the Li ion and the solvent molecules is from Li–O, the primary solvation shell consists of a set of discrete structures based on the number and configuration of solvent molecules. The solvation structure maps for DME and 1,3-dioxolane (DOL) in Fig. 6 consist of unique structures (A and B). Structures are connected by the addition, loss, or rearrangement of solvent molecules, which are represented by the arrows. The energetics of the solvation/desolvation were determined with non-PBC implicit and explicit combined solvation models to predict the thermodynamically favorable structures. The advantage of the solvation structure maps is that the energetics of the primary solvation shells from the AIMD simulations and the energetics can be estimated. Since the SEI molecules have Li–O as the dominant solvation interaction, these solvation maps for Li+ are highly applicable.

FIG. 6.

(a) Solvation map of DME and lithium. The thermodynamically favored solvation structure diagram for Li+ ion with DME. Each letter represents a unique geometry of Li+ ion with DME molecules. DME molecules can be added in either their linear configuration (L) or nonlinear configuration (NL). Rearrangements between linear and nonlinear DME configuration may also occur. All letters surrounded by the green boxes are favored thermodynamically. All letters highlighted in green are favorable solvent additions or rearrangements. (b) Structures corresponding to the solvation map. Adapted from Kamphaus and Balbuena, J. Phys. Chem. C, 121(39), 21105-21117 (2017). Copyright 2017, American Chemical Society.

FIG. 6.

(a) Solvation map of DME and lithium. The thermodynamically favored solvation structure diagram for Li+ ion with DME. Each letter represents a unique geometry of Li+ ion with DME molecules. DME molecules can be added in either their linear configuration (L) or nonlinear configuration (NL). Rearrangements between linear and nonlinear DME configuration may also occur. All letters surrounded by the green boxes are favored thermodynamically. All letters highlighted in green are favorable solvent additions or rearrangements. (b) Structures corresponding to the solvation map. Adapted from Kamphaus and Balbuena, J. Phys. Chem. C, 121(39), 21105-21117 (2017). Copyright 2017, American Chemical Society.

Close modal

1. Solvation maps of ion pairs in solution

In order to characterize the structure of the solvation shells, we first focused on the O coordination number of the lithium ion in each system. This is the number of oxygen atoms that are within coordination distance (<3 Å) at a certain point in time. The average coordination number can be determined by integrating the Li–O RDF, but in this case, we investigated the instantaneous coordination number determined from the trajectories of the Li and O atoms in the system.

Figure 7 shows this coordination number frame by frame for LiF in the DME system. At first, Li starts out with zero close oxygen atoms before rapidly becoming coordinated with three oxygen atoms by 2 ps. Then, the system is stable for the rest of the simulation with some small amount of time spent at coordination numbers of 4 and 2. Given the nature of molecular dynamics, fluctuations will be present, so it is not surprising that the coordination number is not entirely constant. From Fig. 7, we can see that in the first 0.5 ps, DME 3 begins to coordinate followed in quick succession by an O atom from DME 6. Then, around 2 ps, the other O atom from DME 3 begins to coordinate. This is relatively stable until ∼5 ps where the other O atom from DME 6 begins to coordinate pushing the coordination number up to 4 for a short period of time. However, as the second O atom from DME 6 coordinates, it pushes the other O atom on that DME away. Eventually, the originally coordinating O atom from DME 6 retakes its place pushing other O away. Figure 7(b) shows the solvation structure over time for the LiF system. Although the graph looks similar to the coordination number, it is not exactly the same since this graph is based on the solvation structures and not number. The system quickly progresses from structure A to B and then C. This follows the thermodynamics as predicted by our previous publication for Li+. Then, the system progresses to the E structure, which is mostly stable for the rest of the simulation. According to the previous thermodynamic calculations, E is one of the two energy minimum solvation structures for Li (Table S-1). The other energy minimum solvation structure has a coordination number of 6, which is not expected for LiF since the presence of the F counterion will introduce steric or crowding effects on DME. The short-lived transitions from E to F and E to H violate our previous structure predictions based on those pathways having a positive change in free energy. This is likely due to the nature of AIMD not always being at a ground state energy due to the nature of the algorithm or the energetics of the solvation structures being different for reasons such as implicit vs explicit solvent and basis sets. The fact that structures H and F are short lived might indicate that the structures are not as energetically favorable as E in comparison. LiF with the DME system shows that Li–DME solvation via the O atoms is the dominant solvation interaction, which gives a coordination number of 3 with solvation structure E.

FIG. 7.

LiF in DME. (a) Coordinating the center of mass distance of DME and Li with the coordination number. The black points represent the instantaneous coordination number, and the solid color curves represent the center of mass distance between Li and DME molecules that interact closely with the lithium ion in the simulation time. For each DME, there are two possible oxygen atoms that can coordinate. (b) Solvation structures during simulation time corresponding to structures in Fig. 6. (c) Histogram of solvation structures in AIMD simulation corresponding to structures in Fig. 6.

FIG. 7.

LiF in DME. (a) Coordinating the center of mass distance of DME and Li with the coordination number. The black points represent the instantaneous coordination number, and the solid color curves represent the center of mass distance between Li and DME molecules that interact closely with the lithium ion in the simulation time. For each DME, there are two possible oxygen atoms that can coordinate. (b) Solvation structures during simulation time corresponding to structures in Fig. 6. (c) Histogram of solvation structures in AIMD simulation corresponding to structures in Fig. 6.

Close modal

The structural changes for LiF are shown in Fig. 8(a). In the LiOH system [Fig. 8(b)], the total coordination quickly progresses from 2 to 4 and is stable at 4 for long times (structure F) although the coordination number is 3 (structure E) at the end of the simulation. This coordination number (4) is higher than that of LiF, which demonstrates the effect the F anion has compared to the OH anion. With Li2CO3, both the Li+ ions have similar coordination in the 2–3 range [Fig. 8(c)]. From Fig. S-34, we can see that the second Li ion had more dynamic interaction with solvent, which gave the ion a coordination number of 3 occasionally. This demonstrates that both Li ions are not experiencing the same interactions in the simulation. This behavior is also observed in the Li2O simulations [Fig. 8(d)]. LiOH experiences the greatest single lithium ion coordination number, while Li2O has the highest overall coordination number when considering the sum of each lithium ion in the molecule.

FIG. 8.

Solvation map structures over time in (a) LiF, (b) LiOH, (c) Li2CO3, and (d) Li2O in DME AIMD simulations. Structures are named A, B, C, and D, as in Fig. 7.

FIG. 8.

Solvation map structures over time in (a) LiF, (b) LiOH, (c) Li2CO3, and (d) Li2O in DME AIMD simulations. Structures are named A, B, C, and D, as in Fig. 7.

Close modal

As observed in Fig. 8, the solvation structures over time do not follow the exact same trajectory in each system. For LiOH, the system progresses from C to D, which is an interesting development since D is a thermodynamically unfavorable structure. After staying in the D structure for ∼1 ps, the system transitions into the F structure, which is relatively stable until the system moves into the E structure in the last ps of simulation. The presence of D is unexpected according to the calculated energetics in the presence of implicit solvent, but, in general, with explicit solvent and PBC, the structures are less “bound,” which we hypothesized in our previous report. The most commonly observed solvation structure was F, which has a coordination number of 4 (Fig. 9/Fig. S-34). This indicates that LiOH is more easily solvated in comparison to LiF since the solvation is more extensive.

FIG. 9.

Histogram of solvation map structures over time in (a) LiF, (b) LiOH, (c) Li2CO3, and (d) Li2O in DME AIMD simulations.

FIG. 9.

Histogram of solvation map structures over time in (a) LiF, (b) LiOH, (c) Li2CO3, and (d) Li2O in DME AIMD simulations.

Close modal

In the Li2CO3 system, the first Li undergoes a quick transition from structure C to B, which involves losing DME, and then to G for most of the remainder of the simulation. The second Li starts with zero coordination and proceeds to the C solvation structure, which is long lived until it transitions to G/E and oscillates between them. The structural solvation differences between the two Li ions are due to the number of DME. G is a solvation structure with one DME, whereas C and E both involve two DME molecules. Along with the coordination numbers for each Li, this furthers the conclusion that although the two Li+ ions have similar solvation, one undergoes more interaction with DME than the other. For Li2O, we observe that the system transitions from no coordination in structure A to structure B, which is relatively stable until ∼6 ps when structure C becomes more favored. After that, the solvation structure oscillates between C and E. All these transitions are expected and predicted according to the Li+ solvation map energetics. The second Li ion starts with structure B and transitions to C, which is somewhat stable until the simulation moves into the G structure after ∼8 ps and is stable for the rest of the time. As with Li2CO3, the main difference between the two Li ions is that one has more interaction with DME vs the other. The first lithium’s most common structures after 6 ps are E and C both of which involve two DME as opposed to G, which only has one DME for the second Li ion.

Li+ solvation by DME’s O atoms is the primary type of solvation for all the SEI molecules. This type of interaction dictates the structure of the SEI in the liquid state. For Li2O and Li2CO3, the solvation was very similar to the unequal interaction between the two Li ions and the presence of G as a major stable solvation structure (Fig. 9). This type of solvation was different from that of LiF and LiOH, which favored larger structures such as E and F. The most likely explanation for this behavior is crowding or steric effects due to the counterion or the presence of two conflicting solvation shells for two lithium ions. Based on the previous energetic calculations, non-linear or doubly coordinated DME molecules have lower energies than their linear or singly coordinated counterparts if they have the same overall coordination.37 This likely explains why structure G would be more common in the more constrained geometry of the Li2O and Li2CO3 solutions.

2. Solvation maps for separated ions

Many of the solvent-separated molecules were not stable as separated ions for a significant time, making the time-dependent analysis of these systems the same as the original systems. For the ionic systems, Li+ and LiCO3 are the only systems, which contain lithium species that the Li–O coordination analysis can be used for.

Li+ in the DME system offers the best comparison to the solvation structure maps since there are no counterions present. The solvation structure analysis in Fig. 10 aligns with the previously calculated solvation map very well. We can observe that the lithium ion starts unsolvated (structure A) and quickly progresses from A → B → C → E. Then, the system progresses from E → F, which is unfavorable thermodynamically based on the previous calculations, but this transition has been noted many times in this work.37 Then, the system progresses from F to other structures but mostly oscillates between J and F until ∼30 ps where the system starts to oscillate between J and I. Interestingly, I is one of the predicted final solvation structures for Li+ based on providing the minimum energy. Other than the E → F transition, the solvation map provided good predictions for what was observed in the AIMD simulation. This indicates that the different levels of theory can be used to help understand other level of theory simulations and that the discrete solvation maps can be a powerful technique to probe solvation chemistry. For the LiCO3 system, even though there is only one lithium ion in this system, the total coordination number is not significantly higher than Li2CO3, which demonstrates the crowding effect caused by the carbonate ion alone. The most favorable observed solvation structure depends on the time during the simulation. Both C and G are present in much higher amounts than other structures, but they both have regimes of time where the structures are stable.

FIG. 10.

Li+ in DME. (a) Li solvation structure tracking frame by frame based on the previously determined solvation map for Li+ in Fig. 7. (b) Histogram of the corresponding solvation structures for Li.

FIG. 10.

Li+ in DME. (a) Li solvation structure tracking frame by frame based on the previously determined solvation map for Li+ in Fig. 7. (b) Histogram of the corresponding solvation structures for Li.

Close modal

The results for the time independent analysis suggest that the primary factor involved in the solvation structures of the system is the freedom of the lithium atom to interact with oxygen atoms of DME. In the original cases, the solvation structures follow pathways that lower their energy until the steric forces of the counterion make solvation not feasible. When examining the ionic systems, the solvation structures and trends were almost exactly the same except that the lithium ion had more ability to interact with DME. In the cluster simulations, the opposite occurred, where the increased number of lithium ions and counterions encouraged smaller solvation structures.

3. Solvation maps in clusters

We turn to the cluster simulations to explore if the interactions between the SEI molecules change the solvation structures. In the case of the LiF cluster, the two lithium ions are not symmetrically equivalent and this is the root cause of the differences in the coordination numbers. The second Li ion is shared between two fluorine atoms, which significantly increases the steric shielding or crowding of the ion to DME. The natural consequence is that the less shared lithium ion can interact more closely with DME molecules. The coordination number of the first lithium is ∼3–4, whereas the second and more shared lithium is ∼1–2 (Fig. S-35). The solvation structure analysis, Fig. S-36, shows that the first lithium atom progresses to the favored structure E and the second lithium ion’s solvation is mostly structure B, which is a result of the increased crowding of the ion. For lone LiF in the DME case, the favored solvation structure was found to be E, so the addition of another LiF does not change that outcome for the less shared ion. The LiOH simulation leads to similar conclusions as with LiF for the total Li–O coordination. Li that is less shared has a higher total coordination number than the more shared Li. The first lithium ion that has a higher coordination number favors D as the dominant solvation structure. This result is particularly interesting because D is an unfavorable structure, no matter the pathway.37 The other Li ion has B as the most common structure, which is similar to the more crowded Li ion in the LiF cluster simulation.

With the Li2CO3 and Li2O simulations, there are four lithium ions that may be solvated, which can present competitive behavior. In the Li2CO3 simulation, the first lithium ion has little to no coordination vs the third lithium atom having a steady coordination number of 2. We hypothesize that this is because there are only five DME molecules in the system but four lithium ions, which means that lithium ion solvation with DME will not be homogeneous. B and C are the dominant solvation structures, which are unexpected in an unrestricted system, but this simulation is highly crowded with carbonate ions. In the Li2O model, the coordination analysis shows similar behavior to the Li2CO3 cluster system where the coordination numbers are not equivalent and two of the lithium ions even have zero coordination for significant lengths of time. As with the carbonate system, the Li2O system has the same ratio of DME molecules to lithium ions, so there is not quite enough coordination to go around. The most stable solvation structures are B, C, and G (Fig. S-37). The presence of G is interesting since that structure is sterically more challenging to fit around the lithium ion in a crowded environment.

Throughout the discussion so far, the focus has been centered on characterizing the solvation structure of the SEI components in different possible configurations. This answers questions about what would the structure be if the SEI had dissociated or was dissolved as a molecule into the solution. The structural information and conclusions are vital to the overall description and characterization of SEI components that have dissolved. The ionic association and structural details have allowed for a fundamental understanding of different configurations of SEI components in the electrolyte. However, this raises the question of how likely or possible is it for those configurations to exist? This can be evaluated by investigating the energetics, which determines what is thermodynamically favorable. The structure and energetics are closely linked, and together, they provide a thorough characterization of the solution state.

The dissolution of SEI components from the solid state on an electrode into the electrolyte is fundamentally governed by the difference in energy between the two states.

Dissolution behavior is determined by what molecular environment provides a more favorable thermodynamic state. This is best expressed by the difference in free energy of the solid state vs the liquid state, and the energy difference can be directly linked to solvation properties. Fortunately, there are established methods and techniques that make determining the absolute free energy of a system possible. One of the hallmark outputs of a geometry optimization is the ground state energy of the system. This energy corresponds to the absolute energy of the system at 0 K. The absolute zero energy value is an important characteristic, but it is not the free energy of the system. The free energy of the same system would take into energetic contributions such as molecular translation, rotation, and vibration, which are not accounted for in the ground state energy.

For a solid system, translation and rotation are restricted due to the repeated molecular structure pattern of the solid state. This leaves vibration as the sole missing component of the free energy that needs to be accounted for. In order to include it, one must know the vibrational modes of the system in question. Thankfully, there are well-defined algorithms, which are included in most simulation packages that can determine the vibrational modes with relative ease. In the case of ground state optimizations, the derivatives of energy with respect to position were determined numerically in a post-optimization simulation, which determines the unique vibrational modes and their frequencies. The free energy of the system can be then determined from a combination of the ground state energy with the vibrational contributions by the following formula:

G=E0;ρ+i=13N6hνiehνikt1ehνikt+hνi2.

For the liquid system, the same approach could be used, but that would introduce some significant approximations. One is that the previous technique involves a ground state optimization first, which is challenging and questionably accurate to do for a liquid system due to the presence of many local minima. The previous approach also assumed that there were no contributions from translation and rotation, given the molecular structure of a solid, but that is not true for a liquid state. The free energy for the liquid state should have contributions from all three translation, rotation, and vibration degrees of freedom. To obtain these values, a different procedure must be used.

Instead of using a ground state optimization, ab initio molecular dynamics is used to find the equilibrium configuration of the liquid system. The energy of the system will decrease, but the system is not guaranteed to be in a minimum energy state at a particular point in time. Once a system has spent enough time at equilibrium to sample a significant amount of possible state space configurations, the velocity autocorrelation function can be determined. From the velocity autocorrelation function, there is a previously established procedure, which can estimate the translational and rotational components of the free energy for a liquid.27–30 For the vibrational components, instead of finding the vibrational modes through the hessian matrix, the velocity autocorrelation function can provide the vibrational density of states (VDOS). The VDOS is analogous to the difference between discrete orbital energy levels and the electronic density of states (DOS). The statistical thermodynamics equations from the discrete vibrational modes can be used for the VDOS by adapting the formula for continuous functions.

FIG. 11.

Li2O in DME. (a) Total VDOS. The y axis is the product of the VDOS by the speed of light. (b) VDOS per atom type. Black represents the total VDOS, blue represents the H VDOS, red represents the O VDOS, purple represents the Li VDOS, and green represents the C VDOS.

FIG. 11.

Li2O in DME. (a) Total VDOS. The y axis is the product of the VDOS by the speed of light. (b) VDOS per atom type. Black represents the total VDOS, blue represents the H VDOS, red represents the O VDOS, purple represents the Li VDOS, and green represents the C VDOS.

Close modal

Figure 11 shows an example of the VDOS for Li2O with the DME system. These plots show essentially vibrational spectra at different wavenumbers. The height of the peaks is related to the prevalence of that particular mode in the simulated sample. The width of the peaks is due to the close combination of different vibrational modes and a broadening of a particular mode’s frequencies due to interactions with other atoms. The numerical hessian procedure only provides a list of frequencies, which would be equivalent to vertical lines of unit thickness in these plots. From examining Fig. 11(b), we can observe the interactions contributing to each peak. Most of the higher wavenumber peaks are due to C–H vibrations at different frequencies. For instance, the peak in the 3000 wavenumber region has a small contribution of C and O and a large H contribution. This is likely due to C–O bond vibration, but since H atoms are bonded to C, the actual contribution of the VDOS is picked up by the H atoms. Most of the total VDOS arises from the O, C, and H VDOS and not the Li2O molecule. This is expected due to the much greater number of atoms in the system that are DME vs those that are Li2O. In fact, to obtain the “correct” energy for the system, the energy of the bulk DME must be subtracted out from that obtained from that of the SEI molecule with DME AIMD simulations. The corresponding DME simulations and analysis are discussed in the DME variation section of the supplementary material (page S-46). The VDOS for the LiF, LiOH, and Li2CO3 systems are displayed in Fig. S-38.

In order to compare the results of the different methods, we also applied the molecular dynamics version of the technique to the solid-state systems. The rotational and translational components of the free energy were still ignored, but the VDOS technique was used instead of discrete vibrational modes from VASP optimizations. Since AIMD was used to simulate the cells, the corresponding RDFs of the four SEI components can be found in Figs. S-39/40/41/42. The VDOS are also obtained as shown for the Li2O system in Fig. 12. The shape of the VDOS is noticeably different from the corresponding system in the liquid state shown in Fig. 11. The largest difference is the lack of DME vibrational modes that are present in the liquid system, which allows the Li2O vibrations to characterize the VDOS entirely. Another important observation is the shape of the graph in the low frequency regime (∼0 wavenumbers). The VDOS of the solid system has no magnitude at that point unlike the liquid system. This was discussed in previous publications about this technique and is attributed to the fluidicity of the system, which is related to the translational energy of the system.28,38 As expected, the solid system has no translational energy and the VDOS confirms that. The remaining VDOS graphs are shown in Fig. S-43.

FIG. 12.

Li2O solid state. (a) Total VDOS and (b) VDOS per atom type. Black represents the total VDOS, blue represents the H VDOS, red represents the O VDOS, purple represents the Li VDOS, and green represents the C VDOS.

FIG. 12.

Li2O solid state. (a) Total VDOS and (b) VDOS per atom type. Black represents the total VDOS, blue represents the H VDOS, red represents the O VDOS, purple represents the Li VDOS, and green represents the C VDOS.

Close modal

The comparisons of the free energies of solvation/dissolution are shown in Table I. The absolute free energies of the respective solid and liquid systems were calculated and can be found in the supplementary material, Tables S2–S6. The absolute free energies have no meaning except that can be used to calculate ΔG, which tells us the change in free energy between the two states. From the procedure, we can observe that all the tested systems have a positive change in free energy when the system goes from solid to liquid (Gliquid–Gsolid). It is important to note that Gsolid is based on the previous pristine bulk structures. However, as has been identified recently, when dissolution occurs, the energy of the solid state may change, which would subsequently modify the thermodynamics. The calculations reported here do not take this effect into account and should be studied further.39 Based on the way the calculation was defined, the results in Table I imply that the dissolution of an SEI molecule from the solid state is not thermodynamically preferred in any case. Comparing the systems to each other shows that LiOH is the most likely to dissolve followed by Li2CO3, LiF, and Li2O, respectively. There are many possible reasons for this trend, but, in general, if there are any particular chemical interactions that stabilize the energy in the liquid or solid state, then the energy will be favored in that direction. The solvation energy trend indicates that a SEI that consists of Li2O will be more stable regarding dissolution in the electrolyte than a Li2CO3 or LiOH counterpart. The solubility is directly related to the free energy of dissolution by the means of an equilibrium constant expression. The solubility is simply the equilibrium constant that corresponds to ΔGdissolution.

TABLE I.

ΔG of solvation/dissolution calculated as the difference between Gliquid and Gsolid.

CompoundΔGsolv/dissol (eV)
LiF 2.53 
LiOH 1.61 
Li2CO3 2.34 
Li23.38 
CompoundΔGsolv/dissol (eV)
LiF 2.53 
LiOH 1.61 
Li2CO3 2.34 
Li23.38 

It is also easy to test assumptions about the structure of the liquid state. As discussed earlier, the SEI molecules may be dissociated or clustered instead of freely solvating such as in the original simulations. Table II displays the results of the dissociation energy calculation from the associated VDOS (Fig. S-44). The dissociation energy is defined as the free energy of the reaction AB → A + B. In every case, the dissociation energy is positive, which means that the SEI molecules are energetically favored to stay in close-ion pair formation compared to separated ion pairs. LiF is more likely to dissociate than LiCO3, Li2CO3, and LiOH, respectively. This could be due to the solvation of the fluorine ion as observed earlier, but the LiCO3 fragment is solvated by the DME as well, and there is a difference of ∼1 eV. In order to determine a more accurate hypothesis, more simulations would need to be executed. It is also an important note that there are other possible structures the ions could be in, such as aggregate ion pairs and so on, which could be more energetically favorable. Since these results do not thermodynamically favor dissociation, they do not change the previous model for dissolution of a SEI molecule.

TABLE II.

Dissociation energy of different SEI molecules. ΔG was calculated as the difference between GA+B and GAB.

CompoundΔG (eV)
LiF 2.82 
LiOH 3.98 
Li2CO3 3.77 
LiCO3 3.33 
CompoundΔG (eV)
LiF 2.82 
LiOH 3.98 
Li2CO3 3.77 
LiCO3 3.33 

To compare these systems to their non-clustered counterparts, the same procedure can be used and the end result of the energy will be normalized by the number or concentration of SEI species. The VDOS for LiF, LiOH, Li2O, and Li2CO3 can be found in Fig. S-45. For the most part, the shapes and relative magnitudes of the VDOS peaks are the same as what was previously observed, which indicates that any additional inter-cluster vibrations are minimal compared to the molecular vibrational modes.

Cluster formation energies correspond to the reaction nA + A → (A)n+1. The calculated cluster formation energies in Table III demonstrate that it is energetically favorable in every case for the SEI molecules to remain isolated in solution instead of in small clusters. LiF is most likely to form a cluster followed by LiOH, Li2O, and finally Li2CO3. These results indicate that the previous solvation or dissolution model assumption that the SEI would be isolated from other SEI molecules in solution was valid. Clustering is not energetically favored and would not be expected to occur at least in this concentration level.

TABLE III.

Cluster formation energies of different SEI molecules. ΔG was calculated as the difference between GAn+1 and GnA+A.

CompoundΔG (eV)
LiF 1.20 
LiOH 2.44 
Li2CO3 3.98 
Li22.80 
CompoundΔG (eV)
LiF 1.20 
LiOH 2.44 
Li2CO3 3.98 
Li22.80 

All the initial configurations, temperature, and energy evolution are reported for liquid solutions in Figs. S46–S62 and for solid systems in Figs. S63–S67. The properties related to DME simulations are included in Figs. S68–S83. As explained in Sec. II, all our AIMD simulations used the PBE functional. van der Waals empirical corrections were tested for DME, and negligible effects on vibrational spectra were found. All the details are explained in the section of DME properties (see Fig. S74, Table S6, and discussions of the supplementary material). We acknowledge that a few of our systems may have been benefited by the introduction of dispersion corrections, but given the DME results (Fig. S74) for consistency, we did not include them for any of the simulations in this study.

Characterizing the SEI in order to understand its fundamental properties is an important task so that more intelligent design may be carried out. One aspect of the SEI that is not well investigated is the stability of a SEI component in the presence of solvent. Dissolution of the SEI is a possible chemical process that could have a significant impact on the SEI stability and structure, which would lead to poor passivation performance. By using previously existing methods in computational chemistry, we are able to determine the structure of LiF, Li2O, LiOH, and Li2CO3 molecules in solution phase and calculate dissolution energies, dissociation energies, and the thermodynamics of cluster formation.

We find that the SEI molecules interact with a pure DME electrolyte primarily via Li–O interactions. Solvation of the counterions is rare and only observed with F–H interactions in the LiF system. The solvation structures observed in the solution were well characterized following previous work using non-PBC simulations to determine discrete solvation structure maps. From the AIMD simulations, the vibrational density of states was calculated for each system, which allowed for the estimation of the free energies in the solid and liquid phases. This allowed for the calculation of the change in free energy of solvation. We found that in all cases, the dissolution of the SEI was not favored, but LiOH was most likely to dissolve followed by Li2CO3, LiF, and Li2O. Since the results may be offset due to inaccuracies of the methods/models, but the uncertainty about this point exists due to the lack of a direct experimental comparison, we conclude that out of the compared species, Li2O would form the most stable SEI with regard to solvent stability. Given the assumptions in the liquid state system structure, we also evaluated the dissociation energy and cluster formation energy to determine other possible lower energy structures. From the free energies of dissociation, all systems were favored to be in close-ion pair with LiF being the most likely molecule to become dissociated. With the cluster formation energies, we found that all SEI molecules were energetically preferred to stay isolated in solution. These characteristics provide valuable information regarding the properties and characterization of SEI in the solution state, which helps with further understanding and provides design principles for lithium–metal battery technology.

See the supplementary material for bond distance evolutions (Figs. S1–S4, S10–S13, S20–S25, and S27–S29), radial distribution functions (Figs. S5–S8, S14–S18, and S30–S33), initial configurations (Figs. S9 and S19), the snapshot of the cluster (Fig. S26), the coordination number and center of mass-Li distances (Figs. S34 and S35), solvation map evolution (Fig. S36), histograms of solvation structures (Fig. S37), VDOS per atom type (Figs. S38 and S43–S45), radial distribution functions of solid state (Figs. S39–S42), initial configurations of liquid solutions and system variable evolution (Figs. S46–S62) and of solid systems (Figs. S63–S67), DME properties (Figs. S68–S83), solvation energy for transitions between different states (Table S1), and absolute free energies of SEI components in DME (Tables S2–S6).

This work was supported by the U.S. Department of Energy’s Office of Energy Efficiency and Renewable Energy (EERE), as part of the Battery 500 Consortium, Award No. DE-EE0008210. Supercomputer resources from the Texas A&M University High Performance Computer Center and Texas Advanced Computing Center (TACC) are gratefully acknowledged.

The authors have no conflicts to disclose.

The data that support the findings of this study are available from the corresponding author upon reasonable request.

1.
P. G.
Bruce
,
S. A.
Freunberger
,
L. J.
Hardwick
, and
J.-M.
Tarascon
, “
Li–O2 and Li–S batteries with high energy storage
,”
Nat. Mater.
11
(
1
),
19
(
2012
).
2.
P.
Zhai
,
L.
Liu
,
X.
Gu
,
T.
Wang
, and
Y.
Gong
, “
Interface engineering for lithium metal anodes in liquid electrolyte
,”
Adv. Energy Mater.
10
,
2001257
(
2020
).
3.
N.
Khan
,
S.
Dilshad
,
R.
Khalid
,
A. R.
Kalair
, and
N.
Abas
, “
Review of energy storage and transportation of energy
,”
Energy Storage
1
(
3
),
e49
(
2019
).
4.
Y.
Fan
,
X.
Chen
,
D.
Legut
, and
Q.
Zhang
, “
Modeling and theoretical design of next-generation lithium metal batteries
,”
Energy Storage Mater.
16
,
169
193
(
2019
).
5.
X.
Chen
,
T.
Hou
,
K. A.
Persson
, and
Q.
Zhang
, “
Combining theory and experiment in lithium–sulfur batteries: Current progress and future perspectives
,”
Mater. Today
22
,
142
158
(
2019
).
6.
R.
Cao
,
W.
Xu
,
D.
Lv
,
J.
Xiao
, and
J.-G.
Zhang
, “
Anodes for rechargeable lithium‐sulfur batteries
,”
Adv. Energy Mater.
5
(
16
),
1402273
(
2015
).
7.
D.
Lin
,
Y.
Liu
, and
Y.
Cui
, “
Reviving the lithium metal anode for high-energy batteries
,”
Nat. Nanotechnol.
12
(
3
),
194
(
2017
).
8.
W.
Xu
,
J.
Wang
,
F.
Ding
,
X.
Chen
,
E.
Nasybulin
,
Y.
Zhang
, and
J.-G.
Zhang
, “
Lithium metal anodes for rechargeable batteries
,”
Energy Environ. Sci.
7
(
2
),
513
537
(
2014
).
9.
A.
Wang
,
S.
Kadam
,
H.
Li
,
S.
Shi
, and
Y.
Qi
, “
Review on modeling of the anode solid electrolyte interphase (SEI) for lithium-ion batteries
,”
npj Comput. Mater.
4
(
1
),
15
(
2018
).
10.
G.
Bieker
,
M.
Winter
, and
P.
Bieker
, “
Electrochemical in situ investigations of SEI and dendrite formation on the lithium metal anode
,”
Phys. Chem. Chem. Phys.
17
(
14
),
8670
8679
(
2015
).
11.
E. P.
Kamphaus
,
K.
Hight
,
M.
Dermott
, and
P. B.
Balbuena
, “
Model systems for screening and investigation of lithium metal electrode chemistry and dendrite formation
,”
Phys. Chem. Chem. Phys.
22
(
2
),
575
588
(
2020
).
12.
R.
Xu
,
X.-Q.
Zhang
,
X.-B.
Cheng
,
H.-J.
Peng
,
C.-Z.
Zhao
,
C.
Yan
, and
J.-Q.
Huang
, “
Artificial soft–rigid protective layer for dendrite‐free lithium metal anode
,”
Adv. Funct. Mater.
28
(
8
),
1705838
(
2018
).
13.
D.
Wang
,
W.
Zhang
,
W.
Zheng
,
X.
Cui
,
T.
Rojo
, and
Q.
Zhang
, “
Towards high‐safe lithium metal anodes: Suppressing lithium dendrites via tuning surface energy
,”
Adv. Sci.
4
(
1
),
1600168
(
2017
).
14.
F.
Hao
,
A.
Verma
, and
P. P.
Mukherjee
, “
Mechanistic insight into dendrite–SEI interactions for lithium metal electrodes
,”
J. Mater. Chem. A
6
(
40
),
19664
19671
(
2018
).
15.
L. A.
Selis
and
J. M.
Seminario
, “
Dendrite formation in Li-metal anodes: An atomistic molecular dynamics study
,”
RSC Adv.
9
(
48
),
27835
27848
(
2019
).
16.
K.
Tasaki
and
S. J.
Harris
, “
Computational study on the solubility of lithium salts formed on lithium ion battery negative electrode in organic solvents
,”
J. Phys. Chem. C
114
(
17
),
8076
8083
(
2010
).
17.
K.
Tasaki
,
A.
Goldberg
,
J.-J.
Lian
,
M.
Walker
,
A.
Timmons
, and
S. J.
Harris
, “
Solubility of lithium salts formed on the lithium-ion battery negative electrode surface in organic solvents
,”
J. Electrochem. Soc.
156
(
12
),
A1019
(
2009
).
18.
V. J.
Ovejas
and
A.
Cuadras
, “
Impedance characterization of an LCO-NMC/graphite cell: Ohmic conduction, SEI transport and charge-transfer phenomenon
,”
Batteries
4
(
3
),
43
(
2018
).
19.
Y.
Zhou
,
M.
Su
,
X.
Yu
,
Y.
Zhang
,
J.-G.
Wang
,
X.
Ren
,
R.
Cao
,
W.
Xu
,
D. R.
Baer
,
Y.
Du
,
O.
Borodin
,
Y.
Wang
,
X.-L.
Wang
,
K.
Xu
,
Z.
Xu
,
C.
Wang
, and
Z.
Zhu
, “
Real-time mass spectrometric characterization of the solid–electrolyte interphase of a lithium-ion battery
,”
Nat. Nanotechnol.
15
(
3
),
224
230
(
2020
).
20.
H.
Yildirim
,
A.
Kinaci
,
M. K. Y.
Chan
, and
J. P.
Greeley
, “
First-principles analysis of defect thermodynamics and ion transport in inorganic SEI compounds: LiF and NaF
,”
ACS Appl. Mater. Interfaces
7
(
34
),
18985
18996
(
2015
).
21.
F. A.
Soto
,
A.
Marzouk
,
F.
El-Mellouhi
, and
P. B.
Balbuena
, “
Understanding ionic diffusion through SEI components for lithium-ion and sodium-ion batteries: Insights from first-principles calculations
,”
Chem. Mater.
30
(
10
),
3315
3322
(
2018
).
22.
P.
Guan
,
L.
Liu
, and
X.
Lin
, “
Simulation and experiment on solid electrolyte interphase (SEI) morphology evolution and lithium-ion diffusion
,”
J. Electrochem. Soc.
162
(
9
),
A1798
(
2015
).
23.
R.
Mogensen
,
D.
Brandell
, and
R.
Younesi
, “
Solubility of the solid electrolyte interphase (SEI) in sodium ion batteries
,”
ACS Energy Lett.
1
(
6
),
1173
1178
(
2016
).
24.
E. P.
Kamphaus
,
M. S.
Angarita-Gomez
,
X.
Qin
,
M.
Shao
,
M.
Vijayakumar
, and
P. B.
Balbuena
, “
Role of inorganic surface layers on solid electrolyte interphase evolution at Li-metal anodes
,”
ACS Appl. Mater. Interfaces
11
,
31467
(
2019
).
25.
R.
Scott
,
M. P.
Allen
, and
D. J.
Tildesley
, “
Computer simulations of liquids
,”
Math. Comput.
57
,
442
444
(
1991
).
26.
Y.
Yamauchi
,
H.
Nakai
, and
Y.
Okada
, “
Short-time Fourier transform analysis of ab initio molecular dynamics simulation: Collision reaction between NH4+(NH3)2 and NH3
,”
J. Chem. Phys.
121
(
22
),
11098
11103
(
2004
).
27.
M.
Thomas
,
M.
Brehm
,
R.
Fligg
,
P.
Vöhringer
, and
B.
Kirchner
, “
Computing vibrational spectra from ab initio molecular dynamics
,”
Phys. Chem. Chem. Phys.
15
(
18
),
6608
6622
(
2013
).
28.
K.
Wendler
,
M.
Brehm
,
F.
Malberg
,
B.
Kirchner
, and
L.
Delle Site
, “
Short time dynamics of ionic liquids in AIMD-based power spectra
,”
J. Chem. Theory Comput.
8
(
5
),
1570
1579
(
2012
).
29.
T. A.
Pascal
,
K. H.
Wujcik
,
D. R.
Wang
,
N. P.
Balsara
, and
D.
Prendergast
, “
Thermodynamic origins of the solvent-dependent stability of lithium polysulfides from first principles
,”
Phys. Chem. Chem. Phys.
19
(
2
),
1441
1448
(
2017
).
30.
T. A.
Pascal
,
S.-T.
Lin
, and
W. A.
Goddard
 III
, “
Thermodynamics of liquids: Standard molar entropies and heat capacities of common solvents from 2PT molecular dynamics
,”
Phys. Chem. Chem. Phys.
13
(
1
),
169
181
(
2011
).
31.
D. A.
McQuarrie
,
Statistical Mechanics
(
Happer and Row
,
New York
,
1976
).
32.
K.
Leung
and
M.
Marsman
, “
Energies of ions in water and nanopores within density functional theory
,”
J. Chem. Phys.
127
(
15
),
154722
(
2007
).
33.
V. R.
Saunders
,
C.
Freyria-Fava
,
R.
Dovesi
,
L.
Salasco
, and
C.
Roetti
, “
On the electrostatic potential in crystalline systems where the charge density is expanded in Gaussian functions
,”
Mol. Phys.
77
(
4
),
629
665
(
1992
).
34.
A. Y.
Lozovoi
,
A.
Alavi
,
J.
Kohanoff
, and
R. M.
Lynden-Bell
, “
Ab initio simulation of charged slabs at constant chemical potential
,”
J. Chem. Phys.
115
(
4
),
1661
1669
(
2001
).
35.
N. N.
Rajput
,
V.
Murugesan
,
Y.
Shin
,
K. S.
Han
,
K. C.
Lau
,
J.
Chen
,
J.
Liu
,
L. A.
Curtiss
,
K. T.
Mueller
, and
K. A.
Persson
, “
Elucidating the solvation structure and dynamics of lithium polysulfides resulting from competitive salt and solvent interactions
,”
Chem. Mater.
29
(
8
),
3375
3379
(
2017
).
36.
A.
Andersen
,
N. N.
Rajput
,
K. S.
Han
,
H.
Pan
,
N.
Govind
,
K. A.
Persson
,
K. T.
Mueller
, and
V.
Murugesan
, “
Structure and dynamics of polysulfide clusters in a nonaqueous solvent mixture of 1,3-dioxolane and 1,2-dimethoxyethane
,”
Chem. Mater.
31
(
7
),
2308
2319
(
2019
).
37.
E. P.
Kamphaus
and
P. B.
Balbuena
, “
First-principles investigation of lithium polysulfide structure and behavior in solution
,”
J. Phys. Chem. C
121
(
39
),
21105
21117
(
2017
).
38.
W. C.
Swope
,
H. C.
Andersen
,
P. H.
Berens
, and
K. R.
Wilson
, “
A computer simulation method for the calculation of equilibrium constants for the formation of physical clusters of molecules: Application to small water clusters
,”
J. Chem. Phys.
76
(
1
),
637
649
(
1982
).
39.
A.
Wang
,
Z.
Zou
,
D.
Wang
,
Y.
Liu
,
Y.
Li
,
J.
Wu
,
M.
Avdeev
, and
S.
Shi
, “
Identifying chemical factors affecting reaction kinetics in Li-air battery via ab initio calculations and machine learning
,”
Energy Storage Mater.
35
,
595
601
(
2021
).

Supplementary Material