Transport and Mechanical Properties of Model ... - ACS Publications

Jun 30, 2017 - Electrochemistry Branch, Sensor and Electron Devices Directorate, Power and Energy Division, U.S. Army Research Laboratory,. Adelphi ...
0 downloads 0 Views 5MB Size
Article pubs.acs.org/JPCC

Li+ Transport and Mechanical Properties of Model Solid Electrolyte Interphases (SEI): Insight from Atomistic Molecular Dynamics Simulations Dmitry Bedrov,*,† Oleg Borodin,‡ and Justin B. Hooper† †

Department of Materials Science & Engineering, University of Utah, Salt Lake City, Utah 84112, United States Electrochemistry Branch, Sensor and Electron Devices Directorate, Power and Energy Division, U.S. Army Research Laboratory, Adelphi, Maryland 20783, United States



S Supporting Information *

ABSTRACT: A fundamental understanding of solid electrolyte interphase (SEI) properties is critical for enabling further improvement of lithium batteries and stabilizing the anode− electrolyte interface. Mechanical and transport properties of two model SEI components were investigated using molecular dynamics (MD) simulations and a hybrid MD-Monte Carlo (MC) scheme. A many-body polarizable force field (APPLE&P) was employed in all simulations. Elastic moduli and conductivity of model SEIs comprised of dilithium ethylene dicarbonate (Li2EDC) were compared with those comprised of dilithium butylene dicarbonate (Li 2BDC) over a wide temperature range. Both ordered and disordered materials were examined with the ordered materials showing higher conductivity in the conducting plane compared to conductivity of the disordered analogues. Li2BDC was found to exhibit softening and onset of anion mobility at lower temperatures compared to Li2EDC. At 120 °C and below, both SEI model compounds showed single ion conductor behavior. Ordered Li2EDC and Li2BDC phases had highly anisotropic mechanical properties, with the shear modulus of Li2BDC being systematically smaller than that for Li2EDC.

1. INTRODUCTION Progress in developing lighter, safer, and longer lasting energy storage devices critically depends not only on the design of mechanically stable and high capacity electrodes but also on the ability to engineer the interface between electrolyte and electrodes.1 Most electrolyte solvents currently used in Li-ion batteries are comprised of a mixture of linear and cyclic carbonates with additives. Carbonate-based electrolytes are not electrochemically stable at graphite or Si electrodes. They decompose and form a solid electrolyte interphase (SEI) on the negative electrode surface, protecting the electrolyte from further decomposition.2 Unfortunately, the SEI also contributes significantly to the interfacial resistance. The interfacial impedance is believed to have three primary components: (a) resistance due to ion transport through the SEI, (b) resistance due to Li+ desolvation from electrolyte upon insertion into the SEI, and (c) resistance due to the Li+ desolvation from the SEI and intercalation into the anode. It is difficult to separate contributions of those processes to the overall resistance, and hence identify the rate-limiting step in lithium transport. Most experimental efforts have focused on estimating the activation energy for the overall resistance for Li+ transport from bulk electrolyte through the SEI into the graphite electrode, which was found to be on the order of 0.6−0.7 eV when linear carbonate solvent molecules are present in the electrolyte.3−6 © 2017 American Chemical Society

Given the importance of the SEI in lithium batteries, improved molecular level understanding of SEI structure and Li+ transport can help to design novel solvents and/or additives that would form SEIs with enhanced Li+ mobility and yet will be efficient in both screening electrolyte from direct contact with electrode (i.e., protecting from reduction decomposition) as well as possessing desired mechanical characteristics. There have been a number of experimental studies focusing on understanding of the SEI structure and the mechanisms of its formation.1−6 In general, the SEI on graphite, Si, and Li metal anodes in contact with carbonate-based electrolytes is thought to consist of both inner and outer layers. The inner layer is primarily comprised of the doubly reduced compounds, such as Li2CO3, Li2O, and LiF, while the outer SEI layer is thought to be comprised largely of alkyl dicarbonate species7−9 such as dilithium ethylene dicarbonate (CH2OCO2Li)2 (or Li2EDC) and dilithium butylene dicarbonate (CH2CH2OCO2Li)2 (or Li2BDC). In electrolytes containing ethylene carbonate (EC) mixed with linear carbonates, it is believed that EC reduction products are the ones which primarily participate during SEI formation. This could be Received: May 4, 2017 Revised: June 27, 2017 Published: June 30, 2017 16098

DOI: 10.1021/acs.jpcc.7b04247 J. Phys. Chem. C 2017, 121, 16098−16109

Article

The Journal of Physical Chemistry C attributed to the preferential Li+ solvation by cyclic carbonates vs linear carbonates in bulk electrolytes10−15 and at electrode interfaces,16 together with the higher reduction potential for EC compared to the DMC in the cis−cis conformation.17 Li2EDC and Li2BDC have been suggested as major components found in the outer part of the SEI.2,18 FTIR and XPS studies of SEI films formed on lithium metal in EC-based electrolytes reveal the presence of alkyl carbonates,7,19 while FTIR shows the presence of alkyl carbonates in the SEI formed from EC/diethyl carbonate electrolytes on graphite.9 Hydrolysis of alkyl carbonate species formed from reduction of EC on noble metal electrodes leads to the formation of ethylene glycol, which has been taken as evidence for the formation of Li2EDC.20 NMR studies of the surface species formed on graphite in EC electrolytes resulted in spectra consistent with that of alkyl carbonates, and in particular the spectrum of Li 2 EDC. 21 A combination of NMR, XPS, and FTIR examination of cells cycled with LiPF6 in EC or EC/ethyl methyl carbonate (EMC) blends showed that the graphite SEI is ∼50 nm thick after the first full lithiation cycle, and predominantly contains lithium ethylene dicarbonate and LiF.22 A combined TEM/FTIR study of the SEI formed on graphite in EC/LiClO4 electrolyte revealed the formation of alkyl carbonates with chemical compositions of O/C ratio = 1.5 which is consistent with Li2EDC as well as with O/C ratio = 1.0 which is consistent with Li2BDC.23 Michan et al. recently used selective 13C labeling to detect decomposition products of ECand DMC-based electrolytes doped with LiPF6 salts.24 They also reported EC decomposition products being present in higher concentrations dominated by oligomer species with lithium dicarbonates, lithium fluoride, and other lithium carbonate products. Interestingly, signatures of both Li2EDC and Li2BDC were observed, presumably in the outer SEI because alkyl carbonates such as Li2EDC are thermodynamically unstable in contact with lithium metal as concluded from DFT calculations.25 Several classical and ab initio simulation works have provided much needed molecular scale insight into the initial stages of SEI formation.26−29 Molecular dynamics (MD) simulations using reactive force fields (ReaxFF) also have demonstrated that various alkyl carbonates can form from singly reduced electrolyte radicals in the bulk electrolyte and, hence, comprise the outer SEI layer at the anode.30,31 So far, there have been a very limited number of experimental and simulation investigations focused on characterizing Li+ transport in alkylcarbonate based SEIs. The temperature dependence of Li+ coordination and mobility in amorphous Li2EDC and lithium methyl carbonate was probed in MD simulations using a polarizable force field.32 This study showed that the charge transport in the system was primarily due to Li+ mobility, and extrapolation of predicted conductivities to room temperatures gave a value of ∼10−7 S/cm. Subsequent validation of these conductivities with recent experimental measurements as well as a detailed comparison with binding and conformational energies of model compounds/clusters predicted by high-level quantum chemistry (QC) calculations revealed that the original force field required some revisions. Recently, this polarizable force field has been reparameterized against higher-level QC calculations and fit to accurately reproduce conformational properties of alkylcarbonates, electrostatic field around these molecules in a vacuum, molecular dipole moments, and binding energies between alkylcarbonate anions and Li+.33 Using the revised force field, the amorphous and ordered (smectic-like

structure) phases of Li2EDC have been investigated as a function of temperature. Properties obtained from these MD simulations showed a noticeably slower dynamics of Li+ compared to simulations with the previous force field, and extrapolated values of conductivities (∼10−9 S/cm at room temperature) were found to be in good agreement with experimentally measured conductivity.33 Also, the Li2EDC− and Li2BDC−electrolyte interfaces were examined in MD simulations, suggesting that the activation energy for Li+ desolvation is smaller than the activation energy for Li+ diffusion in these SEI compounds.34 Jorn et al.35 performed ambitious MD simulations that not only encompassed both graphite−model SEI and electrolyte in one simulation cell but also applied external potential between electrodes. These simulations provided a molecular level insight into the species present at these interfaces and potential clues to the effect of the SEI on transport. Hybrid MD/Monte Carlo (MC) reaction simulations by Takenaka et al.36 suggested that a dense, ECbased SEI film can protect electrolyte from reduction at the electrode surface but at the same time provides cavities of sufficient size for Li+ passage. Single et al.37 elegantly connected SEI porosity with the interfacial impedance and its growth. Properties of the inner SEI compounds such as LiF and Li2CO3 were examined using DFT identifying the mechanism for Li+ diffusion through these phases38−42 including the voltage dependence43 and the space charge layer effect in the LiF− Li2CO3 composites on the overall conductivity of multicomponent SEI.44 The Li metal/Li2CO3|electrolyte interface was recently examined using the density functional tight binding approach, providing energetic estimates for the Li+ transport through these interfaces.45 Mechanical properties of the SEI impact its stability and ability to suppress the growth of Li dendrites during cycling, yet these properties have not been sufficiently investigated.46 In this paper, using atomistic MD simulations, we predict mechanical properties as a function of temperature and the mechanism of Li motion inside amorphous and ordered phases of alkyl-carbonate-based SEIs using Li2EDC and/or Li2BDC as model compounds.

