Nonequilibrium Electron-Coupled Lithium Ion Diffusion in LiFePO4

Subscriber access provided by ECU Libraries

Article 4

Nonequilibrium Electron Coupled Lithium Ion Diffusion in LiFePO: Nonadiabatic Dynamics with Multi-State Trajectory Approach Guohua Tao J. Phys. Chem. C, Just Accepted Manuscript • DOI: 10.1021/acs.jpcc.5b12676 • Publication Date (Web): 18 Mar 2016 Downloaded from http://pubs.acs.org on March 22, 2016

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

The Journal of Physical Chemistry C is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 45

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Nonequilibrium Electron Coupled Lithium Ion Diffusion in LiFePO4: Nonadiabatic Dynamics with Multi-State Trajectory Approach Guohua Tao*

School of Advanced Materials, Peking University Shenzhen Graduate School, Shenzhen, China 518055

Shenzhen Key Laboratory of New Energy Materials by Design, Peking University, Shenzhen, China 518055

Corresponding Author * To whom correspondence should be addressed. [email protected]

1 / 45

ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

ABSTRACT

The recently developed multi-state trajectory approach is implemented to nonadiabatic molecular dynamics simulations on electron coupled Li ion diffusion in LiFePO4, a promising cathode material for lithium ion batteries. The approach is demonstrated to be highly efficient and numerically stable in treating realistic complex molecular systems including tens of electronic states and hundreds of particles, and to be capable of providing inspiring molecular pictures via a multi-state representation on nonadiabatic quantum dynamics, such as Li ion diffusion and its correlation with electron transport in LiFePO4. The results from our real time nonequilibrium simulations indicate that the electron coupled Li ion diffusion is significantly enhanced upon the formation of small polaron and nonadiabatic couplings further facilitate Li ion diffusion. Therefore we suggest that nonequilibrium diffusion and nonadiabatic couplings could also be responsible for the wide range of Li ion diffusion coefficients obtained from various theoretical calculations and experimental measurements. Our findings may help understand the molecular mechanisms of Li ion diffusion in LiFePO4, and design new cathode materials of lithium ion batteries.

2 / 45


Page 2 of 45

Page 3 of 45

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


I.

Introduction LiFePO4 is a promising cathode material which would promote the large scale

applications of lithium ion batteries.1-6 Although its poor rate capability due to the low ion mobility and electronic mobility7,8 was improved by a variety of treatments, such as chemical doping,9 carbon coating,10 and modifying particles to nanoscales,11 the microscopic mechanism of the intrinsic lithium (Li) ion diffusion and electron transport in LiFePO4 is not fully understood yet. Theoretical and experimental evidences agree on that Li ion diffusion in pure LiFePO4 is basically one dimensional (along the b direction).12-14 However, previous theoretical calculations12,15-21 predicted Li ion diffusion much faster (~4-6 orders of magnitude larger) than the experimental values (D~10-13-10-15cm2s-1).22-26 Early studies also debate on which determines the electrochemical performance of LiFePO4. Morgan et al.12 calculated the activation energy Ea of Li ion diffusion in LiFePO4 to be 0.27 eV (which corresponds to D~10-9cm2s-1) and suggested that electronic conduction rather than the ion mobility is the limiting factor. Amin et al.24,25 measured the electronic mobility, ion mobility, and lithium chemical diffusivity in single crystals of LiFePO4 by using the DC and AC impedance spectrum. They found that the electronic mobility dominates, being orders of magnitude greater than the ion mobility. There are only few works to investigate the coupling of Li ion and electron transport. Maxisch et al.27 predicted low activation energy for either free polaron or free Li vacancy in LiFePO4, but large binding energy of the two, i.e. 500 meV, which implies that the coupling impedes the transport of both free carriers. Furthermore the 3 / 45



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

