Molecular Dynamics of Lithium Ion Transport in a Model Solid Electrolyte Interphase
Ajay Muralidharan Mangesh I. Chaudhari Lawrence R. Pratt Susan B. Rempe
physics.chem-ph

Automatically Extracting Web API Specifications from HTML Documentation
    Jinqiu Yang1,

Erik Wittern2,
Annie T.T. Ying32,
Julian Dolby2,
Lin Tan1
1 University of Waterloo, Waterloo, Canada
2 IBM Research, Yorktown Heights, NY, USA
3 EquitySim, Vancouver, Canada
j223yang@uwaterloo.ca, witternj@us.ibm.com, annie.ying@gmail.com, dolby@us.ibm.com, lintan@uwaterloo.ca

    May 3, 2018
=========================================================================================================================================================================================================================================================================================================



empty
 


§ INTRODUCTION

During charging and discharging cycles of lithium ion batteries, a solid
electrolyte interphase (SEI) layer forms  on the negative electrode due
to decomposition of solvents 

like ethylene carbonate (EC). The SEI
layer is a complex organic material and its composition is not
operationally set.<cit.> Nevertheless, dilithium ethylene
dicarbonate (Li_2EDC) has been identified experimentally as the
primary component of the outer part of the SEI
layer.<cit.> 


Ab initio molecular dynamics (AIMD) simulations, electronic structure
calculations<cit.> and reactive force
field simulations<cit.> on the decomposition of ethylene
carbonate (EC) on anode surfaces concur with those experimental results.
Experimental observations also show that the SEI layer protects the
electrode from further decomposition by blocking electron transport
while simultaneously allowing transport of Li^+ ions between the electrode and
electrolyte solution. A better understanding of the  transport mechanism
 of Li^+ in EDC may lead to modified SEI layers with improved lithium
ion battery performance. 



Molecular dynamics (MD) studies on model SEI layers carried out over
long time-scales may shed new light into the mechanism of transport of
Li^+ ions within the SEI layer.  Borodin, et
al.,<cit.> performed MD
calculations using a specialized polarizable force field to obtain 
transport properties of Li^+ ion in a model SEI layer composed of
ordered and amorphous Li_2EDC. Since polarizable force fields are not
readily available in standard molecular dynamics packages, we have
instead identified non-polarizable force field
parameters<cit.> for simulation of the
Li_2EDC model of the SEI layer. The microsecond time-scales studied
here, longer than earlier work,<cit.> provide additional
insight into structural and transport properties of Li^+ ions in this
model SEI.





§ RESULTS AND DISCUSSION

We summarize the structural and transport properties from MD studies of a model SEI layer with 256 Li_2EDC moieties (Figure <ref>). 
In addition, simulations of a dilute Li^+ ion in EC (single Li^+ ion solvated by 249 EC) provide a perspective for comparison. 


 §.§ Structural data

The radial distribution functions (rdfs) and running coordination
numbers  (Figure <ref>), involving Li^+, carbonyl oxygen
(Oc) and ether oxygen (Oe) of EDC and EC are compared at several temperatures. For EDC, the rdfs become less structured with
increasing temperature, but the peak positions and overall coordination
numbers of the first peak change only slightly. This structural
robustness suggests an amorphous glassy matrix for the SEI. These
results compare well with recent polarizable force field
results,<cit.>  supporting the
applicability of the present non-polarizable force field for these
structural characteristics. In EDC, the peak position for the Li-Oe
distribution is shifted to a slightly larger value than Li-Oc because
the charge on Oe is 40% smaller than that of Oc. 

In the case of EC, the peak height increases with temperature and the
peak position of the Li-Oe is farther out due to strong interactions
between Li^+ and Oc. The structure around Li^+ changes significantly
with temperature, as highlighted by the running coordination number. In contrast to glassy EDC, almost one additional carbonyl oxygen (Oc) atom interacts with Li^+
at lower temperature for EC solvent.   



 §.§ Mean-squared displacements of Li^+

The mean-squared displacements (MSD, Figure <ref>)
of Li^+ ion in EDC and liquid EC distinguish three distinct dynamical regimes:
ballistic at short times, trapping at intermediate times, and diffusive
at long times. Trapping of Li^+in EDC is qualitatively different than in 
 liquid EC, and the trapping times in EDC diminish as temperature increases and the 