2. SIMULATION METHODOLOGY AND SYSTEM DETAILS In this work, we investigated two types of SEI model systems: amorphous and ordered phases. For the former, bulk Li2EDC and Li2BDC as well as a 50/50 Li2EDC/Li2BDC mixture were investigated with each simulation cell containing 256 dilithium alkyl dicarbonates. Initially, each system was set up in a liquid phase and equilibrated for 5 ns at the highest temperature (700 K) and atmospheric pressure using NPT ensemble simulations. Subsequently, the resulting configurations were cooled down to the desired temperatures (600, 500, 450, and 393 K), and equilibration simulations over 5−150 ns in the NPT ensemble were conducted (see Table S1 in the Supporting Information for the summary of equilibration and production run length). Alkyl carbonates form ordered/crystal phases at room temperature. The X-ray study of Li2EDC reported lattice parameters that are consistent with an orthorhombic cell, while the lithium ethylene carbonate lattice was monoclinic.47 The crystal structure, however, has not been resolved. Considering the elongated shape of EDC and BDC molecules and the relatively flexible alkyl linkage between carbonate groups, it is reasonable to assume that a layered (smectic-like) ordering can be energetically favorable at low temperatures. Therefore, in 16099

DOI: 10.1021/acs.jpcc.7b04247 J. Phys. Chem. C 2017, 121, 16098−16109

Article

The Journal of Physical Chemistry C

3. RESULTS AND DISCUSSION 3.1. Structural Properties. We begin our analysis by examining the structural properties of these model SEIs as a function of temperature and composition. Density. Figure 2 compares densities obtained for pure Li2EDC and Li2BDC amorphous melts and ordered phases. As

addition to amorphous phases, we have investigated layered structures comprised of Li2EDC or Li2BDC. In these systems, the initial configurations were generated by shifting a Li2EDC (or Li2BDC) complex in the all-trans conformation of the alkyl segment along the x-, y-, and z-directions to generate a threelayer system in an orthorhombic cell in which alkyl carbonate molecules bridge layers of Li+. Within the layers, there are no long-range structural correlations between molecules; instead, a liquid-like order is observed. The snapshots of amorphous and ordered systems are shown in Figure 1.

Figure 2. Density as a function of temperature obtained from MD simulations of Li2BDC and Li2EDC amorphous and ordered SEIs.

expected, the figure shows that the Li2EDC melt is about 8−9% denser than the Li2BDC melt. At elevated temperatures, the density dependence is almost linear. However, as the temperature drops below 500 K, we observe a reduction in the slope (leveling off) of the density temperature dependence for amorphous systems. This is likely an indication of approaching the glass transition temperature. The ordered phases were not stable at temperatures above 500 K and therefore were only investigated in the 393−500 K temperature range. Densities for these systems are about 1.5−3% higher than the corresponding values in amorphous phases. The ordered systems also show a noticeably weaker (compared to amorphous phases) temperature dependence in the 393−500 K range. Radial Distribution Functions. Next, we analyze the Li+ cation coordination by examining the radial distribution functions (RDFs) g(r). Figure 3 shows g(r) between Li+ and ether (Oe) or carbonyl (Oc) oxygen atoms of BDC as well as Li+−Li+ g(r) obtained from simulations of Li2BDC systems. Several observations can be made. First, as expected, the sharp first peak in the gLi−Oc(r) indicates a strong preference for the

Figure 1. Snapshots of (a) the chemical structure of lithium alkyl dicarbonates, (b) the amorphous phase of Li2BDC, and (c) the ordered phase of Li2BDC.

Simulations have been conducted using the revised version of the polarizable APPLE&P force field discussed in ref 33. The cutoff for nonbonded van der Waals interactions was 12 Å. Long-range electrostatic interactions (charge−charge and charge−induced dipole) were treated using Ewald summation, while induced dipole−induced dipole interactions were calculated directly within the cutoff using a reaction field approach beyond the specified cutoff radius. Simulations were conducted using a multiple time step approach,48 with a 0.5 fs time step for integrating valence degrees of freedom (bonds, bends, torsions, and out-of-plane deformations), a 1.5 fs time step for treating short-range (within cutoff radius of 7 Å) nonbonded interactions, and a 3 fs time step for integration of remaining nonbonded interactions within cutoff radius, calculation of induced dipoles, and the reciprocal part of the Ewald summation for electrostatic interactions. Production runs ranged between 15 and 400 ns depending on composition and temperature. Long (>100 ns) simulation runs were necessary at temperatures below 500 K to ensure proper equilibration. Often in the literature, simulations of similar dynamically slow systems are conducted over much shorter trajectory lengths and are used to extract dynamical correlations/properties.49 Failure to access the necessary length of trajectories in such simulations would likely result in significant overestimation of the lithium transport and ionic conductivity, as the diffusive regime has not been reached. Finally, in addition to standard MD simulations described above, we also conducted combined MD−MC simulations to sample the mechanical properties of investigated SEIs. The protocol of these simulations is discussed in section 3.3.

Figure 3. Radial distribution functions (RDFs) obtained from MD simulations of amorphous Li2BDC at different temperatures. Also shown for comparison are the RDFs for the ordered Li2BDC at 500 K. 16100

DOI: 10.1021/acs.jpcc.7b04247 J. Phys. Chem. C 2017, 121, 16098−16109

Article