strong binding causes the interpretation of the conductivity measurements difficult. Ellis et al.28 proposed that the electron delocalization is coupled with Li disorder based on their Mossbauer and XRD experiments, and also predicted a higher Ea for Li ion diffusion when it is coupled with the small polaron hopping of electron transport. Zaghib et al.29 suggested that the small polaron consists of Fe3+/Li vacancy instead of the e-Li+ excitonic pair,27,28 while the Ea for electron transport that is associated with the intersite transition of Fe2+ -> Fe3+ is not related to the Li+ diffusion although the two are on the same order of magnitude. There is no doubt that previous work27-29 provide valuable physical insights into the coupled Li ion and electron transportation, however an accurate dynamical description is called for the quantitative analysis on Li ion and electron transport in LiFePO4 to fully understand the underlying microscopic mechanisms. Although the coupled e-Li+ diffusion was investigated in TiO2 by using molecular dynamics (MD) simulations combined with phenomenological Marcus theory,30 to the best of our knowledge, real time molecular dynamics simulation incorporating direct nuclear-electronic coupled nonadiabatic dynamics has not been carried out yet, which would definitely provide accurate dynamical information and valuable molecular pictures on the electron coupled Li ion diffusion in LiFePO4. The coupled nuclear-electronic dynamics is however extremely complicated and there are only a few theoretical methodologies available for describing related processes effectively. Ehrenfest method is a mean field approximation and may causes unphysical mixing of multiple electronic states.31 Surface hopping and its 4 / 45


Page 4 of 45

Page 5 of 45

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


variants32,33 are widely applied to simulate the nonadiabatic dynamics involving multiple electronic states due to its simplicity. However it is basically an “ad hoc” method that requires additional information on the treatment of electronic hopping. Alternatively a semiclassical formalism that based on the Meyer-Miller (MM) model34 for the electronic degrees of freedom can provide a consistent and rigorous description on the coupled nuclear-electronic dynamics.31 For example, Miller’s semiclassical methodology can describe nuclear quantum coherence quite accurately,35 which is usually missed in other methods. Recently the symmetric quasiclassical/Meyer-Miller (SQC/MM) approach36-39 has been developed to describe nonadiabatic dynamics in complex molecular systems, such as singlet fission in organic semiconductors.40-42 To perform accurate and efficient nonadiabatic dynamics simulations on realistic molecular systems, recently we proposed a multi-state trajectory (MST) method.43 The MST method associates each state with an individual trajectory therefore may have advantages of describing nuclear motions on different state faithfully. When making the quasiclassical approximation assuming that the nuclear trajectory is classical or the Born-Oppenheimer-like approximation with nuclear forces independent of electronic populations, the resulting MST-QC and MST-BO algorithms are very stable and easy to be implemented to complex systems. It has been shown that the MST-QC and MST-BO methods can provide results in reasonably good agreement with exact quantum calculations on several benchmark systems such as the spin-boson problem. In this work we implement the MST 5 / 45



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 45

approach (MST-QC and MST-BO) to the realistic molecular systems of Li+ diffusion in LiFePO4, with atomic potentials applied. Combined with the conventional equilibrium and nonequlibrium MD simulations, our results provide a dynamical perspective on the understanding of the electron coupled Li ion diffusion in LiFePO4, and on how to reconcile the difference between the theoretical and experimental results. II.

Theory

a. Multi-state trajectory approach To

describe a

nuclear-electronic coupled system,

we

take the

MM

nuclear-electronic Hamiltonian:31,34

(

)

n n 1 2 xi + pi2 − 1 H ii (P, Q) + ∑ ∑ (xi x j + pi p j )H ij (P, Q) i =1 2 i =1 j =i +1 n

H (x, p, Q, P ) = ∑

(1) where n is the number of electronic states, Hii and Hij are the diagonal and off diagonal diabatic Hamiltonian matrix elements to be determined by the nuclear coordinates and momenta Q, P, and xi, pi are the Cartesian variables for the classical mapping of the i-th electronic degree of freedom (DOF), respectively. Classical trajectories for the nuclear and electronic DOF are propagated with respect to time following Hamilton’s equations & = ∂H , P& = − ∂H , x& = ∂H , p& = − ∂H , Q ∂P ∂Q ∂p ∂x