glassy EDC matrix softens.
Diffusivities of  Li^+ are extracted from the slope of the diffusive regime of MSD and then extrapolated to
low temperatures using an Arrhenius fit (see Supplementary Information). The diffusivity of Li^+ in EDC at
333 K is 10^-12 cm^2/s, which is in agreement with Borodin et. al.<cit.> The conductivity of Li_2EDC obtained from the Nernst-Einstein equation (Supplementary Information) is 4.5 x 10^-9 S/cm, which is also in agreement with experiment.<cit.>



 



 §.§ Time correlation functions for Li^+ transport

The van Hove time  correlation function

    G(r,t)=1/N⟨∑_i=1^N∑_j=1^Nδ(r+r_j(0)-r_i(t))⟩
 
for Li^+ ion describes the probability of finding, at time t, a
Li^+ ion displaced by r from its initial position.

This
van Hove function can be split into self and distinct parts, G(r,t)=G_s(r,t)+
G_d(r,t), with
the latter taking the form

    G_d(r,t)=1/N⟨∑_i≠j^N(r+r_j(0)-r_i(t))⟩ .

At t=0, the van Hove function reduces to the static pair correlation function,

    G(r,0)=δ(r) + ρ g(r) .



The self part of the van Hove function provides a jump
probability, and the natural initial approximation is the Gaussian model,

    G_s(r,t)= 3/2π⟨Δ r(t)^2⟩^3/2×exp - (3r^2/2⟨Δr(t)^2⟩) .

For fluids like EC, this Gaussian behavior should be reliable.  In
contrast, G_s(r,t) in the trapping regime of EDC indicates
(Figure <ref>) depletion of probability near the trap
boundaries r> ⟨Δ r(t)^2⟩^1/2, and replacement of
that probability at shorter and longer distances.  Deviation from the
Gaussian behavior can be characterized by the non-Gaussian
parameter<cit.>

    α(t)= 3⟨Δr(t)^4⟩/5⟨Δr(t)^2⟩^2-1 .


The van Hove self-correlation function is accurately Gaussian for liquid
EC, hence α(t)=0. In contrast, α(t) has non-zero values for glassy EDC
(Figure <ref>). For EDC, α(t) has a maximum  that decreases
with increasing temperature.  The mean-squared displacements of the
carbonyl carbons of EDC and EC molecules (Figure <ref>) further
verify the sluggish diffusion of EDC relative to EC. In the case of EDC, the
mean-squared displacement curves are relatively flat for all
temperatures, indicating little diffusion of the EDC matrix.   