The Journal of Physical Chemistry C Li+ to interact with carbonyl oxygen, which is not surprising, since this oxygen atom has a highly negative partial atomic charge (∼ −0.8). The peak position is at about 2.0 Å, while the position of the first minimum at 2.8 Å can be used to define the first coordination shell of Li+. This first coordination shell also contains some ether oxygen atoms, as indicated by a much smaller peak in the gLi−Oe(r). For this pair correlation, the second peak, located at ∼4.0 Å, is slightly higher than the first one. The Li+−Li+ pair correlation shows a peak at 3 Å that is similar to the Li−Li distance in lithium metal and LixS melts.50 Figure 3 compares gLi−x(r) obtained from Li2BDC melt at several temperatures in the 450−700 K range as well as the Li2BDC ordered phase at 500 K. The shape of the pair correlations is remarkably insensitive to the temperature and the phase. Moreover, there are only slight changes in the gLi−x(r) values at the corresponding maxima, indicating that the local coordination of a Li+ is mostly independent of temperature or the ordering of alkyl carbonates. Analysis of gLi−x(r) obtained from Li2EDC melts (not shown) indicated the same features in pair correlation functions but with the magnitudes of the maxima about 20−30% lower than those shown in Figure 3. The latter is due to effectively higher concentration (per unit volume) of the considered atoms (Li, Oc, Oe) in Li2EDC compared to Li2BDC and therefore a different normalization for the pair correlation functions. Li+ Coordination and Distribution. Perhaps a better way to characterize the local structure around Li+ is to compare coordination numbers (CN) by various atoms as a function of the size of the surrounding shell. Figure 4 shows cumulative

coordination shell, the distribution of instantaneous local environments can range from one to eight Oc atoms. This is illustrated in Figure 5, which shows that most of the ions have

Figure 5. Distribution of an instantaneous number of Oc atoms coordinating Li+ (within 2.8 Å) as well as the number of BDC molecules those atoms belong to as obtained from MD simulations of the amorphous Li2BDC system at 450 K.

between three and five neighboring oxygen atoms. When we analyze how many BDC molecules those Oc atoms belong to, it is clear the majority of lithium cations have four different BDCs participating in their coordination (also shown in Figure 5). Similar distributions were obtained in Li2EDC systems as well as in the ordered SEI phases. Finally, Figure 6 illustrates the distribution of Li+ ions in the model amorphous and ordered SEIs. In this figure, the anion

Figure 4. Coordination numbers of Li+ obtained from MD simulations of amorphous SEIs at 500 K.

CN as a function of the shell radius. In Li2BDC melts, within the first coordination shell by oxygen atoms (r < 2.8 Å), each Li+ has about 4.3 Oc and 0.6 Oe atoms, resulting in about 4.9 total oxygen atoms. In this shell, there are also on average 0.7 Li+, while if the shell is extended to 4.2 Å (position of the first minimum in the gLi−Li(r)) then each Li+ has 3.8 cations in the vicinity, confirming the mentioned above tendency for Li+ to form multi-ion complexes/aggregates that include alkylcarbonate divalent anions. Comparison of the Li+ CN in Li2EDC melts showed very similar values, with CN(Oe) increasing to 0.85 and CN(Oc) and CN(Li+) reducing to 4.2 and 0.6, respectively, within the 2.8 Å coordination shell. This further shows that the local Li+ coordination shell is also independent of the composition of alkyl carbonate based SEI. A more detailed analysis of the Li+ coordination reveals that while on average each ion has about 4.3 Oc atoms in its first

Figure 6. Snapshots highlighting the Li+ distribution in the ordered Li2BDC and Li2EDC and the amorphous phases of Li2BDC at 393 K.

molecules are shown as wireframe and transparent, while the Li+ cations are shown as large spheres in order to highlight their arrangements. In the ordered phase, Li+ form 2D layers consistent with the smectic-like ordering of alkyl carbonates. This layering is well-defined in the Li2BDC, while in the Li2EDC the layering structure is less pronounced with a number of Li+ situated between the layers. The top view of one such layer in Li2BDC indicates that the distribution of Li+ in the layers is not homogeneous and represents a percolating 2D 16101

DOI: 10.1021/acs.jpcc.7b04247 J. Phys. Chem. C 2017, 121, 16098−16109

Article

The Journal of Physical Chemistry C

Figure 7. (a) Mean squared displacements of Li+ as obtained from MD simulations of amorphous Li2BDC systems at various temperatures.(b) Master curve obtained using the time−temperature superposition principle (see text).

Figure 8. Non-Gaussianity parameter (defined by eq 1) for amorphous Li2BDC at several temperatures and Li2EDC at 500 K as a function of time (a) and average Li+ MSD (b). The Li2EDC results are from ref 33.

network of regions with increased concentration of Li+. As has been shown by the strong first peak in the Li+−Li+ RDF in Figure 3, in these environments which are highly screened from electrostatic interaction environments the Li+ cations have no problem approaching each other within 3 Å or less. In the amorphous phase, the “clustering” of cations is similar to the ordered phase. However, the percolating network of regions of high Li+ concentrations is now three-dimensional, as shown in Figure 6. 3.2. Dynamical Properties. Next, we focus on the analysis of dynamic and transport correlations in these systems, primarily concentrating on understanding the Li+ transport mechanisms. Mean Squared Displacements. We begin with the analysis of mean squared displacements (MSDs) of Li+ ions obtained from simulations of the Li2BDC systems, as shown in Figure 7a. At higher temperatures, Li+ shows large displacements on the time scale of a few nanoseconds. The mobility of BDC anions is about a factor of 8 smaller (see Figure S1 in the Supporting Information) than that of Li+. Similar behavior was observed in the Li2EDC systems. At elevated temperatures (550−700 K), the MSD(t) are clearly reaching the expected linear dependence and hence the standard Einstein relation can be used to reliably extract self-diffusion coefficients (Di) of the species. However, at lower temperatures (400 K) compared to amorphous Li2BDC SEI, whose glass transition appears to be well below room temperature. Therefore, it brings a question of whether the difference in the glass transition would influence Li ion transport? Examination of the temperature dependence for Li+ diffusion in Figure 9 shows a more complex picture. At higher temperatures, the Li+ diffusion is about an order of magnitude higher than diffusion of the anion. In both systems, the diffusion of Li+ at elevated temperatures qualitatively follows the temperature dependence of anions; i.e., it shows an increasing rate of slowing down with decreasing temperature and the cations in the Li2BDC are becoming noticeably more mobile than the cations in the Li2EDC as the temperature goes down. However, as the diffusion of EDC and BDC anions drops to the range of ∼10−15 m2/s and below, a decoupling of Li+ and anion diffusion is observed. Below 450 K for Li2BDC, the VF-type dependence of diffusion of Li+ transitions to a more Arrhenius-like dependence (i.e., linear in the log(D) vs 1/ T plane) with an activation energy of about 0.6 eV, in agreement with the previously reported51 activation energy for the lithium transport in Li2EDC of 0.64 eV below 450 K for both amorphous and ordered SEI model compounds. Interestingly, another common SEI compound LiF when structurally disordered on the silica surface was also reported to have a similar activation energy of 0.55−0.65 eV, slightly depending on the thickness of the LiF layer.52 Somewhat higher activation energies (0.77 eV) were reported for the Li+ transport in Li2CO3 at 0.23 V and in the bulk,44 further indicating that many SEI components have similar activation energies. Moreover, within the accuracy of our simulations, the diffusion of Li+ in amorphous SEIs at lower temperatures seems to be independent of whether it has EDC or BDC as a counterion. Once the amorphous matrix (formed by anions) is frozen, the mobility of Li+ ions is determined only by their ability to move along the previously discussed percolating regions of high Li+ concentrations. Since the structure of those regions/percolating networks is very similar in the EDC- and BDC-based amorphous SEIs, the diffusion of Li+ is also very similar. Also, as discussed above, while the Li+ is on average coordinated with 4.3 oxygen atoms from anions, for most of the Li+ ions, these oxygen atoms come from four different anions; i.e., each anion is contributing about one oxygen atom to one Li+. Therefore, the energy barrier for a Li+ to exchange a coordinating oxygen from one anion to a coordinating oxygen from another is relatively small. This facilitates the hopping of Li+ along the percolating “channels” shown in Figure 6 and leads to the observed (in Figure 9) decoupling of anion and cation motion. In Figure 10, we compare Li+ diffusion in the amorphous and ordered systems. For the ordered systems, the simulations were limited to the 393−500 K range because at higher temperatures the lamellar-like order was not stable and has “melted” on the time scale of our simulations. Figure 10 illustrates that the diffusion in the ordered phase is higher than in the amorphous phase and that in this phase Li+ diffusion does depend on the type of anion forming the ordered structure. As we can see, the SEI containing BDC anion allows noticeably faster diffusion than the one with EDC anion, although the slopes (hence the activation energy) of the temperature dependence are very