(2)

with the Hamiltonian given by Eq 1, which is a set of highly nonlinear equations and very difficult to be solved numerically. Here we introduce a multi-state trajectory 6 / 45


Page 7 of 45

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


assuming each state is associated with an individual nuclear trajectory43 and adopt the approximations as the follows: (a) Assume the nuclei propagate quasi-classically (MST-QC) or on the corresponding Born-Oppenheimer potential surface (MST-BO);43 (b) The electron dynamics is updated at a shorter time step while the nuclear DOFs are frozen until a normal nuclear time step is reached. The electronic dynamics is treated in the same way as in the SQC method.36,37 And a symmetrical binning window function is applied to determine to which electronic state the MST trajectory belong, i.e.

Wk (nk , N ) =

1  ∆n  h − nk − N  , ∆n  2 

(3)

where h(z) is the Heaviside function, ∆n = 2γ , γ is a parameter, and N (0 or 1) is the occupying number. The nonadiabatic dynamical properties therefore may be evaluated by the ensemble average of the MST trajectories. For example, assuming the initial state of the system is the state i, the time-dependent electronic population of state k is given by the follows:40 Wi (ni (0),1)∏l ≠ i Wl ( nl (0),0)Wk ( nk (t ),1)∏l ≠ k Wl (nl (t ),0) n

Pk (t ) =

∑

n k

n

Wi ( ni (0),1)∏l ≠i Wl ( nl (0),0)Wk ( nk (t ),1)∏l ≠ k Wl (nl (t ),0) n

n

.

(4)

The MST algorithm appears to be stable and efficient.

b. Nuclear-electronic coupled Hamiltonian The energy levels of each exciton state are given by Eii = Eii0 + Eiiint , where

Eii0 , Eiiint are the static zero order energy and the dynamical intra-pair energy of the state i, respectively. For the ferric-ferrous electron transfer in this work, Eii0 are the 7 / 45



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

same for all state i, i.e. Eii0 = 0 , and we take Eiiint = −

Page 8 of 45

N

q Fe3+ ( i ) q J

J ≠ Fe3+

rFe3+ ( i ) − rJ

∑

where the sum

is taken over all other ions J, N is the number of ions, and q, r are the atomic charge and coordinate. The inter-state couplings are electronic couplings between Fe3+ and Fe2+, and may be approximated by the interaction of “transition charge”,44 i.e.

Eij = −

q * Fe3+ ( i ) q *Fe 2+ ( j ) rFe3+ ( i ) − rFe2 + ( j )

.

(5)

To account for the screening effect, we approximate the Columbic coupling of the transition charges by that of screened partial atomic charges:

Eij = −

qFe3+ (i ) qFe2+ ( j ) rFe3+ (i ) − rFe2+ ( j )

(

erfc α rFe3+ (i ) − rFe2+ ( j )

)

,

(6)

where erfc() is the complementary error function, and α is a screening factor. An example of the energy levels and couplings are listed in Table S1.

III.

Computational Details The simulation box consists of a 2x3x4 supercell (20.66Åx18.03Åx18.77Å) of

LiFePO4, which includes 672 atoms (96 formula units) in total. The force field is taken from references.17,45 Columbic interactions are calculated by using the Ewald summation, and the periodical boundary condition (PBC) is applied. We carry out the NVE simulation with the velocity Verlet algorithm for the integration of nuclear equations of motion. The system is equilibrated for about 600 ps after rescaling the temperature to the desired values. The subsequent 6 ns trajectories are used for the equilibrium calculations. The numerical stability is checked and the ∆E/E Fe3+ + Fe2+. The nonadiabatic MST simulations are performed to describe electron transfer and electron coupled Li ion diffusion in nonequilibrium dynamics in LiFePO4. For the MST-QC calculations in this work, we consider the closest 11 Fe2+ neighbors to the central Fe3+ as possible electron transfer sites, i.e. 12 electronic states included in total. The simplified MST-BO method may treat more complex systems such as the 24-state system, an illustration for which is shown in Figure 1. 9 / 45



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 45