Vineyard's convolution approximation,<cit.>

    G_d(r,t)≈∫^3 r ^' g(r ^') G_s(|r⃗-r⃗ ^'|,t),

provides an initial characterization for the distinct part of the
Li^+- Li^+ van Hove function (Figure <ref>). Here,
the convolution of the radial distribution function is made with the
self part of the van Hove function that is generated by dynamics of
Li^+ in EDC. This approximation is consistent with the idea that
G_d(r,t) is a dynamical counterpart to g(r), the radial
distribution function. The non-zero population in the core region
surrounding  r ≈ 0 describes correlation of Li^+ jumps;

that is, refilling a hole left by a Li^+ ion with a neighboring
Li^+ ion.

 
















 §.§ Vibrational power spectra of Li^+

The vibrational power spectra are obtained by spectral decomposition of the velocity autocorrelation (Figure <ref>, left). 
Since we are interested in Li^+ transport, this analysis is carried out  for Li^+ atoms exclusively. 
This fact distinguishes our results from the FTIR spectrum of EDC molecule that was reported previously.<cit.> 
The spectral distribution for Li^+ in EDC (Figure <ref>, right) is bi-modal, 
with a temperature dependence at the higher frequency mode.
Electronic structure calculations, using structures sampled from MD trajectories, confirm these modes and provide a molecular assignment (Figure <ref>).
The lower frequency mode (near 400 cm^-1) corresponds to motion of a
Li^+ ion trapped in a cage formed by its nearest neighbors. The higher frequency mode (near 700
cm^-1) corresponds to Li^+ ion picking up the scissoring motion of
a neighboring carbonate group. 


The Einstein frequency (ν _e) is obtained as the coefficient in the quadratic approximation to the velocity autocorrelation function at short times, 

    ⟨v⃗(t)·v⃗(0)⟩/⟨ v^2⟩≈ 1 - (ν _e  
    t)^2
 In a simple Einstein model, all atoms vibrate with a single frequency. Fittingly, for the bi-modal spectra, ν _e lies in-between the two significant 
modes.

  

 





§ CONCLUSIONS

Non-polarizable, classical force field parameters were used to study
 transport characteristics  of Li^+ in a model SEI layer
composed of EDC.  The structural results in EDC are consistent with
prior studies that use polarizable force fields. An advantage of
non-polarizable force fields is their ready availability in standard
simulation packages and accessibility to  MD studies over microsecond
timescales.  Thus, the dynamical characteristics presented here lay a
basis for careful molecular-scale examination of the mechanism of
transport of Li^+ ions in the SEI.

These observations  over microsecond simulation times provide new physical
insights.

Specifically, the results compare the glassy behavior of the
ethylene dicarbonate SEI matrix with the fluid behavior of liquid
ethylene carbonate (Figure <ref>). Further, the Li^+ MSDs
examined in the nanosecond time intervals distinguish Li^+ ion
trapping in cages formed by the EDC matrix. Our results establish the
sizes of the cages and the trapping lifetimes
(Figure <ref>), and also the dynamical motions of
the Li^+ ions when trapped (Figure <ref>).  The
vibrational frequency of a trapped ion (about 440 cm^-1) is
confirmed by electronic structure calculations
(Figure <ref>). Our results invalidate a naive Einstein
model of trapped ions that would be plausible otherwise.  The van Hove
correlation functions (Figures <ref> and
<ref>) provide information for analysis of the
correlation of Li^+ jumps.



§ METHODS

Li_2EDC (Figure <ref>) is known to be a dominant component of
the SEI layer in lithium ion batteries involving carbonate solvents.
Although Li_2EDC is synthesized in crystalline form, its structure at
the SEI layer is unknown.<cit.> We
constructed a system of 256 Li_2EDC moieties for our initial SEI
studies. This system size is identical to previous molecular simulations
performed using polarizable force fields.<cit.> For
comparison, we also simulated a single Li^+ ion solvated by 249 EC
(Figure <ref>). The GROMACS molecular dynamics simulation
package<cit.> was used for all simulations, and the necessary
topology files for EDC and EC were created using non-polarizable
all-atom optimized potentials for liquid simulations (OPLS-AA) force
field parameters.<cit.> The partial charges on EC atoms were
adjusted down to 80% to match experimental transport properties for
EC.<cit.>

The EDC  and Li^+ ions were randomly placed into the  simulation cell
and MD simulations were performed at 700, 500 and 333 K. Since EDC ions
are sluggish, configurations from the highest temperature calculation
were used to obtain starting points for further simulations;

cooled down to 500 K and subsequently to 333 K to study moderate
temperature phenomena. Thus, the results presented here are based on
amorphous configurations of the EDC/SEI layer. Although it is unclear
that Li_2EDC is crystalline at SEI layer, we have simulated ordered
layers and found that the solvent density and radial distribution
functions are not substantially changed compared with amorphous Li_2EDC. The ordered Li_2EDC is more
conductive compared to amorphous Li_2EDC,<cit.> but
formation of an ordered SEI structure is unlikely. Therefore, we here
discuss only results of amorphous Li_2EDC. 

Periodic boundary conditions  mimic the bulk environment in these
calculations. A Nose-Hoover thermostat<cit.> and a
Parrinello-Rahman<cit.> barostat were utilized to achieve
equilibration in the NpT ensemble at 1 atm pressure.  A 200 ns
production run at 700 K was carried out after initial energy
minimization and equilibration steps, then a 250 ns calculation at
500 K. Finally, a 1 μs trajectory at 333 K temperature was
constructed. Configurations were saved after each 1 ps of those
production runs. A separate 1 ns simulation with a sampling rate of 1 fs
was carried out at each temperature to calculate the time-independent
pair correlation functions discussed below. 







§ ACKNOWLEDGEMENT
 Sandia National Laboratories is a
multi-mission laboratory managed and operated by National Technology and
Engineering Solutions of Sandia, LLC., a wholly owned subsidiary of
Honeywell International, Inc., for the U.S. Department of Energy's
National Nuclear Security Administration under contract DE-NA-0003525.
This work is supported by the Assistant Secretary for Energy Efficiency
and Renewable Energy, Office of Vehicle Technologies of the U.S.
Department of Energy under Contract No. DE-AC02-05CH11231, Subcontract
No. 7060634 under the Advanced Batteries Materials Research (BMR)
Program. This work was performed, in part, at the Center for Integrated
Nanotechnologies (CINT), an Office of Science User Facility operated for
the U.S. DOE's Office of Science by Los Alamos National Laboratory
(Contract DE-AC52-06NA25296) and SNL.



§ AUTHOR CONTRIBUTION STATEMENT

AM and MIC performed simulations, analysis and prepared all the figures. LRP and SBR wrote the manuscript and all authors reviewed the manuscript. 




§ ADDITIONAL INFORMATION

The author(s) declare no competing financial interests.



§ SUPPLEMENTARY INFORMATION

Diffusion constants, conductivity and a topology file for EDC are reported as  supplementary information. 




10
url<#>1urlprefixURL doiprefixDOI 


NatureSEI2015
authorPiper, D. M. et al.
journaltitleStable silicon-ionic liquid
  interface for next-generation lithium-ion batteries.
Nature Comm. volume6
  (year2015).

identify_edc
authorZhuang, G., authorXu, K., authorYang,
  H., authorJow, T. & authorRoss, P. J.
journaltitleLithium ethylene dicarbonate
  identified as the primary product of chemical and electrochemical reduction
  of EC in 1.2 M LiPF6/EC:EMC electrolyte.
J. Phys. Chem. B
  volume109, pages17567–73
  (year2005).

Zorba:2012hc
authorZorba, V., authorSyzdek, J.,
  authorMao, X., authorRusso, R. E. &
  authorKostecki, R.
journaltitleUltrafast laser induced breakdown
  spectroscopy of electrode/electrolyte interfaces.
Applied Physics Letters
  volume100, pages234101
  (year2012).

Xu:2007ie
authorXu, K., authorLam, Y., authorZhang,
  S. S., authorJow, T. R. & authorCurtis, T. B.
journaltitleSolvation Sheath of Li +in
  Nonaqueous Electrolytes and Its Implication of Graphite/Electrolyte Interface
  Chemistry.
The Journal of Physical Chemistry C
  volume111, pages7411–7421
  (year2007).

Leung:2016jq
authorLeung, K., authorSoto, F.,
  authorHankins, K., authorBalbuena, P. B. &
  authorHarrison, K. L.
journaltitleStability of Solid Electrolyte
  Interphase Components on Lithium Metal and Reactive Anode Material
  Surfaces.
The Journal of Physical Chemistry C
  volume120, pages6302–6313
  (year2016).

BalbuenaP
authorMartinez de la Hoz, J. M., authorSoto, F. A. &
  authorBalbuena, P. B.
journaltitleEffect of the Electrolyte
  Composition on SEI Reactions at Si Anodes of Li-Ion Batteries.
J. Phys. Chem. C
  volume119, pages7060–68
  (year2015).

LeungK
authorLeung, K. & authorBudzien, J. L.
journaltitleAb initio molecular dynamics
  simulations of the initial stages of solid–electrolyte interphase formation
  on lithium ion battery graphitic anodes.
Phys. Chem. Chem. Phys.
  volume12, pages6583–6586
  (year2010).

Borodin:2016gm
authorBorodin, O. et al.
journaltitleCompetitive lithium solvation of
  linear and cyclic carbonates from quantum chemistry.
Physical Chemistry Chemical Physics
  volume18, pages164–175
  (year2016).

Kim:2011dn
authorKim, S.-P., authorDuin, A. C. T. v. &
  authorShenoy, V. B.
journaltitleEffect of electrolytes on the
  structure and evolution of the solid electrolyte interphase (SEI) in Li-ion
  batteries: A molecular dynamics study.
Journal of Power Sources
  volume196, pages8590–8597
  (year2011).

Borodin:2013bq
authorBorodin, O., authorZhuang, G. V.,
  authorRoss, P. N. & authorXu, K.
journaltitleMolecular Dynamics Simulations and
  Experimental Study of Lithium Ion Transport in Dilithium Ethylene
  Dicarbonate.
J. Phys. Chem. C
  volume117, pages7433–7444
  (year2013).

Bedrov:2017ci
authorBedrov, D., authorBorodin, O. &
  authorHooper, J. B.
journaltitleLi+ Transport and Mechanical
  Properties of Model Solid Electrolyte Interphases (SEI): Insight from
  Atomistic Molecular Dynamics Simulations.
The Journal of Physical Chemistry C
  volume121, pages16098–16109
  (year2017).

Jorn:2013kq
authorJorn, R., authorKumar, R.,
  authorAbraham, D. P. & authorVoth, G. A.
journaltitleAtomistic Modeling of the
  Electrode–Electrolyte Interface in Li-Ion Energy Storage Systems:
  Electrolyte Structuring.
The Journal of Physical Chemistry C
  volume117, pages3747–3761
  (year2013).

muralidharan2017molecular
authorMuralidharan, A., authorChaudhari, M.,
  authorRempe, S. & authorPratt, L. R.
journaltitleMolecular dynamics simulations of
  lithium ion transport through a model solid electrolyte interphase (sei)
  layer.
ECS Transactions
  volume77, pages1155–1162
  (year2017).

Borodin:2014gq
authorBorodin, O. & authorBedrov, D.
journaltitleInterfacial Structure and Dynamics
  of the Lithium Alkyl Dicarbonate SEI Components in Contact with the Lithium
  Battery Electrolyte.
J. Phys. Chem. C
  volume118, pages18362–18371
  (year2014).

Vorselaars:2007dg
authorVorselaars, B., authorLyulin, A. V.,
  authorKaratasos, K. & authorMichels, M. A. J.
journaltitleNon-Gaussian nature of glassy
  dynamics by cage to cage motion.
Phys. Rev. E volume75,
  pages011504 (year2007).

Stecki:1970df
authorStecki, J. & authorNarbutowicz, M. A.
journaltitleA generalization of vineyard's
  convolution approximation.
Chem. Phys. Letts.
  volume5, pages345–346
  (year1970).

YeomansReyna:2000jd
authorYeomans-Reyna, L., authorAcuña-Campa, H. &
  authorMedina-Noyola, M.
journaltitleVineyard-like approximations for
  colloid dynamics.
Phys. Rev. E volume62,
  pages3395–3403 (year2000).

g09
authorFrisch, M. J. et al.
titleGaussian 09 Revision A.1.
noteGaussian Inc. Wallingford CT 2009.

Rempe:1998
authorRempe, S. B. & authorJónsson, H.
journaltitleA computational exercise
  illustrating molecular vibrations and normal modes.
Chem. Ed. volume3,
  pages1–17 (year1998).

Gromacs
authorVan Der Spoel, D. et al.
journaltitleGROMACS: Fast, Flexible, and
  Free.
J. Comput. Chem.
  volume26, pages1701–1718
  (year2005).

oplsaa
authorJorgensen, W. L. & authorMaxwell, D. S.
journaltitleDevelopment and Testing of the
  OPLS All-atom Force Field on Conformational Energetics and Properties of
  Organic Liquids.
J. Am. Chem. Soc.
  volume118, pages11225–11236
  (year1996).

jctc:2016
authorChaudhari, M. I. et al.
journaltitle Scaling Atomic Partial Charges of
  Carbonate Solvents for Lithium ion (Li^+) Solvation and Diffusion.
J. Chem. Theory Comp.
  volume12, pages5709–5718
  (year2016).

Nose
authorNosé, S.
journaltitleA Molecular Dynamics Method for
  Simulations in the Canonical Ensemble.
Mol. Phys. volume52,
  pages255–268 (year1984).

Hoover
authorHoover, W. G.
journaltitleCanonical Dynamics: Equilibrium
  Phase-space Distributions.
Phys. Rev. A volume31,
  pages1695–1697 (year1985).

barostat
authorParrinello, M. & authorRahman, A.
journaltitlePolymorphic Transitions in Single
  Crystals: A New Molecular Dynamics Method.
J. App. Phys. volume52,
  pages7182–7190 (year1981).


	