becomes Gaussian on this time scale. As the temperature decreases, the peak in α(t) both grows and moves to longer times. At 450 K, the peak has a value of about 2.5 and is shifted to time scales of several nanoseconds, indicating a significant dynamic heterogeneity of Li+ motion on this scale. The shifting of the peak to longer times with decreasing temperature is due to slowing down of the overall dynamics in the system and the intrinsically longer time scales necessary for the ion to break out of its local “cage”. This can be clearly seen from Figure 8b where the same non-Gaussianity parameter is plotted as a function of average MSD(t) of Li+ ions at given t. In this case, the peaks for all temperatures are located between 1.0 and 2.0 Å2 which defines the displacements that are necessary for Li+ ions to escape their local cages. Finally, at 500 K, we also compare the α(t) and α(MSD(t)) for Li2BDC and Li2EDC systems. Consistent with Li2EDC being intrinsically slower compared to Li2BDC, while the locations of the maxima of the non-Gaussianity parameter (i.e., the cage size) are almost the same, the magnitude of the maxima (i.e., the extent of dynamic heterogeneity) is almost a factor of 3 larger in the Li2EDC system at 500 K. Self-Diffusion Coefficients. Figure 9 shows the self-diffusion coefficients for each component extracted from MD

Figure 9. Diffusion coefficients of Li+ and alkyl dicarbonates obtained from MD simulations of amorphous Li2EDC and Li2BDC systems. The dashed lines illustrate divergent Vogel−Fulcher-type behavior for EDC and BDC anions, while the dotted green line highlights the implied Arrhenius dependence for Li+.

simulations of amorphous Li2EDC and Li2BDC SEIs obtained from the analysis of MD simulations as discussed above. Several observations can be drawn from these data: First, the temperature dependence of diffusion coefficients of EDC and BDC anions is noticeably different. At 700 K, DEDC and DBDC are very similar. However, as the temperature decreases to 500 K, significant differences in diffusion can be observed. At 500 K, DEDC is about 20 times smaller than DBDC. As the temperature continues to lower, we expect those differences to increase even further, with the increasing rate of anion slowing down pointing to the behavior consistent with approaching a glass transition. To further illustrate this behavior, we have fit the inverse of diffusion coefficients with the Vogel−Fulcher (VF) equation log(1/D) = A + B /(T − T0)

(2)

which is typically used to describe the divergence of relaxation times in systems approaching the glass transition. Figure 9 shows that the obtained VF fit nicely describes the available data for both anions. The T0 value, which is typically somewhat 16103

DOI: 10.1021/acs.jpcc.7b04247 J. Phys. Chem. C 2017, 121, 16098−16109

Article

The Journal of Physical Chemistry C

Kohlrausch−Williams−Watts (KWW) stretched exponential function53 Hij(t )Hij(0) = A exp[−(t /τKWW )β ]

(3)

The stretching exponent β varies between 0.47 and 0.56, indicating a relatively pronounced dynamic heterogeneity of this process. In Figure 11b, we compare the residence autocorrelation functions for Li2BDC and Li2EDC amorphous melts at 500 K. Consistent with the diffusion data shown in Figure 9, at this temperature, the relaxation of the autocorrelation function in amorphous Li2BDC is about a factor of 8 faster than that in amorphous Li2EDC. However, as the temperature decreased further, the decay of autocorrelation functions becomes comparable (not shown) consistent with observed diffusion of Li+ at lower temperatures (Figure 9). 3.3. Mechanical Properties. In this section, we analyze mechanical properties of the amorphous and ordered Li2BDC and Li2EDC at the temperatures where anion does not significantly diffuse and materials could be considered to be in the essentially glass or smectic crystal phases. Methodology. The mechanical properties of a simulated system may be obtained by characterizing the fluctuations of dimensions of a fully flexible simulation cell,54 with the fluctuations providing a spontaneous expression of instantaneous strain for a simulation cell, which can be related to the compliance tensor of the material at equilibrium as

Figure 10. The Li+ cation diffusion coefficients for ordered and amorphous SEI model compounds Li2BDC and Li2EDC from MD simulations. The Li2EDC results are from ref 33.

similar. Note that, in the ordered phase, the motion of Li+ ions is basically confined to a 2D diffusion within the layer. In the Li2EDC, the spacing of layers is 8.0 Å which is close enough to allow occasional transfers of Li+ from one layer to another (as can be seen from the snapshot in Figure 6, where several Li+ are located in between the main layers). In the ordered Li2BDC, the spacing of Li+ layers is 10.0 Å and no transitions of Li+ across the layers were observed (also consistent with the snapshot in Figure 6 showing well-defined layers for this system). Li+ Residence Time. Finally, we characterize the residence time for Li+ to be coordinated by the carbonate groups of anions. On the basis of RDFs shown in Figure 3, the first coordination shell for Oc atoms is defined as a sphere with a 2.8 Å radius. We define a function H(t) which keeps track of atoms coordinating each Li+. If at time t the Oc atom j is within the coordinating shell of the ith Li+, then Hij(t) = 1.0; otherwise, Hij(t) = 0.0. Calculating the time dependent autocorrelation function Hij(t)Hij(0) allows us to determine the time scale over which the coordination of Li+ ions is changing; i.e., the Li+ changes the oxygen atoms coordinating it. Figure 11a shows the decay of this autocorrelation function for amorphous Li2BDC at several different temperatures. At high temperatures, Li+ changes its coordination shell on the time scale of a couple nanoseconds. As the temperature goes down and the system approaches the glass transition, this time scale increases with scaling similar to that of the Li+ self-diffusion coefficient shown in Figure 9. As can be seen from Figure 11a, the obtained autocorrelation functions can be fit reasonably well with the

̂ = ⟨εij̅ εkl̅ ⟩ Sijkl

V kBT

(4)

Under the appropriate choice of contracted stress/strain vectors,55 the rank 4 compliance tensor of eq 4 may be compactly represented in Voigt notation, yielding a rank 2 tensor which is the matrix inverse of the elastic stiffness tensor, allowing for calculation of standard elastic constants (Cij), as well as the Voigt56 and Reuss57 estimates58 of the overall bulk (K) and shear (G) moduli of the material. To calculate the elastic stiffness coefficients of the simulation cell, it is therefore necessary to adequately sample the fluctuations of strain vectors for a system that is in equilibrium or in stationary condition. These strain vectors, in turn, are derived from the spontaneous thermodynamic fluctuations of the simulation cell, as has previously been outlined for a variety of molecular crystalline systems.59−61 To allow accurate sampling of the underlying fluctuations of the unit cell, we rely on a combined molecular-dynamics/ Monte Carlo (MD/MC) simulation approach.59,61 In this

Figure 11. (a) Residence autocorrelation function for Oc atoms in the first coordination shell of Li+ in Li2BDC amorphous melts (symbols). Solid lines show the KWW fits. (b) Comparison of this autocorrelation function for Li2BDC and Li2EDC amorphous melts at 500 K. 16104

DOI: 10.1021/acs.jpcc.7b04247 J. Phys. Chem. C 2017, 121, 16098−16109

Article

The Journal of Physical Chemistry C