Figure 1. A sample configuration of LiFePO4 with a small polaron in nonadiabatic MD simulations. Fe3+ (blue, not to scale), Li vacancy (light blue, removed), Fe2+ (golden, labeled for different polaron states for which Fe2+ is oxidized to Fe3+), Li (green), P (purple), and O (red).

The nonadiabatic mean square displacement (MSD) function is calculated by averaging the MSD function for each individual electronic state weighted by the time-dependent state population, i.e. n

MSD(t ) = ∑ Pk (t ) k =1

1 N

N

∑ [r

ki

(t ) − rki (0)]

2

,

(7)

i =1

where rki(t) is the coordinate of the ith Li ion at time t for the electronic state k, Pk(t) is the time dependent state population, N and n are the number of ions and electronic states, respectively. The overbar here denotes the nonequilibrium average. The nonadiabatic electron MSD function is given as the follows:

10 / 45


Page 11 of 45

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


[

n

MSD(t ) = ∑ Pk (t ) rk ,Fe( III ) (t ) − rk ,Fe( III ) (0)

]

2

,

(8)

k =1

here rk ,Fe ( III ) (t ) is the time dependent coordinate of the Fe3+ ion for the electronic state k. The diffusion constant is given by the Einstein equation:

D=

d [ MSD(t )] 1 lim , t → ∞ 2F dt

(9)

where F is the dimension of diffusion motions. In this work, we obtain the diffusion constant from the slope of the linear fit to MSD functions for t >1 ps. The activation energy Ea is estimated by the following relation: D = D0 exp[− Ea / k BT ] ,

(10)

here kB is the Boltzmann constant, T the temperature and D0 a constant. The pair radial distribution function (RDF) for particle A and B is given by

g A− B ( r ) =

1

ρ AρB

N A NB

∑∑ δ (r

ij

i

− r)

,

(11)

j

where ρ , N, rij are the number density, number of particles and the distance between particle i and j, respectively. In the nonequilibrium nonadiabatic

dynamics simulations, 96 nuclear

configurations are used, for each of which 20 electronic configurations are sampled for the final statistical averaging for correlation functions if not specified otherwise. A multiple time step treatment for electronic and nuclear DOFs is used, i.e. dte = 2.0 a.u. and dtn = 50.0 a.u. The screening factor α=0.5 is picked up such that the interstate electronic couplings between nearest-neighbor Fe2+/Fe3+ pairs are comparable with those in the FeO system.46 Molecular visualizations were generated using the VESTA program.47 11 / 45



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

IV.

Results and Discussions

a. Lithium ion diffusion First the equilibrium MSD functions of Li ions in pure bulk LiFePO4 have been calculated (see Figure S1). Starting from zero, the very short-time rise in the MSD functions represents the translational libration of Li ions in local potential wells. The magnitude of the short-time rise is temperature dependent, indicating that Li ions can move around in a larger domain in space as temperature increases. The long-time part of the MSD functions is contributed by Li+ diffusion. Apparently Li+ diffusion is preferentially one dimensional in the b direction, which is consistent with previous studies.12-14

Figure 2. MSD function of Li+ diffusion along different directions at various temperatures. Plotted are the results from the equilibrium simulation (green), nonequilibrium single state (blue) and nonadiabatic (MST-QC, black; MST-BO, red) simulation upon the formation of a small polaron. X 12 / 45


Page 12 of 45

Page 13 of 45

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(dotted), Y (solid), and Z (dashed) denote the directions along [100], [010], and [001], respectively. (a) 1000 K; (b) 900 K; (c) 800 K; and (d) 700 K. The screening factor α = 0.5 and γ = 0.5 in the nonadiabatic simulations.

Figure 2 illustrates the effect of the formation of small polaron on Li+ diffusion. Plotted are the results at different temperatures T = 1000, 900, 800, and 700 K (see Figure 2a-d, respectively) for the equilibrium, nonequilibrium single state (the initial state relaxation with no electronic couplings), and nonequilibrium nonadiabatic simulations. After the small hole polaron/Li vacancy is generated, the short-time dynamics is nearly identical to the equilibrium ones, which implies that the inertial motions of Li+ in local environment does not change substantially for the first ~200fs, and the linear response governs the nonequilibrium relaxation.48,49 On the other hand, the long-time Li+ diffusion is greatly accelerated even if no electron hopping is involved. The enhancement on the Li+ diffusion appears mainly one dimensional along the y axis. Introducing nonadiabatic couplings among multiple states that correspond to that the Fe3+ located on different iron ions could further speed up Li+ diffusion, which is in contrast to the conclusions from previous studies on that the nuclear-electronic coupling impedes Li+ diffusion.27,28 Note that Maxisch et al.27 performed their static DFT calculations on the adiabatic potential surfaces, while Ellis et al.28 attributed the coincidence of the onset of small polaron hopping with the Li disorder at elevated temperatures upon the formation of solid solution without specifying the nature of the coupled electron-lithium dynamics. As temperature decreases, the MSD of Li ion in all directions reduce and diffusion slows down (note that the different scales used for Figure 2a-b from that for 2c-d). The more 13 / 45



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

approximate MST-BO method predicts a slightly slower Li+ diffusion than what the MST-QC does. The enhancement of Li+ diffusion due to the electronic-nuclear coupling may be illustrated as the follows: once a Li vacancy is generated together with the formation of a small hole polaron, nearby Li ions are more easily to move into the vacancy in comparison with the otherwise filled site, resulting in facile diffusion. This is supported by our MD simulations. Figure 3 displays the snapshots from a representative individual trajectory for the state 1 (initial state) in the MST-BO simulation on the 12-state system. At the time t = 0, a small polaron is formed with one of Fe2+ replaced by Fe3+ (blue in Figure 3a) and its closest neighboring Li ion removed to leave a vacancy (represented by a white sphere). The subsequent dynamics is mapped out at times t = 0.121, 0.605, 1.21, 6.05, and 12.1 ps following the nonequilibrium small polaron formation (see Figure 3b-3f, respectively). Particularly Figure 3c indicates that Li ion may cross over adjacent Li cages (see the pink one on the bottom is now in between two adjacent equilibrium Li sites). In case of that nonadiabatic transitions occur, a multi-state picture has to be taken, i.e. individual trajectories for other states also make contributions to the dynamics, and the overall behavior is an ensemble average of individual trajectories associated with each involved state, the system hops from the initial state to another one (see Figure 4) during which electron transfers from a nearby Fe2+ to the original Fe3+ (see for example Figure 3a vs. Figure 4a) resulting in a switch of Fe2+ (in brown polyhedral) and Fe3+ (in blue) in coupled with Li ion diffusion. Figure 4a-f present the 14 / 45


Page 14 of 45

Page 15 of 45

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


snapshots at the corresponding nonequlibrium time instants of an individual trajectory for the state 2 in the same MST simulation as those shown in Figure 3. The electron coupled Li ion diffusion is in general faster than the single-state diffusion since the small polaron formed in other states span a wide range on the separation of the Fe3+/Li vacancy, which may induce long-range Li ion diffusion in the nonequilibrium relaxation of molecular structures. For example, a nearby Li attempts to move into the vacancy in state 1 (see Figure 3c) but no successful crossover has been made at least until 12.1 ps (Figure 3f) due to the strong binding energy of the Fe3+/Li vacancy pair. Fe3+ in the other state may prefer a close Li vacancy again due to the larger binding energy, pushing the nearby Li ion away into the initial vacancy site. Figure 4 captures such a molecular event: at t = 1.21 ps (Figure 4d), the topmost Li ion in pink is trying to move upward to create a vacancy close to Fe3+ in this state; at t = 6.05 ps (Figure 4e), this Li ion has successfully moved into its neighboring cage where the black one occupied (note that here PBC is applied) and at the same time, the black Li ion has been pushed into the original Li vacancy and remains there for a while, at least until t = 12.1 ps (Figure 4f) whereas the topmost pink Li ion seems switching back and forth around the new vacancy close to Fe3+. Therefore the large binding energy of the small polaron and Li vacancy, which was predicted from the previous static DFT calculations,27 could actually be the driving force for facile Li+ diffusion, i.e. the fast electron transfer promotes Li ion diffusion due to strong electronic-nuclear couplings in a dynamical picture.