distribution of cell dimensions and then, using this combined distribution, calculate Cij, K, and G. In Figure 13, we show several Cii components as well as the bulk and shear moduli for the ordered Li2EDC at 393 K extracted from simulations using the two analysis approaches discussed above. In the figure, we show properties obtained from the analysis of fluctuations in each replica (symbols) as well as properties where fluctuations from all replicas were used to create a master distribution before calculating mechanical properties (lines). Figure 13a clearly shows that properties extracted from individual replicas show a noticeable scatter between different replicas, but more importantly, they are systematically higher than the values extracted using the master distribution. The latter can be explained by the relatively short MD/MC sampling of each replica. Starting from a specific initial configuration, the 3 ns MD + 20,000 MC steps simulation of a replica is not sufficient to allow the system to explore all possible simulation cell dimensions/shapes that would adequately sample the corresponding stress/strain coupling in the system. Instead, each replica samples a relatively narrow portion of that distribution. Therefore, analysis of individual replicas results in predicting a “stiffer” mechanical response of studied SEIs. In contrast, if we use all replica data for the analysis, we obtain a more realistic sampling of possible stress/ strain coupling. The initial configurations for replicas are distributed along a very long MD simulation trajectory and therefore allow the sampling over a much broader range of stress/strain correlations. In this way, the underlying long-time MD simulation provides the long-time component of the sample cell fluctuations, while the MD/MC simulations provide shorter time, higher frequency fluctuations around the longtime correlated sample configurations, allowing for efficient sampling of a far larger subset of expected cell fluctuations. The difference between two approaches is schematically illustrated in Figure 12. The results in Figure 13 further illustrate the necessity of long trajectories or multiple independent realizations for glassy systems to reliably extract mechanical properties from simulations. Elastic Constants and Moduli. Using the second described above approach, i.e., simultaneous analysis of cell fluctuations from all replicas, we have analyzed the mechanical properties of investigated SEIs in the 333−450 K temperature range. In this temperature range, both the amorphous and ordered phases are stable (at least on the time scale accessible to our simulations) and hence their mechanical properties can be compared with each other. Table 1 summarizes the obtained values for elastic constants as well as the bulk and shear moduli. Values reported in Table 1 have error bars around 5% which were estimated by monitoring the change in elastic constants and moduli with addition of each replica to comprise the master distribution. After inclusion of five replicas, the extracted properties were not changing more than 5% with inclusion of additional replicas to the statistics, indicating the convergence of sampled distributions of cell fluctuations. First, examination of the ordered Li2EDC system indicates that elastic constants describing the stress/strain coupling in the direction perpendicular to Li+ layers (C11) is noticeably larger than the constants corresponding to two other directions (C22 and C33), indicating that the direction perpendicular to Li+ layering has the largest stiffness. The elastic constants responsible for the shear deformations (C44−C66) are almost an order of magnitude lower, indicating that formed SEIs have a relatively low shear modulus. Indeed, Table 1 shows that,

approach, the system evolution is accomplished using standard MD simulations with fixed shape/dimensions of the simulation cell (i.e., in the NVT ensemble). Periodically (e.g., every 1.5 ps), this evolution is interrupted with several (in this work 10) MC attempts to randomly perturb both the simulation cell lattice vector lengths and angles between them relative to the current simulation cell, followed by center-of-mass projection of the molecule distribution within the cell based on the new cell geometry. This MD/MC scheme allows more efficient sampling of the unit cell dimension distributions compared to standard fully flexible cell MD simulations in the NPT ensemble, as has been demonstrated in our previous investigation of organic crystals. However, the modeling of dynamically sluggish or glassy systems with low symmetry or amorphous ordering, such as the SEIs investigated here, introduces an additional challengethe necessity of sampling unit cell fluctuations around a representative set of system configurations/realizations. To combat this issue, we utilized multiple (up to 8) initial configurations extracted from a long-length trajectory obtained using standard MD simulation and separated from each other by Δt0 time intervals. For example, from a 250 ns production run of the Li2EDC system at 393 K,33 we have extracted six system configurations separated from each other by Δt0 = 30 ns. Subsequently, each of those configurations was subjected to MD/MC simulation over ΔtMD−MC = 3 ns of MD steps with 10 MC attempts to change cell dimensions every 1.5 ps (i.e., a total of 20,000 MC trials for each initial configuration). For 333 and 450 K systems, we used the same initial configurations as those for 393 K with the addition of a 5 ns annealing run (using standard MD simulation in the NPT ensemble) at the desired temperature before conducting the MD/MC simulations. This simulation protocol is schematically illustrated in Figure 12.

Figure 12. Schematic illustration of the combined MD−MC sampling employed to predict mechanical properties of model SEI compounds.

When the MD/MC simulations of all replicas have been performed, we then have two possible approaches to analyze the cell dimension fluctuations and extract corresponding mechanical properties. These two approaches are also illustrated in Figure 12. First, one can analyze fluctuations (i.e., cell dimensions distributions) for each replica, extract corresponding Cij, K, and G, and then simply average the extracted properties over all replicas. In the second approach, one can use fluctuations from all replicas to create a master 16105

DOI: 10.1021/acs.jpcc.7b04247 J. Phys. Chem. C 2017, 121, 16098−16109

Article

The Journal of Physical Chemistry C

Figure 13. Elastic constants (a) and bulk and shear moduli (b) obtained from the analysis of MD/MC simulations of the ordered Li2EDC at 393 K.

significantly higher than the plateau modulus (∼105−107 Pa) and comparable to values of the glassy modulus (couple GPa) for typical bulk polymers.62 The shear modulus of the Li2O SEI component and solid electrolytes (LLZO garnet, LLTIO perovskite) is an order of magnitude higher than the shear modulus of Li2BDC and Li2EDC at 393 K.63 When these results are considered in the context of the ability of SEI and electrolyte in resisting the growth of Li metal dendrites at the electrode surface during battery cycling, they indicate that Li2BDC and Li2EDC SEI components are not sufficiently effective in preventing dendrite growth. According to the Newman and Monroe model, an SEI with a shear modulus larger than 7.0 GPa is necessary to push dendrites back and prevent dendrite growth.64 Clearly, SEIs formed out of dilithium alkyldicarbonates have a lower shear modulus. Considering that SEIs can also contain entrapped solvent molecules or longer oligomers, it is unlikely that the outer SEI region, comprised mostly of dilithium alkyldicarbonates, provides sufficient mechanical strength to prevent dendrite growth based only on consideration of mechanical properties. This is consistent with poor lithium cycling in carbonate solvents.65 In contrast, recent experimental works showed that polymeric systems with low shear modulus are capable of significant reduction of dendrite growth in Li-ion batteries.66 A more systematic understanding of coupling between SEI structure and dendrite growth mechanisms is therefore needed.

Table 1. Elastic Constants and Bulk and Shear Moduli (in GPa) for the Investigated Model SEI Compounds from MD/ MC Simulations Using All Constructed Replicas Simultaneously Li2EDC-ordered C11 C22 C33 C44 C55 C66 KVoigt KReuss GVoigt GReuss

Li2BDC-ordered

Li2BDC-amorphous

333 K

393 K

450 K

393 K

393 K

21.6 18.0 17.9 4.5 5.4 4.5 13.2 12.5 4.7 4.2

19.0 13.7 13.3 2.8 2.3 2.6 10.6 10.0 3.0 2.6

15.5 11.3 11.2 1.9 1.8 2.2 8.3 8.0 2.5 2.0

10.3 7.6 8.0 3.4 1.5 0.8 7.0 6.1 1.7 1.0

10.1 9.2 10.8 2.9 3.1 2.9 7.1 6.7 2.7 1.5