15 / 45



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 3. Snapshots of a representative trajectory for state 1 in the MST-BO simulation. The initial nonequilibrium configuration consists of Fe3+ (blue) and Li vacancy (white), and other Li ions in the same b channel with the Li vacancy are in pink or black. Color schemes for other particles are O (red), P (purple), Fe2+ (golden) and Li+ (green). Shown are the a-b plane projections at different times following the small polaron formation, namely (a) 0 ps; (b) 0.121 ps; (c) 0.605 ps; (d) 1.21 ps; (e) 6.05 ps; and (e) 12.1 ps. T = 700K.

16 / 45


Page 16 of 45

Page 17 of 45

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 4. Snapshots of a representative trajectory for state 2 in the MST-BO simulation. Notations are the same as for Figure 3, except for the brown polyhedral which represents for the location of Fe3+ corresponding to the state 1 in the initial nonequilibrium configuration.

17 / 45



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 5. Snapshots of a representative trajectory for state 1 in the MST-QC simulation. Initial configuration and notations are the same as for Figure 3.

18 / 45


Page 18 of 45

Page 19 of 45

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 6. Snapshots of a representative trajectory for state 2 in the MST-QC simulation. Initial configuration and notations are the same as for Figure 4.

19 / 45



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

In the MST-BO simulations, nuclear dynamics for each state is not determined by electronic variables, but electronic dynamics is modulated by nuclear motions, therefore the electronic-nuclear couplings are not fully accounted and the enhancement of Li ion diffusion is attributed mainly to the introduction of multi-state contributions. By contrast, nuclear and electronic dynamics are directly coupled with each other in the more accurate MST-QC approach, and the enhancement of Li ion diffusion due to the electron transfer may be illustrated by the snapshots of representative MST-QC trajectories (see Figure 5 and 6) starting from the same initial conditions as in the MST-BO simulation (Figure 3 and 4, respectively). Figure 5a-f take the snapshots of the individual trajectory for state 1 in the MST-QC simulation on the nonequilibrium dynamics. Molecular events for occupying the Li vacancy involving Li ions from both sides in the b channel can be seen in Figure 5e and 5f. Figure 6a-f display the dynamics for state 2 at the same selective time instants, which reveals that the new vacancy may not be bounded tightly to the Fe3+ (see Figure 6e and 6f), indicating a long-time nonequlibrium situation. The inclusion of quantum back-reaction from electronic dynamics in the MST-QC simulations therefore fully accounts for the nuclear electron couplings, which allows more freedom for Li ion diffusion. The MST approach therefore provides a clear semiclassical view on nonadiabatic quantum dynamics in terms of an MST representation (see Figure 7). Although direct quantum trajectory is hard to be presented straightforwardly, nonadiabatic systems involving multiple states may be described by the MST, i.e. multiple classical-like 20 / 45


Page 20 of 45

Page 21 of 45

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


trajectories, for example, nuclear dynamics shown in Figure 3-6. Together with the time dependent electronic population functions, a semiclassical molecular picture on nonadiabatic dynamics may be obtained from the real time MST simulations, which would be extremely helpful for realistic applications such as materials design.

Figure 7. A schematic presentation for the MST representation for nonadiabatic quantum dynamics for a three-state model. All individual trajectories (represented by wave packets) associated with its corresponding state evolve simultaneously. The wave packet in solid color indicates that the system is on the corresponding state, while those for other states are semi-transparent. The initial state is on state 1, followed by a transition to state 2, and then to state 3. To better understand the time evolution of molecular structures during the nonequiliubrium relaxation following the small polaron formation, we plot the RDFs of Fe3+- Li+, Fe3+- Fe2+, and Fe3+- O for state 1 (initial state) and state 2, shown in Figure 8 and 9, respectively. In comparison with the equilibrium results, the creation of Li vacancy is clearly seen in Figure 8a, and the subsequent evolution results in a partial refill of this vacancy. Note that the initial increase in the first peak of the RDF of Fe3+- Li+ at t = 121 fs in Figure 9a is largely due to the repulsion between the new formed Fe3+ and its neighboring Li+, which is followed by the long time genuine refill 21 / 45



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