while the bulk modulus of the ordered Li2EDC varies between 8 and 13 GPa in the 333−450 K temperature range, the shear modulus ranges between 2.0 and 4.7 GPa in the same temperatures range. A comparison of Li2EDC and Li2BDC ordered phases at 393 K indicates that the latter is a softer material with a bulk modulus around 6.5 GPa (compared to 10.3 GPa for Li2EDC) and a shear modulus around 1.4 GPa (compared to 2.8 GPa for Li2EDC). Note that in the Li2BDC there is a much larger difference in the shear components of elastic constants, which is likely related to a much more ordered structure of this SEI compared to Li2EDC (compare the snapshots for these systems in Figure 6). The well-defined layering of Li+ in the Li2BDC results in one shear direction (parallel to the layers) that has a very small shear modulus (∼0.5 GPa). A comparison of ordered and amorphous Li2BDC phases at 393 K indicates that, while the bulk modulus and most of the elastic constants are similar, the shear modulus (and the corresponding elastic constants) are noticeably larger for the amorphous phase. The average shear modulus for amorphous Li2BDC is 2.1 GPa (compared to 1.4 GPa for the ordered Li2BDC). This is consistent with the fact that there is no layering (i.e., “easy” shear planes) in the amorphous phase where Li+ form a bicontinuous phase, as shown in Figure 6. The predicted mechanical properties for Li2BDC and Li2EDC presented above, to the best of our knowledge, are the first direct estimate of mechanical characteristics of dilithium alkyldicarbonate SEIs as a function of temperature. Our results indicate that the bulk modulus of these materials is

4. CONCLUSIONS An extensive investigation of structural, dynamic, and mechanical properties of Li2EDC and Li2BDC model SEIs was conducted as a function of temperature utilizing atomistic MD simulations. At temperatures below 450 K, these materials are glassy and the mobility of Li+ cations becomes decoupled from anions. At lower temperatures, the anion comprised matrix is frozen, while Li+ move through the matrix following an Arrhenius-like temperature dependence with an activation energy of about 0.6 eV for both amorphous and ordered materials. This activation energy is noticeably larger than the activation energy for Li+ transition from carbonate electrolyte to SEI which has been estimated to be around 0.42−0.46 eV in our previous work using the same models for SEI components and EC:DMC/LiPF6 electrolyte. However, this value is in good agreement with experimentally estimated activation energy for the overall interfacial resistance at the anode side. The lithium transport in Li2EDC and Li2BDC has an activation energy slightly lower than the other SEI compound Li2CO3 but 16106

DOI: 10.1021/acs.jpcc.7b04247 J. Phys. Chem. C 2017, 121, 16098−16109

Article

The Journal of Physical Chemistry C

computing resources. The authors also would like to thank Dr. Kang Xu (ARL) for productive discussions.

significantly higher than in the state-of-the-art ceramic electrolytes that showed a much lower activation energy of 0.2−0.3 eV.67 At lower temperatures, the Li2EDC and Li2BDC SEIs exist either in amorphous or ordered (smectic-like) phases. Slow SEI formation is expected to allow sufficient time for Li2EDC to form ordered structures, in agreement with experimentally observed orthorhombic crystals, as discussed in ref 46. Fast SEI formation, however, would lead to more disordered and glasslike SEI. The higher conductivity of the ordered phase compared to the amorphous phase predicted by our MD simulations is consistent with the observed better performance of batteries when the SEI formation cycle is performed at a very slow rate. Increasing the length of the alkyl spacer from two (EDC) to four (BDC) carbons between carbonate groups on the anion results in the formation of more defined layers of Li+, leading to an increased mobility of Li+ inside the layers. However, at the same time, it leads to the formation of a shearing plane, resulting in relatively low shear strength along the plane of the formed layers. The bulk modulus is on the order of 7−13 GPa, while the shear modulus is between 1.5 and 4.7 GPa. This shear modulus is significantly higher than the shear modulus of the state-of-the-art block copolymer rubbery electrolytes68 but lower than the Li2O SEI component and many of the solid electrolytes. Comparison of SEIs comprised of BDC or EDC anions indicated that longer alkyl spacers result in lower SEI modulus, i.e., making the materials softer and hence reducing their ability to prevent dendrite growth during lithium plating. Therefore, the presence of the Li2BDC components in SEI should be less desirable for battery operation.





(1) Gauthier, M.; Carney, T. J.; Grimaud, A.; Giordano, L.; Pour, N.; Chang, H.-H.; Fenning, D. P.; Lux, S. F.; Paschos, O.; Bauer, C.; Maglia, F.; Lupart, S.; Lamp, P.; Shao-Horn, Y. Electrode−Electrolyte Interface in Li-Ion Batteries: Current Understanding and New Insights. J. Phys. Chem. Lett. 2015, 6, 4653−4672. (2) Xu, K. Nonaqueous Liquid Electrolytes for Lithium-Based Rechargeable Batteries. Chem. Rev. 2004, 104, 4303−4417. (3) Zhang, S. S.; Xu, K.; Jow, T. R. Charge and Discharge Characteristics of a Commercial LiCoO2-Based 18650 Li-Ion Battery. J. Power Sources 2006, 160, 1403−1409. (4) Jow, T. R.; Marx, M. B.; Allen, J. L. Distinguishing Li+ Charge Transfer Kinetics at Nca/Electrolyte and Graphite/Electrolyte Interfaces, and NCA/Electrolyte and LFP/Electrolyte Interfaces in Li-Ion Cells. J. Electrochem. Soc. 2012, 159, A604−A612. (5) Ogumi, Z. Interfacial Reactions of Lithium-Ion Batteries. Electrochemistry 2010, 78, 319−324. (6) Xu, K.; von Wald Cresce, A. Li+ Solvation/Dissolvation Dictates Intephasial Processes on Graphitic Anode in Li Ion Cells. J. Mater. Res. 2012, 27, 2327−2341. (7) Zorba, V.; Syzdek, J.; Mao, X. L.; Russo, R. E.; Kostecki, R. Ultrafast Laser Induced Breakdown Spectroscopy of Electrode/ Electrolyte Interfaces. Appl. Phys. Lett. 2012, 100, 234101. (8) Aurbach, D.; Ein-Ely, Y.; Zaban, A. The Surface Chemistry of Lithium Electrodes in Alkyl Carbonate Solutions. J. Electrochem. Soc. 1994, 141, L1−3. (9) Aurbach, D.; Zaban, A.; Schecheter, A.; Ein-Eli, Y.; Zinigrad, E.; Markovsky, B. The Study of Electrolyte Solutions Based on Ethylene and Diethyl Carbonates for Rechargeable Li Batteries. I. Li Metal Anodes. J. Electrochem. Soc. 1995, 142, 2873−2882. (10) Borodin, O.; Smith, G. D. Quantum Chemistry and Molecular Dynamics Simulation Study of Dimethyl Carbonate: Ethylene Carbonate Electrolytes Doped with LiPF6. J. Phys. Chem. B 2009, 113, 1763−1776. (11) Seo, D. M.; Reininger, S.; Kutcher, M.; Redmond, K.; Euler, W. B.; Lucht, B. L. Role of Mixed Solvation and Ion Pairing in the Solution Structure of Lithium Ion Battery Electrolytes. J. Phys. Chem. C 2015, 119, 14038−14046. (12) Bogle, X.; Vazquez, R.; Greenbaum, S.; Cresce, A. v. W.; Xu, K. Understanding Li+−Solvent Interaction in Nonaqueous Carbonate Electrolytes with 17O NMR. J. Phys. Chem. Lett. 2013, 4, 1664−1668. (13) von Cresce, A.; Xu, K. Preferential Solvation of Li+ Directs Formation of Interphase on Graphitic Anode. Electrochem. Solid-State Lett. 2011, 14, A154−A156. (14) Morita, M.; Asai, Y.; Yoshimoto, N.; Ishikawa, M. A Raman Spectroscopic Study of Organic Electrolyte Solutions Based on Binary Solvent Systems of Ethylene Carbonate with Low Viscosity Solvents Which Dissolve Different Lithium Salts. J. Chem. Soc., Faraday Trans. 1998, 94, 3451−3456. (15) Borodin, O.; Olguin, M.; Ganesh, P.; Kent, P. R. C.; Allen, J. L.; Henderson, W. A. Competitive Lithium Solvation of Linear and Cyclic Carbonates From Quantum Chemistry. Phys. Chem. Chem. Phys. 2016, 18, 164−175. (16) Vatamanu, J.; Borodin, O.; Smith, G. D. Molecular Dynamics Simulation Studies of the Structure of a Mixed Carbonate/LiPF6 Electrolyte Near Graphite Surface as a Function of Electrode Potential. J. Phys. Chem. C 2012, 116, 1114−1121. (17) Borodin, O.; Olguin, M.; Spear, C. E.; Leiter, K.; Knap, J. Towards High Throughput Screening of Electrochemical Stability of Battery Electrolytes. Nanotechnology 2015, 26, 354003. (18) Dudney, N. J.; Hagaman, E. W.; Veith, G. M.; Gill, L. W.; Sacci, R. L. Lithiated and Passivated Lithium Ion Battery Anodes. Patent US20160218351 A1, 2016. (19) Aurbach, D.; Weismann, I.; Schechter, A. X-ray Photoelectron Spectroscopy Studies of Lithium Surfaces Prepared in Several Important Electrolyte Solutions. A Comparison with Previous Studies

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcc.7b04247. The length of simulation trajectories and the mean square displacements of anion (PDF)



REFERENCES

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. ORCID

Dmitry Bedrov: 0000-0002-3884-3308 Oleg Borodin: 0000-0002-9428-5291 Notes

The authors declare no competing financial interest.



ACKNOWLEDGMENTS The authors gratefully acknowledge support for the project sponsored by the Army Research Laboratory under Cooperative Agreement Number W911NF-12-2-0023. The views and conclusions contained in this document are those of the authors and should not be interpreted as representing the official policies, either expressed or implied, of ARL or the U.S. Government. The U.S. Government is authorized to reproduce and distribute reprints for Government purposes notwithstanding any copyright notation herein. D.B. and J.B.H. would like to acknowledge the Center of High Performance Computing at the University of Utah for technical support and allocation of 16107

DOI: 10.1021/acs.jpcc.7b04247 J. Phys. Chem. C 2017, 121, 16098−16109

Article

The Journal of Physical Chemistry C by Fourier Transform Infrared Spectroscopy. Langmuir 1996, 12, 3991−4007. (20) Aurbach, D.; Gofer, Y.; Ben-Zion, M.; Aped, P. The Behavior of Lithium Electrodes in Propylene and Ethylene Carbonate: Te Major Factors That Influence Li Cycling Efficiency. J. Electroanal. Chem. 1992, 339, 451−471. (21) Xu, K.; Lam, Y.; Zhang, S. S.; Jow, T. R.; Curtis, T. B. Solvation Sheath of Li+ in Nonaqueous Electrolytes and Its Implication of Graphite/Electrolyte Interface Chemistry. J. Phys. Chem. C 2007, 111, 7411−7421. (22) Nie, M. Y.; Chalasani, D.; Abraham, D. P.; Chen, Y. J.; Bose, A.; Lucht, B. L. Lithium Ion Battery Graphite Solid Electrolyte Interphase Revealed by Microscopy and Spectroscopy. J. Phys. Chem. C 2013, 117 (3), 1257−1267. (23) Naji, A.; Ghanbaja, J.; Humbert, B.; Willmann, P.; Billaud, D. Electroreduction of Graphite in LiClO4-Ethylene Carbonate Electrolyte. Characterization of the Passivating Layer by Transmission Electron Microscopy and Fourier-Transform Infrared Spectroscopy. J. Power Sources 1996, 63, 33−39. (24) Michan, A. L.; Leskes, M.; Grey, C. P. Voltage Dependent Solid Electrolyte Interphase Formation in Silicon Electrodes: Monitoring the Formation of Organic Decomposition Products. Chem. Mater. 2016, 28, 385−398. (25) Leung, K.; Soto, F.; Hankins, K.; Balbuena, P. B.; Harrison, K. L. Stability of Solid Electrolyte Interphase Components on Lithium Metal and Reactive Anode Material Surfaces. J. Phys. Chem. C 2016, 120, 6302−6313. (26) Wang, Y.; Nakamura, S.; Ue, M.; Balbuena, P. B. Theoretical Studies To Understand Surface Chemistry on Carbon Anodes for Lithium-Ion Batteries: Reduction Mechanisms of Ethylene Carbonate. J. Am. Chem. Soc. 2001, 123, 11708. (27) Leung, K.; Budzien, J. L. Ab Initio Molecular Dynamics Simulations of the Initial Stages of Solid-Electrolyte Interphase Formation on Lithium Ion Battery Graphitic Anodes. Phys. Chem. Chem. Phys. 2010, 12, 6583−6586. (28) Leung, K.; Rempe, S. B.; Foster, M. E.; Ma, Y.; Martinez del la Hoz, J. M.; Sai, N.; Balbuena, P. B. Modeling Electrochemical Decomposition of Fluoroethylene Carbonate on Silicon Anode Surfaces in Lithium Ion Batteries. J. Electrochem. Soc. 2014, 161, A213−A221. (29) Leung, K. Electronic Structure Modeling of Electrochemical Reactions at Electrode/Electrolyte Interfaces in Lithium Ion Batteries. J. Phys. Chem. C 2013, 117, 1539−1547. (30) Bedrov, D.; Smith, G. D.; van Duin, A. C. T. Reactions of Singly-Reduced Ethylene Carbonate in Lithium Battery Electrolytes. A Molecular Dynamics Study Using the ReaxFF. J. Phys. Chem. A 2012, 116, 2978−2985. (31) Islam, M. M.; van Duin, A. C. T. Reductive Decomposition Reactions of Ethylene Carbonate by Explicit Electron Transfer from Lithium: An e-ReaxFF Molecular Dynamics Study. J. Phys. Chem. C 2016, 120, 27128−27134. (32) Borodin, O.; Smith, G. D.; Fan, P. Molecular Dynamics Simulations of Lithium Alkyl Carbonates. J. Phys. Chem. B 2006, 110, 22773. (33) Borodin, O.; Zhuang, G. V.; Ross, P. N.; Xu, K. Molecular Dynamics Simulations and Experimental Study of Lithium Ion Transport in Dilithium Ethylene Dicarbonate. J. Phys. Chem. C 2013, 117, 7433−7444. (34) Borodin, O.; Bedrov, D. Interfacial Structure and Dynamics of the Lithium Alkyl Dicarbonate SEI Components in Contact with the Lithium Battery Electrolyte. J. Phys. Chem. C 2014, 118, 18362−18371. (35) Jorn, R.; Kumar, R.; Abraham, D. P.; Voth, G. A. Atomistic Modeling of the Electrode−Electrolyte Interface in Li-Ion Energy Storage Systems: Electrolyte Structuring. J. Phys. Chem. C 2013, 117, 3747−3761. (36) Takenaka, N.; Suzuki, Y.; Sakai, H.; Nagaoka, M. On Electrolyte-Dependent Formation of Solid Electrolyte Interphase Film in Lithium-Ion Batteries: Strong Sensitivity to Small Structural