of Li vacancy and partial recovery of the initial expansion of the other Li neighbors around Fe3+. In the RDFs of Fe3+- O, the short time shift of the first peak reflects the strong attractions between Fe3+ and its neighboring O atoms, as shown in Figure 8b and 9b. Whereas the small decease in the second peak is attributed to the broken Li-O bonds surrounding the Li vacancy due to the removal of Li+, therefore these former Li-connecting O atoms effectively move further away from the Fe3+. On the other hand, it is accompanied by the contraction of the Fe-O bonds for these O atoms and the O-connecting Fe2+ getting closer to Fe3+, indicated by the small left shift of the third peak in the RDF of Fe3+- Fe2+ (see Figure 8c and 9c). These structure changes alter the local environment for Li ion diffusion, which induces nonequilibrium relaxations in structures and dynamics.

Figure 8. Time evolution of the molecular structure during the nonequilibrium relaxation upon the 22 / 45


Page 22 of 45

Page 23 of 45

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


formation of small polaron in LiFePO4. Shown are pair radial distribution functions of the iron ion with other particles calculated at different times (t0) for the diabatic state 1, together with the corresponding equilibrium results. T = 700K.

Figure 9. The same as Figure 8 but for state 2.

Figure 10. Lithium ion diffusion coefficients for the [010] direction. Plotted are the results obtained from the equilibrium simulation (black squares), the nonequilibrium single state (red squares) and nonadiabatic (blue squares) simulation upon the formation of small polaron. Experimental values are plotted for comparison, taken from Ref 22 (golden spheres), Ref 23 (pink spheres) and Ref 24, 25 (green spheres), respectively.

23 / 45



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 10 displays the Li+ diffusion coefficients, which are calculated from the corresponding MSD functions by the Einstein equation, and the activation energy for diffusion obtained from a linear fit to the temperature dependent diffusion constants in the semi-log plot. The Li+ diffusion activation energy in pure bulk LiFePO4 is 811 meV, and the estimated diffusion coefficient at room temperature is about 10-16 cm2s-1 from the extrapolation. When small polaron is introduced, Li+ diffusion is enhanced, and the nonadiabatic couplings make it even faster. The corresponding activation energies are 304 meV, and 229 meV (MST-QC), for the diabatic, and nonadiabatic electron coupled Li ion diffusion, respectively. These values are in good agreement with previous static DFT calculations by Dathar et al.16 (290 meV) and by Maxisch et al.27 (215 meV). For comparison, representative experimental results from Franger et al.,22 Prosini et al.,23 and from Amin et al.24,25 are also plotted in Figure 10. They are well distributed in the region set by the extrapolation lines of our equilibrium and nonequilibrium nonadiabatic simulation results. Many other reported experimental measurements also fall into this region. The existence of anti-site defects,12,14,50 interfaces,16 and partial delithiathed phases11,14,29 were suggested to be responsible for the significant differences between theoretical and experimental results for Li+ diffusion. Other complicated factors normally exist in experiments, such as external field, random distributed crystalline orientations and so on. Here our MD simulations demonstrate that even without these factors, Li+ diffusion in LiFePO4 may be modulated significantly by nonequilibrium relaxation and nonadiabatic couplings. Therefore our results provide unambiguous 24 / 45


Page 24 of 45

Page 25 of 45

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


references to interpret experimental results and to understand the underlying mechanisms.

Figure 11. Convergence check for the MST nonadiabatic dynamics simulations. (a) Dependence on the number of electronic (Ne) and nuclear (Nn) trajectories for the MST-QC simulations. (b) Dependence on the number of states involved in the MST-BO simulations. T = 700 K, α = 0.5, and γ = 0.5. X (dotted), Y (solid), and Z (dashed) denote the directions along [100], [010], and [001], respectively.