Difference of Electrolyte Molecules. J. Phys. Chem. C 2014, 118, 10874−10882. (37) Single, F.; Horstmann, B.; Latz, A. Dynamics and Morphology of Solid Electrolyte Interphase (SEI). Phys. Chem. Chem. Phys. 2016, 18, 17810−17814. (38) Andreev, O. L.; Raskovalov, A. A.; Larin, A. V. A Molecular Dynamics Simulation of Lithium Fluoride: Volume Phase and Nanosized Particle. Russ. J. Phys. Chem. A 2010, 84, 48−52. (39) Chen, Y. C.; Ouyang, C. Y.; Song, L. J.; Sun, Z. L. Electrical and Lithium Ion Dynamics in Three Main Components of Solid Electrolyte Interphase from Density Functional Theory Study. J. Phys. Chem. C 2011, 115, 7044−7049. (40) Iddir, H.; Curtiss, L. A. Li Ion Diffusion Mechanisms in Bulk Monoclinic Li2CO3 Crystals from Density Functional Studies. J. Phys. Chem. C 2010, 114, 20903−20906. (41) Shi, S. Q.; Lu, P.; Liu, Z. Y.; Qi, Y.; Hector, L. G.; Li, H.; Harris, S. J. Direct Calculation of Li-Ion Transport in the Solid Electrolyte Interphase. J. Am. Chem. Soc. 2012, 134, 15476−15487. (42) Shang, S. L.; Hector, L. G.; Shi, S. Q.; Qi, Y.; Wang, Y.; Liu, Z. K. Lattice Dynamics, Thermodynamics and Elastic Properties of Monoclinic Li2CO3 from Density Functional Theory. Acta Mater. 2012, 60, 5204−5216. (43) Pan, J.; Cheng, Y. T.; Qi, Y. General Method to Predict VoltageDependent Ionic Conduction in a Solid Electrolyte Coating on Electrodes. Phys. Rev. B: Condens. Matter Mater. Phys. 2015, 91, 134116. (44) Shi, S. Q.; Qi, Y.; Li, H.; Hector, L. G. Defect Thermodynamics and Diffusion Mechanisms in Li2CO3 and Implications for the Solid Electrolyte Interphase in Li-Ion Batteries. J. Phys. Chem. C 2013, 117, 8579−8593. (45) Li, Y.; Leung, K.; Qi, Y. Computational Exploration of the LiElectrode|Electrolyte Interface in the Presence of a Nanometer Thick Solid-Electrolyte Interphase Layer. Acc. Chem. Res. 2016, 49, 2363− 2370. (46) Monroe, C.; Newman, J. The Impact of Elastic Deformation on Deposition Kinetics at Lithium/Polymer Interfaces. J. Electrochem. Soc. 2005, 152, A396−A404. (47) Xu, K.; Zhuang, G. V.; Allen, J. L.; Lee, U.; Zhang, S. S.; Ross, P. N.; Jow, T. R. Syntheses and Characterization of Lithium Alkyl Monoand Dicarbonates as Components of Surface Films in Li-Ion Batteries. J. Phys. Chem. B 2006, 110, 7708−7719. (48) Martyna, G. J.; Tuckerman, M. E.; Tobias, D. J.; Klein, M. L. Explicit Reversible Integrators for Extended Systems Dynamics. Mol. Phys. 1996, 87, 1117−1157. (49) Tasaki, K.; Goldberg, A.; Lian, J. J.; Walker, M.; Timmons, A.; Harris, S. J. Solubility of Lithium Salts Formed on the Lithium-Ion Battery Negative Electrode Surface in Organic Solvents. J. Electrochem. Soc. 2009, 156, A1019−A1027. (50) Islam, M. M.; Ostadhossein, A.; Borodin, O.; Yeates, A. T.; Tipton, W. W.; Hennig, R. G.; Kumar, N.; van Duin, A. C. T. ReaxFF Molecular Dynamics Simulations on Lithiated Sulfur Cathode Materials. Phys. Chem. Chem. Phys. 2015, 17, 3383−3393. (51) Borodin, O. Molecular Modeling of Electrolytes. In Electrolytes for Lithium and Lithium-Ion Batteries; Jow, T. R., Xu, K., Borodin, O., Ue, M., Eds.; Springer: New York, 2014; Vol. 58, pp 371−401. (52) Li, C. L.; Gu, L.; Maier, J. Enhancement of the Li Conductivity in LiF by Introducing Glass/Crystal Interfaces. Adv. Funct. Mater. 2012, 22, 1145−1149. (53) Williams, G.; Watts, D. C. Non-symmetrical Dielectric Relaxation Behaviour Arising from a Simple Empirical Decay Function. Trans. Faraday Soc. 1970, 66, 80−85. (54) Parrinello, M.; Rahman, A. Strain Fluctuations and Elastic Constants. J. Chem. Phys. 1982, 76, 2662−2669. (55) Tsai, S. W. AFML Technical Report No. AFMLTR-66-149, Nov 1966. (56) Voigt, W. Theoretical Studies in the Elastic Behavior of Crystals I. Abh. d. Kon. Ges. d. Wiss. Gottingen 1892, 34, 3−52. 16108

DOI: 10.1021/acs.jpcc.7b04247 J. Phys. Chem. C 2017, 121, 16098−16109

Article

The Journal of Physical Chemistry C (57) Reuss, A. Berechnung der Fließgrenze von Mischkristallen auf Grund der Plastizitätsbedingung für Einkristalle. Z. Angew. Math. Mech. 1929, 9, 49−58. (58) Hill, R. The Elastic Behaviour of a Crystalline Aggregate. Proc. Phys. Soc., London, Sect. A 1952, 65, 349−355. (59) Sewell, T. D.; Menikoff, R.; Bedrov, D.; Smith, G. D. Molecular Dynamics Simulation Study of Elastic Properties of HMX. J. Chem. Phys. 2003, 119, 7417−7426. (60) Bedrov, D.; Borodin, O.; Smith, G. D.; Sewell, T. D.; Dattelbaum, D. M.; Stevens, L. L. A Molecular Dynamics Simulation Study of Crystalline 1,3,5-triamino-2,4,6-trinitrobenzene as a Function of Pressure and Temperature. J. Chem. Phys. 2009, 131, 224703. (61) Hooper, J. B.; Borodin, O.; Schneider, S. Insight into Hydrazinium Nitrates, Azides, Dicyanamide, and 5-Azidotetrazolate Ionic Materials from Simulations and Experiments. J. Phys. Chem. B 2011, 115, 13578−13592. (62) Shaw, M. T.; MacKnight, W. J. Introduction to Polymer Viscoelasticity, 3rd ed.; John Wiley & Sons: Hoboken, NJ, 2005. (63) Deng, Z.; Wang, Z.; Chu, I.-H.; Luo, J.; Ong, S. P. Elastic Properties of Alkali Superionic Conductor Electrolytes from First Principles Calculations. J. Electrochem. Soc. 2016, 163, A67−A74. (64) Monroe, C.; Newman, J. The Impact of Elastic Deformation on Deposition Kinetics at Lithium/Polymer Interfaces. J. Electrochem. Soc. 2005, 152, A396−A404. (65) Qian, J.; Henderson, W. A.; Xu, W.; Bhattacharya, P.; Engelhard, M.; Borodin, O.; Zhang, J.-G. High Rate and Stable Cycling of Lithium Metal Anode. Nat. Commun. 2015, 6, 6362. (66) Khurana, R.; Schaefer, J. L.; Archer, L. A.; Coates, G. W. Suppression of Lithium Dendrite Growth Using Cross-Linked Polyethylene/Poly(ethylene oxide) Electrolytes: A New Approach for Practical Lithium-Metal Polymer Batteries. J. Am. Chem. Soc. 2014, 136, 7395−7402. (67) Wang, Y.; Richards, W. D.; Ong, S. P.; Miara, L. J.; Kim, J. C.; Mo, Y.; Ceder, G. Design Principles for Solid-State Lithium Superionic Conductors. Nat. Mater. 2015, 14, 1026−1031. (68) Schauser, N. S.; Harry, K. J.; Parkinson, D. Y.; Watanabe, H.; Balsara, N. P. Lithium Dendrite Growth in Glassy and Rubbery Nanostructured Block Copolymer Electrolytes. J. Electrochem. Soc. 2015, 162, A398−A405.

16109

DOI: 10.1021/acs.jpcc.7b04247 J. Phys. Chem. C 2017, 121, 16098−16109