The convergence check is performed for the nonadiabatic dynamics simulations using the MST approach, and the results are shown in Figure 11. When increasing the number of electronic trajectories from 20 to 100, or the number of nuclear trajectories from 96 to 480, the MSD functions of Li ion are virtually unchanged for both MST-QC (Figure 11a) and MST-BO (not shown) simulations. However the increase in the number of electronic trajectories indeed greatly improves the convergence of electronic population and electronic MSD function, therefore we use 100 electronic 25 / 45



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

trajectories (times 96 nuclear trajectories) in the following corresponding calculations. The dependence on the number of electronic states has been examined for the MST-BO simulations (Figure 11b). During the time window considered, only weak dependence in the y direction is noticed. We do not do the same check for the MST-QC approach since the results for the system with more states are not available yet, and this will be done in the future work. The dependence of nonadiabatic Li ion diffusion on the choice of the screening factor α and the parameter γ has also been examined. Figure 12 compares the MST-QC results for the 12-state system at temperatures T = 1000 K (Figure 12a) and T = 700 K (Figure 12b), respectively, by using α = 0.5 and γ = 0.5 (from Figure 2), with those for a different value of γ, (γ = 0.366, the preferred one in the SQC calculations for model systems37) and with those for a different screening factor α = 0.64, for which the interstate couplings are about 10 times smaller. The nonadiabatic Li+ diffusion on the x and the z directions are largely insensitive to the choice of α and γ. For the y direction, the dependence of Li ion diffusion on γ is weak whereas the small electronic couplings (large α) results in slow Li+ diffusion. In particular at low temperature T = 700 K (see Figure 12b), the nonadiabatic calculations with large α are about the same as the single state results, which is consistent with our previous analysis that nonadiabatic couplings enhance Li+ diffusion. The nonadiabatic enhancement becomes more appreciable as temperature increases (see Figure 12a), which reflects the fact that the dynamics is nuclear-electron coupled.

26 / 45


Page 26 of 45

Page 27 of 45

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 12. MSD function of Li+ diffusion along different directions at various temperatures. Plotted are the nonadiabatic MST-QC results for the 12-state system. The parameters α = 0.5, γ = 0.5 (black, from Figure 2), α = 0.5, γ = 0.366 (red), and α = 0.64, γ = 0.5 (blue). The results for equilibrium (green), and nonequilibrium single state (golden) are also shown for comparison. X (dotted), Y (solid), and Z (dashed) denote the directions along [100], [010], and [001], respectively. (a) 1000 K; (b) 700 K.

Similar comparisons on the results obtained for different values of α and γ in the MST-BO calculations are made in Figure 13. For both high (T = 1000 K, Figure 13a) and low temperatures (T = 700 K, Figure 13b), the results for the x and the z directions are virtually the same irrespective to the changes in the screening factor and γ. Results for the Li+ diffusion along the y axis for different values of γ are nearly overlapped with each other, and the dependence on the screening factor α is much weaker than that seen in the MST-QC simulations.

27 / 45



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 13. MSD function of Li+ diffusion along different directions at various temperatures. Plotted are the nonadiabatic MST-BO results for the 24-state system with the parameters α = 0.5, γ = 0.5 (black, taken from Figure 2), α = 0.5, γ = 0.366 (red), and α = 0.64, γ = 0.5 (blue). X (dotted), Y (solid), and Z (dashed) denote the directions along [100], [010], and [001], respectively. (a) 1000 K; (b) 700 K.

We note that the high rate capability of cathode materials containing LiFePO4 has been attributed to the formation of a nonequilibrium solid solution phase51,52 (in contrast to the thermally equilibrium two-phase state of lithiated and delithiated phases) in a certain Li concentration regime (0.05

Nonequilibrium Electron-Coupled Lithium Ion Diffusion in LiFePO4

Recommend Documents