- Email: [email protected]
Ab initio and kinetic Monte Carlo simulation study of lithiation in crystalline and amorphous silicon
Journal of Power Sources 272 (2014) 1010e1017
Contents lists available at ScienceDirect
Journal of Power Sources journal homepage: www.elsevier.com/locate/jpowsour
Ab initio and kinetic Monte Carlo simulation study of lithiation in crystalline and amorphous silicon Janghyuk Moon a, Byeongchan Lee b, Maenghyo Cho a, Kyeongjae Cho a, c, * a
WCU Multiscale Mechanical Design Division, Department of Mechanical and Aerospace Engineering, Seoul National University, Seoul 151-742, Republic of Korea b Department of Mechanical Engineering, Kyung Hee University, Yongin 446-701, Republic of Korea c Department of Materials Science and Engineering and Department of Physics, The University of Texas at Dallas, Richardson, TX 75080, USA
h i g h l i g h t s DFT and KMC study of energetics and kinetics of c-Si and a-Si under Li insertion. Phase transformation of c-Si to a-Si under Li insertion via SieSi bond breaking. Li diffusion in a-Si predicted by environment dependent KMC method. KMC shows larger Li diffusivity in a-Si than c-Si at room temperature.
a r t i c l e i n f o
a b s t r a c t
Article history: Received 6 February 2014 Received in revised form 19 August 2014 Accepted 2 September 2014 Available online 16 September 2014
Energetics and kinetics of Li insertion into c-Si and a-Si systems are investigated using the density functional (DFT) theory calculations and kinetic Monte Carlo (KMC) simulations. DFT formation energies show the mechanism of phase separation between crystalline silicon and amorphous lithium silicide. Both crystalline and amorphous Si show similar trends in volume expansion and phase transition under lithiation, and kinetics of Li diffusion in bulk silicon (from DFT and KMC) shows a big difference between c-Si and a-Si. The Li migration barrier is 0.6 eV in c-Si, and quickly decreases to 0.4 eV under increasing Li concentration or Si volume expansion. To simulate Li diffusion in amorphous silicon using KMC, we have developed a formulation for environment dependent migration energy barriers of Li in a-Si using a volume dependent function. KMC simulations are performed for Li diffusion in both c-Si and a-Si, and the diffusion coefficient of Li in a-Si is an order of magnitude larger than in c-Si. These studies help to understand mechanisms of lithiation with atomic scale details and elucidate the phase separation between c-Si and lithium silicide. © 2014 Elsevier B.V. All rights reserved.
Keywords: Ab initio method Kinetic Monte Carlo simulation Silicon Li diffusion Phase separation
1. Introduction The lithiation of silicon has attracted a lot of research interest due to the possible application of nano-structured silicon as a high capacity anode material for Li-ion batteries. [1] Compared to graphite (372 mAh g1), silicon has a large theoretical charge capacity (4200 mAh g1) which is more than 10 times higher than that of graphite. However the lithiation and delithiation characteristics of Si are qualitatively different from those of graphitic
* Corresponding author. Department of Materials Science and Engineering and Department of Physics, The University of Texas at Dallas, Richardson, TX 75080, USA. E-mail address: [email protected] (K. Cho). http://dx.doi.org/10.1016/j.jpowsour.2014.09.004 0378-7753/© 2014 Elsevier B.V. All rights reserved.
carbon anode. In contrast to intercalation compounds (graphite anode and transition metal oxide cathode materials), Si anode undergoes phase transformation involving SieSi bond breaking as the Li concentration increases in LixSi. For x 1, diverse crystalline phases (LiSi, Li12Si7, Li7Si3, Li13Si4, Li15Si4, Li22Si5) are known to be thermodynamically stable. However, room temperature lithiation and delithiation studies have shown a crystalline-to-amorphous phase transition that takes place during the first charge cycle of crystalline Si anode. [2e4] Furthermore, the nonconventional phase change characteristics of the lithiation of Si introduces a volume changes up to 400%, an elastic softening, and a phase separation during insertion of Li [5e8]. The lithiation of Si results in substantial atomic structure changes, in contrast to transition metal oxides or graphite electrode materials where Li ion intercalates into existing layers with limited
J. Moon et al. / Journal of Power Sources 272 (2014) 1010e1017
structural distortions. Specifically, Si transforms into various crystalline LixSi phases at high temperature (400 C), but amorphous lithium silicides (LixSi) are formed during lithiation at room temperature. [9e11] The formation of crystalline or amorphous phases during the lithiation of Si is known to be associated with capacity fade. [5,11,12] Previous studies have attempted to explain the performance degradation of Si anode resulting from the complex phase change during lithiation. The strong coupling between diffusion and large mechanical deformation would induce the corresponding mechanical stress, [13e16] crack formation [11,17], and mechanical weakness due to the formation of weak LieSi bonds. [18e20] Even though these previous studies have explained certain aspects of Si anode degradation mechanisms, quantitative nature of atomic structure evolution during the lithiation of Si is not well understood yet. Specifically, the initial lithiation of crystalline silicon is known to produce amorphous LixSi compound, but the detailed nature of amorphous phase formation is not understood yet. Furthermore, how Li diffuses in a-LixSi is not well understood, either. In recent studies, [21,22] many investigations on failure mechanism of Si anode have focused on the issues of phase separation and diffusion kinetics of Li because they significantly affect the stress state at interfaces or lithiated domains during lithiation cycles (Fig. 1). Fig. 1b schematically shows a phase separation between LixSi and Si, and large interfacial stresses are expected to induce mechanical failure modes. Furthermore, the Li diffusion is a key parameter in determining how fast a battery can be cycled (i.e., power density). [2,9,15,23e25] Current experimental and theoretical reports on the failure mechanisms of LieSi systems are not fully consistent with each other. Specifically, several DFT-based studies of LieSi systems have been reported, focusing on Li behavior studies in bulk or nanostructure, atomic and electronic properties of LieSi structure based on high temperature crystalline compounds, or aspects of amorphous model structure generated by abinitio molecular dynamic (AIMD) simulation. Two-phase (Si and LixSi) separation (shown in Fig. 1b) is generally known as the key failure mechanism of Si anode which arises from a thermodynamic driving force, diffusion induced Li accumulation, and a high activation energy of phase transition on the SieLixSi interface. [26] However, a quantitative theoretical basis to explain the phase separation phenomenon is yet to be developed even though many authors have reported different aspect of it. It has been reported in experimental studies that the thermodynamic driving force of two-phase separation and activation energy of phase transition is high, but previous theoretical studies are not consistent with these finding. [27,28] Consequently, it is necessary to develop more detailed theoretical understanding based on DFT and kinetic Monte Carlo (KMC) simulation studies to elucidate the phase separation between c-Si and a-LixSi. Current DFT and KMC studies in this paper show complex phenomena of initial phase instability and Li accumulation due to slow Li diffusion in c-Si as we discuss later.
1011
Kang et al. have experimentally shown that crystalline LiSi phases are present during battery operation cycles, although primarily amorphous LixSi states are present for most Li composition during charging and discharging. [29] The maximum capacity of Si anode at room temperature has been attributed to the crystalline Li15Si4 phase. [30,31] In this work, we focus on initial stage of Li insertion into c-Si before the occurrence of crystalline LixSi (i.e., x < 1 in LixSi) and investigate the energetics and kinetics of initial lithiation. 2. Methodology 2.1. Computational methods The ab initio calculations are performed using Vienna ab initio simulation package (VASP) for the density functional theory (DFT) study of LixSi model systems within the local density approximation. [32e34] The projector augmented wave method is chosen to describe the interaction between ion core and valence electrons. [35] For Li atom, the Li-sv pseudo-potential treated the semi-core 1s states as valence electrons for accurate description of Liþ ionic states. A kepoint mesh in the Monkhorst-Pack scheme is set by 3 3 3 for 2 2 2 supercell of Si. The plane wave basis cut off energy is 500 eV. The structures are fully relaxed for both internal atomic coordinates and supercell shapes. The geometry optimizations are carried out until the forces on each atom become smaller than 0.01 eV Å1. For ab initio molecular dynamics simulation, the kepoint is sampled only at the gamma point. To build amorphous structures, the MD simulations are performed at 1200 K to randomize their initial crystalline configurations. A time step of 3.0 fs and overall simulation steps of 2000 are used. Insertion of Li atoms is done in such a way to achieve a homogenous distribution. At least five samples are generated at each Li concentration (x) in LixSi models. The Delaunary triangulation method (previously used to generate amorphous LieSi bulk) is chosen for Li insertion. [36] To study Li kinetics in Si environment, diffusion pathways and barriers were determined using the climbing-image nudged elastic band (NEB) method. Macroscopic diffusion coefficients of Li in c-Si and aSi are obtained using Arrhenius equation and kinetic Monte Carlo (KMC) simulations. The detailed KMC algorithm for simulating the time evolution follows the conventional procedure [37,38]. Results of the KMC simulations are statistically averaged to obtain macroscopic behavior. In this study, the standard KMC algorithm is used to simulate an activated random-walk process of Li interstitial defects in Si. Movements of the Li defects have followed the general KMC procedures: (i) Identify all possible events from the current atomic configuration. In this study, we consider first-nearest-neighbor (<3.5 Å) Li sites from sites of tetrahedral Si. (ii) Determine a list of the rates (ki) of all possible transition events in the system. The rates are proportional to exp(eEm/kT), where the migration energy (Em) is obtained from the ab-initio calculation and generalized to a-
Fig. 1. Schematic of Li diffusivity migration and stress distribution in Si: (a) one-phase diffusion, (b) two-phase diffusion.
1012
J. Moon et al. / Journal of Power Sources 272 (2014) 1010e1017
P Si configurations. (iii) Calculate a cumulative function Ki ¼ ij¼1 kj for i ¼ 1, …, N at each current position of the Li. (iv) Generate a pseudorandom number u between 0 and 1. (v) Advance the statistical time of each step (Dt) by eln(u)/KN. (vi) Find one event to carry out i by finding the i for which Ki1 < uK < Ki. (vii) Reconfigure the system according to the chosen event. (viii) Update new position of the Li and simulation time. (ix) Return to the step (i). The self-diffusion coefficient of Li is
D ¼ lim
t/∞
〈r 2 〉 6t
(1)
where t is the time, calculated as the sum of all Dt of each jump and is the mean-squared displacement. 3. Results of DFT and KMC calculations In the following, we analyze the c-Si/a-LixSi two-phase separation under the insertion of Li into c-Si according to three possible mechanisms: thermodynamics of phase separation, diffusion kinetic and phase transformation kinetics from crystalline to amorphous structures. On the basis of the calculations at initial c-Si, we discuss the initial lithiation mechanisms in c-Si and a-Si. The results will give some insights on the mechanisms of two-phase separation by a fast diffusion kinetics at Li rich states as discussed in the following sections. 3.1. Li insertion in c-Si and a-Si We first compute the formation energies of lithiated c-Si and aSi for the optimized LieSi configurations. Assuming that the insertion of Li into Si is homogenous, we put Li atoms randomly at tetrahedral site of Si and make four additional samples (5 random samples for each x in LixSi) to reduce the inhomogeneity and sampling error. For the amorphous structures, both ab initio molecular dynamics simulation and optimization with relaxation of cell volume, shape and atomic position are performed. We have performed separate cell volume only optimization, and subsequent full optimization of the crystalline structure to compare phase transformation kinetics as shown in the next section. The formation energy is defined as:
Ef ðxÞ ¼ ELixSi ðxELi þ ESi Þ
(2)
where E is the energy of each species and x is the molar ratio of Li and Si atoms in the supercell. Comparison of the Li insertion in c-Si and a-Si is shown in Fig 2. The formation energies of c-LixSi (0 < x < 0.4) are positive. With increasing x, the associated formation energy increases and subsequently decreases monotonically from x ¼ 0.3. At x ¼ 0.4 where the curves of formation energy of c-LixSi meets (within error bars) that of a-LixSi, amorphization of c-LixSi spontaneously happens. We note that previous modeling works have presented similar data on the formation energy of Li insertion in c-Si systems, and our DFT results are consistent with them. [27,36,39,40] The blue solid line (coexistence of c-Si and c-LiSi) and green solid line (coexistence of a-Si and a-LiSi) drawn from x ¼ 0 to x ¼ 1 in Fig. 2 indicate that there are different thermodynamic forces for the phase separation between crystalline and amorphous Si structures. The initial incorporation of Li in c-Si is unfavorable (shown by blue solid line below the calculated formation energy values) leading to phase separation of c-Si and c-LiSi (or a-LiSi). However, for the lithiation of a-Si, amorphous LieSi phase would grows stably leading to homogenous phase evolution of LixSi (shown by green solid line above the calculated formation energy values).
Fig. 2. (a) Formation energies as function of Li concentration in c-Si and a-Si. The colored solid line are drawn for comparison of thermodynamic stability. (b) voltagecomposition (Vx) curve for lithiated c-Si and a-Si. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
From the thermodynamic point of view, the phase separation would be expected for incorporation of Li in c-Si by the profile of formation energy. However, a phase separation also requires fast enough kinetics to happen. When c-Si is lithiated under slow charge rate, a phase separation of c-Si and LiSi would not necessarily happen if it is still faster than the phase transformation kinetics (rates of Li diffusion and SieSi bond breaking). To develop a quantitative understanding on the phase transformation kinetics, we examine the mechanisms and kinetics of Li diffusion and SieSi bond breaking in the following sections. As we discuss in the next section, SieSi bond breaking is facilitated by increasing number of Li atoms near the bond, and the bond breaking becomes favorable with 3 or 4 Li atoms (Fig. 3). However, the initial lithiation in c-Si (Figs. 4 and 5) shows repulsive interaction (~0.2 eV) among the Li atoms leading to uniform distribution rather than Li clustering around a SieSi bond. The energy plot in Fig. 2a shows increasingly unfavorable Li formation energy for x < 0.3 which is also indicated by negative voltage
J. Moon et al. / Journal of Power Sources 272 (2014) 1010e1017
1013
should experience complex phase transformations involving SieSi bond breaking. Since it is difficult to define the overall mechanism and the corresponding model systems, we first focus on single SieSi bonding breaking to gain an elementary understanding. We compare the formation energy of fully relaxed c-Si case (shown in Fig. 2) with that of c-Si allowed expanding with fixed shape (i.e., maintaining, c-Si bonding configuration) on the given LieSi configurations. The difference between two energy profiles in Fig. 3a indicates that the crystalline structure of silicon is increasing less stable than the fully relaxed amorphous configuration (energy difference shown as arrows in Fig. 3a) as increasing Li concentration. Consequently, Li insertion is increasing the thermodynamic driving force of phase transformation in Si system. To explain the detailed mechanisms of phase transition, we examine the SieSi bonding elongation with changing Li configurations near a selected SieSi bond. We put Li atoms at Td sites around a particular SieSi bond and impose a gradual elongation of the SieSi bond as shown in Fig 3b, c. Two Si atoms are fixed and the others are fully relaxed to provide the optimized system energy at the given SieSi bond length. Fig. 3d shows the energy profiles as a function of SieSi bonding elongation. Activation energy of SieSi bond breaking is highly dependent on the nearby Li concentration. The incorporation of Li atoms around the Si bond decreases the activation energy of SieSi bond breaking. The activation energy is about 0.1e0.2 eV in lithium rich state and the final state becomes energetically more stable. This finding provides the underlying mechanism of the increasing driving force for phase transformation as the lithium concentration increases (shown in Fig. 3a). The formation energy of c-Si optimized with shape constraint further decreases during the full relaxation. The difference becomes larger as the Li concentration increases indicating that Li atoms become more likely to break SieSi bonds with correspondingly lower energy as shown Fig. 3d. Therefore, the crystalline Si structure transforms to amorphous phase for x > 0.3 with increasingly large thermodynamic force and faster kinetics (0.1e0.2 eV barrier). 3.3. Li kinetic behavior in bulk c-Si To investigate the relationship between Li kinetics and phase transformation, we first examine the Li kinetic behavior in c-Si. We simulate Li in a Si 2 2 2 supercell which contains 64 Si atoms. One or two Li atoms are inserted in the supercell representing the initial charging state. We have performed the DFT simulations of Li migration as a function of the supercell volume change and the distance between two Li atoms.
Fig. 3. (a) Formation energies as function of Li concentration in fully relaxed c-Si and c-Si relaxed with fixed shape constraint (b) The cubic cell of c-Si is outlined by four equivalent tetrahedron sites (labeled 1, 2, 3, 4) surrounding SieSi bond. (c) The configuration of SieSi bonding elongation for [Li 1, 2, 3, 4] case, in units of 2.339 Å. (d) The total energy as a function of the SieSi bond elongation for different numbers of Li atoms surrounding SieSi bond.
in Fig. 2b. As the Li concentration increases (x > 0.3), the Li formation energy decreases and increasing number of Li ions will be present around any Si atom and SieSi bonding configuration leading to bond breaking as analyzed in the next section. 3.2. SieSi bond breaking in LixSi We now investigate the phase transformation kinetic from c-Si to aLixSi with increasing Li concentration. Unlike intercalation electrode materials (e.g., graphite, LiFePO4, LiCoO2), LieSi alloying process
3.3.1. Single Li kinetics In Si diamond structure (a ¼ 5.436 Å in Fig. 4a), we have studied three Li interstitial configurations: tetrahedral (Td), hexagonal (Hex) and bond-center sites. The Td sites are more stable than the Hex and the bond-center sites. When a Li atom is located at a Td site, nearby SieSi bond lengths are increased from 2.339 to 2.376 Å (shown in Fig 4b) as the anti-bonding states are partially occupied due to charge transfer from Li to Si. The Li interstitial defect has only one kind of migration pathway through a Td / Hex / Td trajectory as shown in Fig. 4c. The calculated result shows the migration energy barrier of 0.602 eV through hexagonal ring in [111] direction plotted in Fig. 4d. This result agrees with other theoretical and experimental results. [41e43] Microscopic diffusion constant can be calculated by the Arrhenius equation.
D ¼ 〈r 2 〉=6*f n exp½ Em =ðkTÞ
(3)
The predicted value of Li diffusion constant is about 1012e1013 cm2 s1 at room temperature. To examine the
1014
J. Moon et al. / Journal of Power Sources 272 (2014) 1010e1017
Fig. 4. (a) The unit cell of Si bulk diamond structure and (b) Si bonding distances (in Å) are labeled for both (bulk silicon) and Li interstitial case. The optimized distances between the Li interstitial and its 1st and 2nd neighboring Si atoms are also labeled in (b). (c) The Li interstitial atom in tetrahedral pore of silicon and its neighboring one. (d) The Li atom transition to the neighboring tetrahedron sites with migration energy barrier of 0.6 eV.
quantitative electronic properties, we have also calculated the amount of electron transfers from the inserted Li atom to the Si matrix by the grid-based Bader charge analysis. The amount of charge transfer is 0.83e, and the transferred charges partially fill the anti-bonding sp3 orbital of Si which makes SieSi bonds weaker. This analysis explains
Fig. 5. Li diffusion energy barriers as a function of the distance between two Li atoms. One Li is migrating with the other interstitial Li atom fixed.
the mechanism of progressive weakening of SieSi bonds with increasing Li concentration as discussed in the previous section.
3.3.2. Concentration effects on Li kinetics As the Li concentration increases, Li ions (with positive charge of 0.83e) would interact with each other through screened Coulomb repulsion. Such repulsion is quantitatively described in the total energy increase as an inverse function of LieLi distance. This repulsive interaction subsequently leads to the reduction of the Li migration energy barriers from 0.6 eV to 0.4 eV as illustrated in Fig. 5. The range of repulsive interaction between Li atoms is about 4 Å, and the interaction strength diminishes with increasing distance. The observed repulsion between Li atoms suggests that Li atoms are likely to spread evenly in the Si anode at low Li concentration.
3.3.3. Volume change effects on Li kinetics We next study the volume change effects on Li migration to understand how Li kinetics changes while Si undergoes a large volume expansion during lithiation. In the unstrained c-Si, the Li migration energy barrier is 0.6 eV, and Li atom migration strongly depends on volume change as follows. With three types of strain condition: hydrostatic (εx ¼ εy ¼ εz s 0), biaxial (εx ¼ εy s 0, εz ¼ 0), and uniaxial strain (εx s 0, εy ¼ εz ¼ 0), Li migration energy barriers are calculated. It is found that the migration barrier is determined by the volume change rather than the details of specific strain states. Fig. 6 shows an overlapping linear correlation between migration barrier energy and the volume change for all three strain conditions. If Si is compressed, the Li migration energy barrier increases leading to slower diffusion. In contrast, Li diffusion would be faster when Si is expanded. From these results, we may expect
J. Moon et al. / Journal of Power Sources 272 (2014) 1010e1017
Fig. 6. The schematic of applied strain to c-Si (top insets) and the calculated migration energy barrier as a function of volume change depending on the stretching and compressing deformation.
that the Li diffusivity would be larger in LieSi system which has large volume than in bulk-Si. 3.4. Li kinetic behavior in amorphous silicon The Arrhenius equation (Eqn. (3)) is difficult to use for amorphous Si due to the random nature of local bonding configurations and the corresponding variation of Li migration barriers. However, it is important to note that bulk a-Si maintains the local bonding structure of Si as deformed tetrahedral configuration with the correspondingly deformed interstitial space shapes. From this observation, the Li migration dependence on the interstitial space volume (Fig. 6) enables us to develop an accurate formula to describe Li migration energy as a function of the interstitial volume in the amorphous silicon. Using the volume effects on Li kinetics, we can apply KMC to predict Li diffusion constant in a-Si. To predict Li migration energy barriers in arbitrary amorphous system, we develop the following equation based on local volume change
ln E ri ; rj ¼ a1 $Vi þ a2 $Vj þ a3 $dij þ a4 $Vunit
(4)
where E is the Li migration energy from ri to rj site; Vi is the local volume of Si tetrahedral structure shown in Fig. 4c; d is the migration distance of Li and Vunit is the total volume of amorphous system. Even though there are conflicting opinions on how to describe the kinetics in Si random networks, [44,45] our massive KMC study in a-Si would provide an important insight on Li kinetics and a quantitative determination of the Li diffusion coefficient in aSi. Model systems of a-Si are generated by using ab initio MD at 1200 K over 3 ps and then fully relaxing the AIMD structures at zero temperature. We have also examined the amorphous structure using radial distribution function as shown in Fig. 7. Since it is not feasible to calculate macroscopic diffusion coefficient in this amorphous system using Arrhenius equation, we have performed KMC simulations using the environment dependent migration
1015
barriers (Eqn. (4)). For the confirmation of KMC simulations, we have calculated the Li diffusion coefficient in c-Si and obtained the value, DLi~1013 cm2 s1 at room temperature, which agrees well with the value from the Arrhenius equation (Eqn. (3)) as compared in Fig. 8. To provide input values for KMC simulation of Li diffusion in a-Si, we have calculated migration energy barriers using Eqn. (4) with respect to local volume of Li interstitial sites and migration distance from one Li site to the other site. The coefficients (a1, a2, a3 and a4) in Eqn. (4) are 4.35, 6.30, 0.87, 0.16, respectively. The average value of DLi with 5 sample KMC simulations in 10 10 10 super cell of a-Si is about 1011 cm2 s1. From these calculations, we would expect that the diffusivity of Li is at least ten times higher in a-Si than in c-Si (as shown in Fig. 8). To provide a fundamental understanding on the fast diffusion mechanism in a-Si, we have examined the local volume of Li interstitial sites as shown in Fig. 9. The local volume of an ideal c-Si tetrahedral site is 2.19 Å3. Compared to the local volume in c-Si, the probability of larger volume tetrahedral sites in a-Si is higher than that of smaller volume tetrahedral sites in a-Si. Overall, the migration energy barriers in randomly generated pathways would be quite smaller than 0.6 eV of c-Si due to larger local volumes. The averaged migration energy barriers from KMC simulations are 0.445 eV for a-Si and 0.576 eV (comparable to DFT results 0.60 eV) for c-Si. Although there are relatively few pathways with a high energy barrier in amorphous networks, they are hindered by an unfavorable thermodynamic of Li insertion, and the diffusion is dominated by percolated low energy pathways. Through Sections 3.1e3.4, we have studied the thermodynamics of LixSi systems over 0 < x < 1 and the kinetics of Li atom in c-Si and a-Si, and the phase separation behaviors can be summarized as follows. In the a-Si case, one can expect that the insertion of Li in aSi would not induce phase separation due to both low formation energy (thermodynamic stability in Fig. 2) and fast diffusion kinetics of Li (Fig. 8). Phase transformation kinetics (bond breaking) would become faster (0.1e0.2 eV barrier) than diffusion kinetics (0.44 eV for a-Si and 0.6 eV for c-Si) as Li concentration increases. On the other hands, the phase separation of c-Si and LiSi is expected to occur during the initial lithiation of c-Si system driven by thermodynamic instability (Fig. 2), slow Li kinetics in c-Si and fast SieSi bond breaking kinetics under Li accumulation. The delithiation of a-LixSi is unlikely to form c-Si as x approaches zero, and the a-Si remains in subsequent cycles. 4. Conclusion We have elucidated the atomic scale mechanisms of the phase separation between c-Si and a-LixSi during lithiation using ab-initio and KMC calculations. The formation energies of LixSi systems have shown the presence of two-phase regions for lithiated c-Si and single phase lithiation of a-Si. Thermodynamic and kinetic properties of c-Si and LixSi are expected to separate Si and LixSi phases and induce phase transition at the interface due to low activation energy. The DFT and KMC studies of Li kinetic in LixSi (0 < x < 1) have shown that Li atom diffusion is affected by the local geometry and concentration of Li, and that the Li diffusivity is higher in a-Si than in c-Si since the bonding network of Si in amorphous structure produces percolating pathways with relatively low energy barriers. This finding provides an insight that the slow diffusivity of Li in c-Si is the primary, kinetic bottleneck and introduces thermodynamic force to drive the phase separation. The large volume difference at the phase boundaries and the corresponding interfacial strain and stress are responsible for mechanical failure of Si anode during the initial lithiation. However, the fast Li diffusion and mixing phase stability in a-Si system would allow continuous Li concentration
1016
J. Moon et al. / Journal of Power Sources 272 (2014) 1010e1017
Fig. 7. Comparison of the amorphous and diamond cubic crystals of Si based on structure figures (a, b) and radial distribution functions (c).
Fig. 8. Temperature dependence of the lithium diffusion coefficient in crystalline and amorphous silicon using Arrhenius equation and kinetic Monte Carlo simulations.
Fig. 9. Statistical distribution of a local volume of Si tetrahedral in crystalline silicon (red dot line) and amorphous silicon (blue bar plot). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
J. Moon et al. / Journal of Power Sources 272 (2014) 1010e1017
without forming any sharp boundary. Eventually, an interfacial fracture might develop at the boundary within c-Si leading to rapid capacity loss. [21,22,46e48] Our DFT and KMC studies provide quantitative and mechanistic understanding on the mechanical failure of the Si anode, and these findings would provide promising future research directions to improve the silicon anode performance suitable for lithium ion battery applications. Acknowledgments This work was supported by the Industrial Strategic Technology Development Program (10041589, Development of Web based Multi-scale Simulation Platform for the Efficient Design of Energy Nano Materials) funded by the Ministry of Knowledge Economy (MKE, Korea), and by the National Research Foundation of Korea (NRF) grant funded by the Korea Government (MSIP) (No. 2012R1A3A2048841). References [1] J.M. Tarascon, M. Armand, Nature 414 (2001) 359e367. [2] P. Limthongkul, Y.I. Jang, N.J. Dudney, Y.M. Chiang, J. Power Sources 119 (2003) 604e609. [3] T.D. Hatchard, J.R. Dahn, J. Electrochem. Soc. 151 (2004) A838eA842. [4] M.N. Obrovac, L. Christensen, Electrochem. Solid St. 7 (2004) A93eA96. [5] C.K. Chan, H.L. Peng, G. Liu, K. McIlwrath, X.F. Zhang, R.A. Huggins, Y. Cui, Nat. Nanotechnol. 3 (2008) 31e35. [6] V.L. Chevrier, J.R. Dahn, J. Electrochem. Soc. 156 (2009) A454eA458. [7] V.B. Shenoy, P. Johari, Y. Qi, J. Power Sources 195 (2010) 6825e6830. [8] J. Moon, K. Cho, M. Cho, Int. J. Precis. Eng. Man. 13 (2012) 1191e1197. [9] B. Key, R. Bhattacharyya, M. Morcrette, V. Seznec, J.-M. Tarascon, C.P. Grey, J. Am. Chem. Soc. 131 (2009) 9239e9249. [10] V.L. Chevrier, J.W. Zwanziger, J.R. Dahn, J. Alloy Compd. 496 (2010) 25e36. [11] M.J. Chon, V.A. Sethuraman, A. McCormick, V. Srinivasan, P.R. Guduru, Phys. Rev. Lett. 107 (2011) 045503. [12] V.A. Sethuraman, V. Srinivasan, A.F. Bower, P.R. Guduru, J. Electrochem. Soc. 157 (2010) A1253eA1261. [13] J. Christensen, J. Newman, J. Solid State Electr. 10 (2006) 293e319. [14] S. Golmon, K. Maute, S.H. Lee, M.L. Dunn, Appl. Phys. Lett. 97 (2010). [15] A.F. Bower, P.R. Guduru, V.A. Sethuraman, J. Mech. Phys. Solids 59 (2011) 804e828. [16] Z.W. Cui, F. Gao, J.M. Qu, J. Mech. Phys. Solids 60 (2012) 1280e1295.
́
1017
[17] R. Grantab, V.B. Shenoy, J. Electrochem. Soc. 159 (2012) A584eA591. [18] K. Zhao, W.L. Wang, J. Gregoire, M. Pharr, Z. Suo, J.J. Vlassak, E. Kaxiras, Nano Lett. 11 (2011) 2962e2967. [19] K.J. Zhao, M. Pharr, S.Q. Cai, J.J. Vlassak, Z.G. Suo, J. Am. Ceram. Soc. 94 (2011) S226eS235. [20] K.J. Zhao, M. Pharr, J.J. Vlassak, Z.G. Suo, J. Appl. Phys. 109 (2011). [21] J.W. Wang, Y. He, F.F. Fan, X.H. Liu, S.M. Xia, Y. Liu, C.T. Harris, H. Li, J.Y. Huang, S.X. Mao, T. Zhu, Nano Lett. 13 (2013) 709e715. [22] M.T. McDowell, S.W. Lee, J.T. Harris, B.A. Korgel, C.M. Wang, W.D. Nix, Y. Cui, Nano Lett. 13 (2013) 758e764. [23] C.J. Wen, R.A. Huggins, J. Solid State Chem. 37 (1981) 271e278. [24] K. Yoshimura, J. Suzuki, K. Sekine, T. Takamura, J. Power Sources 174 (2007) 653e657. [25] R. Chandrasekaran, A. Magasinski, G. Yushin, T.F. Fuller, J. Electrochem. Soc. 157 (2010) A1139eA1151. [26] V.L. Chevrier, H.M. Dahn, J.R. Dahn, J. Electrochem. Soc. 158 (2011) A1207eA1213. [27] J. Rohrer, K. Albe, J. Phys. Chem. C 117 (2013) 18796e18803. [28] M. Pharr, K.J. Zhao, X.W. Wang, Z.G. Suo, J.J. Vlassak, Nano Lett. 12 (2012) 5039e5047. [29] Y.-M. Kang, S.-B. Suh, Y.-S. Kim, Inorg. Chem. 48 (2009) 11631e11635. [30] Y. Kubota, M.C.S. Escano, H. Nakanishi, H. Kasai, J. Appl. Phys. 102 (2007). [31] Y.H. Xu, G.P. Yin, P.J. Zuo, Electrochim. Acta 54 (2008) 341e345. [32] G. Kresse, J. Hafner, Phys. Rev. B 48 (1993) 13115e13118. €hnle, Phys. Rev. Lett. 77 (1996) 518e521. [33] W. Frank, U. Breier, C. Els€ asser, M. Fa [34] M.L. Tiago, S. Ismail-Beigi, S.G. Louie, Phys. Rev. B 69 (2004). [35] P.E. Blochl, Phys. Rev. B 50 (1994) 17953e17979. [36] V.L. Chevrier, J.R. Dahn, J. Electrochem. Soc. 157 (2010) A392eA398. [37] M.I. Larsson, R.F. Sabiryanov, K. Cho, B.M. Clemens, J. Appl. Phys. 94 (2003) 3470e3484. [38] P. Ramanarayanan, K.J. Cho, B.M. Clemens, J. Appl. Phys. 94 (2003) 174e185. [39] S.C. Jung, Y.K. Han, Electrochim. Acta 62 (2012) 73e76. [40] H. Kim, C.-Y. Chou, J.G. Ekerdt, G.S. Hwang, J. Phys. Chem. C 115 (2011) 2514e2521. [41] E.M. Pell, Phys. Rev. 119 (1960) 1222e1225. [42] H. Kim, K.E. Kweon, C.-Y. Chou, J.G. Ekerdt, G.S. Hwang, J. Phys. Chem. C 114 (2010) 17942e17946. [43] W. Wan, Q. Zhang, Y. Cui, E. Wang, J. Phys. Condens. Matter 22 (2010) 415501. [44] G.A. Tritsaris, K. Zhao, O.U. Okeke, E. Kaxiras, J. Phys. Chem. C 116 (2012) 22212e22216. [45] A.S. Fedorov, Z.I. Popov, A.A. Kuzubov, S.G. Ovchinnikov, Jetp Lett.þ 95 (2012) 143e147. [46] M.K.Y. Chan, C. Wolverton, J.P. Greeley, J. Am. Chem. Soc. 134 (2012) 14362e14374. [47] S.C. Jung, J.W. Choi, Y.K. Han, Nano Lett. 12 (2012) 5342e5347. [48] X.H. Liu, J.W. Wang, S. Huang, F.F. Fan, X. Huang, Y. Liu, S. Krylyuk, J. Yoo, S.A. Dayeh, A.V. Davydov, S.X. Mao, S.T. Picraux, S.L. Zhang, J. Li, T. Zhu, J.Y. Huang, Nat. Nanotechnol. 7 (2012) 749e756.
Contents lists available at ScienceDirect
Journal of Power Sources journal homepage: www.elsevier.com/locate/jpowsour
Ab initio and kinetic Monte Carlo simulation study of lithiation in crystalline and amorphous silicon Janghyuk Moon a, Byeongchan Lee b, Maenghyo Cho a, Kyeongjae Cho a, c, * a
WCU Multiscale Mechanical Design Division, Department of Mechanical and Aerospace Engineering, Seoul National University, Seoul 151-742, Republic of Korea b Department of Mechanical Engineering, Kyung Hee University, Yongin 446-701, Republic of Korea c Department of Materials Science and Engineering and Department of Physics, The University of Texas at Dallas, Richardson, TX 75080, USA
h i g h l i g h t s DFT and KMC study of energetics and kinetics of c-Si and a-Si under Li insertion. Phase transformation of c-Si to a-Si under Li insertion via SieSi bond breaking. Li diffusion in a-Si predicted by environment dependent KMC method. KMC shows larger Li diffusivity in a-Si than c-Si at room temperature.
a r t i c l e i n f o
a b s t r a c t
Article history: Received 6 February 2014 Received in revised form 19 August 2014 Accepted 2 September 2014 Available online 16 September 2014
Energetics and kinetics of Li insertion into c-Si and a-Si systems are investigated using the density functional (DFT) theory calculations and kinetic Monte Carlo (KMC) simulations. DFT formation energies show the mechanism of phase separation between crystalline silicon and amorphous lithium silicide. Both crystalline and amorphous Si show similar trends in volume expansion and phase transition under lithiation, and kinetics of Li diffusion in bulk silicon (from DFT and KMC) shows a big difference between c-Si and a-Si. The Li migration barrier is 0.6 eV in c-Si, and quickly decreases to 0.4 eV under increasing Li concentration or Si volume expansion. To simulate Li diffusion in amorphous silicon using KMC, we have developed a formulation for environment dependent migration energy barriers of Li in a-Si using a volume dependent function. KMC simulations are performed for Li diffusion in both c-Si and a-Si, and the diffusion coefficient of Li in a-Si is an order of magnitude larger than in c-Si. These studies help to understand mechanisms of lithiation with atomic scale details and elucidate the phase separation between c-Si and lithium silicide. © 2014 Elsevier B.V. All rights reserved.
Keywords: Ab initio method Kinetic Monte Carlo simulation Silicon Li diffusion Phase separation
1. Introduction The lithiation of silicon has attracted a lot of research interest due to the possible application of nano-structured silicon as a high capacity anode material for Li-ion batteries. [1] Compared to graphite (372 mAh g1), silicon has a large theoretical charge capacity (4200 mAh g1) which is more than 10 times higher than that of graphite. However the lithiation and delithiation characteristics of Si are qualitatively different from those of graphitic
* Corresponding author. Department of Materials Science and Engineering and Department of Physics, The University of Texas at Dallas, Richardson, TX 75080, USA. E-mail address: [email protected] (K. Cho). http://dx.doi.org/10.1016/j.jpowsour.2014.09.004 0378-7753/© 2014 Elsevier B.V. All rights reserved.
carbon anode. In contrast to intercalation compounds (graphite anode and transition metal oxide cathode materials), Si anode undergoes phase transformation involving SieSi bond breaking as the Li concentration increases in LixSi. For x 1, diverse crystalline phases (LiSi, Li12Si7, Li7Si3, Li13Si4, Li15Si4, Li22Si5) are known to be thermodynamically stable. However, room temperature lithiation and delithiation studies have shown a crystalline-to-amorphous phase transition that takes place during the first charge cycle of crystalline Si anode. [2e4] Furthermore, the nonconventional phase change characteristics of the lithiation of Si introduces a volume changes up to 400%, an elastic softening, and a phase separation during insertion of Li [5e8]. The lithiation of Si results in substantial atomic structure changes, in contrast to transition metal oxides or graphite electrode materials where Li ion intercalates into existing layers with limited
J. Moon et al. / Journal of Power Sources 272 (2014) 1010e1017
structural distortions. Specifically, Si transforms into various crystalline LixSi phases at high temperature (400 C), but amorphous lithium silicides (LixSi) are formed during lithiation at room temperature. [9e11] The formation of crystalline or amorphous phases during the lithiation of Si is known to be associated with capacity fade. [5,11,12] Previous studies have attempted to explain the performance degradation of Si anode resulting from the complex phase change during lithiation. The strong coupling between diffusion and large mechanical deformation would induce the corresponding mechanical stress, [13e16] crack formation [11,17], and mechanical weakness due to the formation of weak LieSi bonds. [18e20] Even though these previous studies have explained certain aspects of Si anode degradation mechanisms, quantitative nature of atomic structure evolution during the lithiation of Si is not well understood yet. Specifically, the initial lithiation of crystalline silicon is known to produce amorphous LixSi compound, but the detailed nature of amorphous phase formation is not understood yet. Furthermore, how Li diffuses in a-LixSi is not well understood, either. In recent studies, [21,22] many investigations on failure mechanism of Si anode have focused on the issues of phase separation and diffusion kinetics of Li because they significantly affect the stress state at interfaces or lithiated domains during lithiation cycles (Fig. 1). Fig. 1b schematically shows a phase separation between LixSi and Si, and large interfacial stresses are expected to induce mechanical failure modes. Furthermore, the Li diffusion is a key parameter in determining how fast a battery can be cycled (i.e., power density). [2,9,15,23e25] Current experimental and theoretical reports on the failure mechanisms of LieSi systems are not fully consistent with each other. Specifically, several DFT-based studies of LieSi systems have been reported, focusing on Li behavior studies in bulk or nanostructure, atomic and electronic properties of LieSi structure based on high temperature crystalline compounds, or aspects of amorphous model structure generated by abinitio molecular dynamic (AIMD) simulation. Two-phase (Si and LixSi) separation (shown in Fig. 1b) is generally known as the key failure mechanism of Si anode which arises from a thermodynamic driving force, diffusion induced Li accumulation, and a high activation energy of phase transition on the SieLixSi interface. [26] However, a quantitative theoretical basis to explain the phase separation phenomenon is yet to be developed even though many authors have reported different aspect of it. It has been reported in experimental studies that the thermodynamic driving force of two-phase separation and activation energy of phase transition is high, but previous theoretical studies are not consistent with these finding. [27,28] Consequently, it is necessary to develop more detailed theoretical understanding based on DFT and kinetic Monte Carlo (KMC) simulation studies to elucidate the phase separation between c-Si and a-LixSi. Current DFT and KMC studies in this paper show complex phenomena of initial phase instability and Li accumulation due to slow Li diffusion in c-Si as we discuss later.
1011
Kang et al. have experimentally shown that crystalline LiSi phases are present during battery operation cycles, although primarily amorphous LixSi states are present for most Li composition during charging and discharging. [29] The maximum capacity of Si anode at room temperature has been attributed to the crystalline Li15Si4 phase. [30,31] In this work, we focus on initial stage of Li insertion into c-Si before the occurrence of crystalline LixSi (i.e., x < 1 in LixSi) and investigate the energetics and kinetics of initial lithiation. 2. Methodology 2.1. Computational methods The ab initio calculations are performed using Vienna ab initio simulation package (VASP) for the density functional theory (DFT) study of LixSi model systems within the local density approximation. [32e34] The projector augmented wave method is chosen to describe the interaction between ion core and valence electrons. [35] For Li atom, the Li-sv pseudo-potential treated the semi-core 1s states as valence electrons for accurate description of Liþ ionic states. A kepoint mesh in the Monkhorst-Pack scheme is set by 3 3 3 for 2 2 2 supercell of Si. The plane wave basis cut off energy is 500 eV. The structures are fully relaxed for both internal atomic coordinates and supercell shapes. The geometry optimizations are carried out until the forces on each atom become smaller than 0.01 eV Å1. For ab initio molecular dynamics simulation, the kepoint is sampled only at the gamma point. To build amorphous structures, the MD simulations are performed at 1200 K to randomize their initial crystalline configurations. A time step of 3.0 fs and overall simulation steps of 2000 are used. Insertion of Li atoms is done in such a way to achieve a homogenous distribution. At least five samples are generated at each Li concentration (x) in LixSi models. The Delaunary triangulation method (previously used to generate amorphous LieSi bulk) is chosen for Li insertion. [36] To study Li kinetics in Si environment, diffusion pathways and barriers were determined using the climbing-image nudged elastic band (NEB) method. Macroscopic diffusion coefficients of Li in c-Si and aSi are obtained using Arrhenius equation and kinetic Monte Carlo (KMC) simulations. The detailed KMC algorithm for simulating the time evolution follows the conventional procedure [37,38]. Results of the KMC simulations are statistically averaged to obtain macroscopic behavior. In this study, the standard KMC algorithm is used to simulate an activated random-walk process of Li interstitial defects in Si. Movements of the Li defects have followed the general KMC procedures: (i) Identify all possible events from the current atomic configuration. In this study, we consider first-nearest-neighbor (<3.5 Å) Li sites from sites of tetrahedral Si. (ii) Determine a list of the rates (ki) of all possible transition events in the system. The rates are proportional to exp(eEm/kT), where the migration energy (Em) is obtained from the ab-initio calculation and generalized to a-
Fig. 1. Schematic of Li diffusivity migration and stress distribution in Si: (a) one-phase diffusion, (b) two-phase diffusion.
1012
J. Moon et al. / Journal of Power Sources 272 (2014) 1010e1017
P Si configurations. (iii) Calculate a cumulative function Ki ¼ ij¼1 kj for i ¼ 1, …, N at each current position of the Li. (iv) Generate a pseudorandom number u between 0 and 1. (v) Advance the statistical time of each step (Dt) by eln(u)/KN. (vi) Find one event to carry out i by finding the i for which Ki1 < uK < Ki. (vii) Reconfigure the system according to the chosen event. (viii) Update new position of the Li and simulation time. (ix) Return to the step (i). The self-diffusion coefficient of Li is
D ¼ lim
t/∞
〈r 2 〉 6t
(1)
where t is the time, calculated as the sum of all Dt of each jump and
Ef ðxÞ ¼ ELixSi ðxELi þ ESi Þ
(2)
where E is the energy of each species and x is the molar ratio of Li and Si atoms in the supercell. Comparison of the Li insertion in c-Si and a-Si is shown in Fig 2. The formation energies of c-LixSi (0 < x < 0.4) are positive. With increasing x, the associated formation energy increases and subsequently decreases monotonically from x ¼ 0.3. At x ¼ 0.4 where the curves of formation energy of c-LixSi meets (within error bars) that of a-LixSi, amorphization of c-LixSi spontaneously happens. We note that previous modeling works have presented similar data on the formation energy of Li insertion in c-Si systems, and our DFT results are consistent with them. [27,36,39,40] The blue solid line (coexistence of c-Si and c-LiSi) and green solid line (coexistence of a-Si and a-LiSi) drawn from x ¼ 0 to x ¼ 1 in Fig. 2 indicate that there are different thermodynamic forces for the phase separation between crystalline and amorphous Si structures. The initial incorporation of Li in c-Si is unfavorable (shown by blue solid line below the calculated formation energy values) leading to phase separation of c-Si and c-LiSi (or a-LiSi). However, for the lithiation of a-Si, amorphous LieSi phase would grows stably leading to homogenous phase evolution of LixSi (shown by green solid line above the calculated formation energy values).
Fig. 2. (a) Formation energies as function of Li concentration in c-Si and a-Si. The colored solid line are drawn for comparison of thermodynamic stability. (b) voltagecomposition (Vx) curve for lithiated c-Si and a-Si. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
From the thermodynamic point of view, the phase separation would be expected for incorporation of Li in c-Si by the profile of formation energy. However, a phase separation also requires fast enough kinetics to happen. When c-Si is lithiated under slow charge rate, a phase separation of c-Si and LiSi would not necessarily happen if it is still faster than the phase transformation kinetics (rates of Li diffusion and SieSi bond breaking). To develop a quantitative understanding on the phase transformation kinetics, we examine the mechanisms and kinetics of Li diffusion and SieSi bond breaking in the following sections. As we discuss in the next section, SieSi bond breaking is facilitated by increasing number of Li atoms near the bond, and the bond breaking becomes favorable with 3 or 4 Li atoms (Fig. 3). However, the initial lithiation in c-Si (Figs. 4 and 5) shows repulsive interaction (~0.2 eV) among the Li atoms leading to uniform distribution rather than Li clustering around a SieSi bond. The energy plot in Fig. 2a shows increasingly unfavorable Li formation energy for x < 0.3 which is also indicated by negative voltage
J. Moon et al. / Journal of Power Sources 272 (2014) 1010e1017
1013
should experience complex phase transformations involving SieSi bond breaking. Since it is difficult to define the overall mechanism and the corresponding model systems, we first focus on single SieSi bonding breaking to gain an elementary understanding. We compare the formation energy of fully relaxed c-Si case (shown in Fig. 2) with that of c-Si allowed expanding with fixed shape (i.e., maintaining, c-Si bonding configuration) on the given LieSi configurations. The difference between two energy profiles in Fig. 3a indicates that the crystalline structure of silicon is increasing less stable than the fully relaxed amorphous configuration (energy difference shown as arrows in Fig. 3a) as increasing Li concentration. Consequently, Li insertion is increasing the thermodynamic driving force of phase transformation in Si system. To explain the detailed mechanisms of phase transition, we examine the SieSi bonding elongation with changing Li configurations near a selected SieSi bond. We put Li atoms at Td sites around a particular SieSi bond and impose a gradual elongation of the SieSi bond as shown in Fig 3b, c. Two Si atoms are fixed and the others are fully relaxed to provide the optimized system energy at the given SieSi bond length. Fig. 3d shows the energy profiles as a function of SieSi bonding elongation. Activation energy of SieSi bond breaking is highly dependent on the nearby Li concentration. The incorporation of Li atoms around the Si bond decreases the activation energy of SieSi bond breaking. The activation energy is about 0.1e0.2 eV in lithium rich state and the final state becomes energetically more stable. This finding provides the underlying mechanism of the increasing driving force for phase transformation as the lithium concentration increases (shown in Fig. 3a). The formation energy of c-Si optimized with shape constraint further decreases during the full relaxation. The difference becomes larger as the Li concentration increases indicating that Li atoms become more likely to break SieSi bonds with correspondingly lower energy as shown Fig. 3d. Therefore, the crystalline Si structure transforms to amorphous phase for x > 0.3 with increasingly large thermodynamic force and faster kinetics (0.1e0.2 eV barrier). 3.3. Li kinetic behavior in bulk c-Si To investigate the relationship between Li kinetics and phase transformation, we first examine the Li kinetic behavior in c-Si. We simulate Li in a Si 2 2 2 supercell which contains 64 Si atoms. One or two Li atoms are inserted in the supercell representing the initial charging state. We have performed the DFT simulations of Li migration as a function of the supercell volume change and the distance between two Li atoms.
Fig. 3. (a) Formation energies as function of Li concentration in fully relaxed c-Si and c-Si relaxed with fixed shape constraint (b) The cubic cell of c-Si is outlined by four equivalent tetrahedron sites (labeled 1, 2, 3, 4) surrounding SieSi bond. (c) The configuration of SieSi bonding elongation for [Li 1, 2, 3, 4] case, in units of 2.339 Å. (d) The total energy as a function of the SieSi bond elongation for different numbers of Li atoms surrounding SieSi bond.
in Fig. 2b. As the Li concentration increases (x > 0.3), the Li formation energy decreases and increasing number of Li ions will be present around any Si atom and SieSi bonding configuration leading to bond breaking as analyzed in the next section. 3.2. SieSi bond breaking in LixSi We now investigate the phase transformation kinetic from c-Si to aLixSi with increasing Li concentration. Unlike intercalation electrode materials (e.g., graphite, LiFePO4, LiCoO2), LieSi alloying process
3.3.1. Single Li kinetics In Si diamond structure (a ¼ 5.436 Å in Fig. 4a), we have studied three Li interstitial configurations: tetrahedral (Td), hexagonal (Hex) and bond-center sites. The Td sites are more stable than the Hex and the bond-center sites. When a Li atom is located at a Td site, nearby SieSi bond lengths are increased from 2.339 to 2.376 Å (shown in Fig 4b) as the anti-bonding states are partially occupied due to charge transfer from Li to Si. The Li interstitial defect has only one kind of migration pathway through a Td / Hex / Td trajectory as shown in Fig. 4c. The calculated result shows the migration energy barrier of 0.602 eV through hexagonal ring in [111] direction plotted in Fig. 4d. This result agrees with other theoretical and experimental results. [41e43] Microscopic diffusion constant can be calculated by the Arrhenius equation.
D ¼ 〈r 2 〉=6*f n exp½ Em =ðkTÞ
(3)
The predicted value of Li diffusion constant is about 1012e1013 cm2 s1 at room temperature. To examine the
1014
J. Moon et al. / Journal of Power Sources 272 (2014) 1010e1017
Fig. 4. (a) The unit cell of Si bulk diamond structure and (b) Si bonding distances (in Å) are labeled for both (bulk silicon) and Li interstitial case. The optimized distances between the Li interstitial and its 1st and 2nd neighboring Si atoms are also labeled in (b). (c) The Li interstitial atom in tetrahedral pore of silicon and its neighboring one. (d) The Li atom transition to the neighboring tetrahedron sites with migration energy barrier of 0.6 eV.
quantitative electronic properties, we have also calculated the amount of electron transfers from the inserted Li atom to the Si matrix by the grid-based Bader charge analysis. The amount of charge transfer is 0.83e, and the transferred charges partially fill the anti-bonding sp3 orbital of Si which makes SieSi bonds weaker. This analysis explains
Fig. 5. Li diffusion energy barriers as a function of the distance between two Li atoms. One Li is migrating with the other interstitial Li atom fixed.
the mechanism of progressive weakening of SieSi bonds with increasing Li concentration as discussed in the previous section.
3.3.2. Concentration effects on Li kinetics As the Li concentration increases, Li ions (with positive charge of 0.83e) would interact with each other through screened Coulomb repulsion. Such repulsion is quantitatively described in the total energy increase as an inverse function of LieLi distance. This repulsive interaction subsequently leads to the reduction of the Li migration energy barriers from 0.6 eV to 0.4 eV as illustrated in Fig. 5. The range of repulsive interaction between Li atoms is about 4 Å, and the interaction strength diminishes with increasing distance. The observed repulsion between Li atoms suggests that Li atoms are likely to spread evenly in the Si anode at low Li concentration.
3.3.3. Volume change effects on Li kinetics We next study the volume change effects on Li migration to understand how Li kinetics changes while Si undergoes a large volume expansion during lithiation. In the unstrained c-Si, the Li migration energy barrier is 0.6 eV, and Li atom migration strongly depends on volume change as follows. With three types of strain condition: hydrostatic (εx ¼ εy ¼ εz s 0), biaxial (εx ¼ εy s 0, εz ¼ 0), and uniaxial strain (εx s 0, εy ¼ εz ¼ 0), Li migration energy barriers are calculated. It is found that the migration barrier is determined by the volume change rather than the details of specific strain states. Fig. 6 shows an overlapping linear correlation between migration barrier energy and the volume change for all three strain conditions. If Si is compressed, the Li migration energy barrier increases leading to slower diffusion. In contrast, Li diffusion would be faster when Si is expanded. From these results, we may expect
J. Moon et al. / Journal of Power Sources 272 (2014) 1010e1017
Fig. 6. The schematic of applied strain to c-Si (top insets) and the calculated migration energy barrier as a function of volume change depending on the stretching and compressing deformation.
that the Li diffusivity would be larger in LieSi system which has large volume than in bulk-Si. 3.4. Li kinetic behavior in amorphous silicon The Arrhenius equation (Eqn. (3)) is difficult to use for amorphous Si due to the random nature of local bonding configurations and the corresponding variation of Li migration barriers. However, it is important to note that bulk a-Si maintains the local bonding structure of Si as deformed tetrahedral configuration with the correspondingly deformed interstitial space shapes. From this observation, the Li migration dependence on the interstitial space volume (Fig. 6) enables us to develop an accurate formula to describe Li migration energy as a function of the interstitial volume in the amorphous silicon. Using the volume effects on Li kinetics, we can apply KMC to predict Li diffusion constant in a-Si. To predict Li migration energy barriers in arbitrary amorphous system, we develop the following equation based on local volume change
ln E ri ; rj ¼ a1 $Vi þ a2 $Vj þ a3 $dij þ a4 $Vunit
(4)
where E is the Li migration energy from ri to rj site; Vi is the local volume of Si tetrahedral structure shown in Fig. 4c; d is the migration distance of Li and Vunit is the total volume of amorphous system. Even though there are conflicting opinions on how to describe the kinetics in Si random networks, [44,45] our massive KMC study in a-Si would provide an important insight on Li kinetics and a quantitative determination of the Li diffusion coefficient in aSi. Model systems of a-Si are generated by using ab initio MD at 1200 K over 3 ps and then fully relaxing the AIMD structures at zero temperature. We have also examined the amorphous structure using radial distribution function as shown in Fig. 7. Since it is not feasible to calculate macroscopic diffusion coefficient in this amorphous system using Arrhenius equation, we have performed KMC simulations using the environment dependent migration
1015
barriers (Eqn. (4)). For the confirmation of KMC simulations, we have calculated the Li diffusion coefficient in c-Si and obtained the value, DLi~1013 cm2 s1 at room temperature, which agrees well with the value from the Arrhenius equation (Eqn. (3)) as compared in Fig. 8. To provide input values for KMC simulation of Li diffusion in a-Si, we have calculated migration energy barriers using Eqn. (4) with respect to local volume of Li interstitial sites and migration distance from one Li site to the other site. The coefficients (a1, a2, a3 and a4) in Eqn. (4) are 4.35, 6.30, 0.87, 0.16, respectively. The average value of DLi with 5 sample KMC simulations in 10 10 10 super cell of a-Si is about 1011 cm2 s1. From these calculations, we would expect that the diffusivity of Li is at least ten times higher in a-Si than in c-Si (as shown in Fig. 8). To provide a fundamental understanding on the fast diffusion mechanism in a-Si, we have examined the local volume of Li interstitial sites as shown in Fig. 9. The local volume of an ideal c-Si tetrahedral site is 2.19 Å3. Compared to the local volume in c-Si, the probability of larger volume tetrahedral sites in a-Si is higher than that of smaller volume tetrahedral sites in a-Si. Overall, the migration energy barriers in randomly generated pathways would be quite smaller than 0.6 eV of c-Si due to larger local volumes. The averaged migration energy barriers from KMC simulations are 0.445 eV for a-Si and 0.576 eV (comparable to DFT results 0.60 eV) for c-Si. Although there are relatively few pathways with a high energy barrier in amorphous networks, they are hindered by an unfavorable thermodynamic of Li insertion, and the diffusion is dominated by percolated low energy pathways. Through Sections 3.1e3.4, we have studied the thermodynamics of LixSi systems over 0 < x < 1 and the kinetics of Li atom in c-Si and a-Si, and the phase separation behaviors can be summarized as follows. In the a-Si case, one can expect that the insertion of Li in aSi would not induce phase separation due to both low formation energy (thermodynamic stability in Fig. 2) and fast diffusion kinetics of Li (Fig. 8). Phase transformation kinetics (bond breaking) would become faster (0.1e0.2 eV barrier) than diffusion kinetics (0.44 eV for a-Si and 0.6 eV for c-Si) as Li concentration increases. On the other hands, the phase separation of c-Si and LiSi is expected to occur during the initial lithiation of c-Si system driven by thermodynamic instability (Fig. 2), slow Li kinetics in c-Si and fast SieSi bond breaking kinetics under Li accumulation. The delithiation of a-LixSi is unlikely to form c-Si as x approaches zero, and the a-Si remains in subsequent cycles. 4. Conclusion We have elucidated the atomic scale mechanisms of the phase separation between c-Si and a-LixSi during lithiation using ab-initio and KMC calculations. The formation energies of LixSi systems have shown the presence of two-phase regions for lithiated c-Si and single phase lithiation of a-Si. Thermodynamic and kinetic properties of c-Si and LixSi are expected to separate Si and LixSi phases and induce phase transition at the interface due to low activation energy. The DFT and KMC studies of Li kinetic in LixSi (0 < x < 1) have shown that Li atom diffusion is affected by the local geometry and concentration of Li, and that the Li diffusivity is higher in a-Si than in c-Si since the bonding network of Si in amorphous structure produces percolating pathways with relatively low energy barriers. This finding provides an insight that the slow diffusivity of Li in c-Si is the primary, kinetic bottleneck and introduces thermodynamic force to drive the phase separation. The large volume difference at the phase boundaries and the corresponding interfacial strain and stress are responsible for mechanical failure of Si anode during the initial lithiation. However, the fast Li diffusion and mixing phase stability in a-Si system would allow continuous Li concentration
1016
J. Moon et al. / Journal of Power Sources 272 (2014) 1010e1017
Fig. 7. Comparison of the amorphous and diamond cubic crystals of Si based on structure figures (a, b) and radial distribution functions (c).
Fig. 8. Temperature dependence of the lithium diffusion coefficient in crystalline and amorphous silicon using Arrhenius equation and kinetic Monte Carlo simulations.
Fig. 9. Statistical distribution of a local volume of Si tetrahedral in crystalline silicon (red dot line) and amorphous silicon (blue bar plot). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
J. Moon et al. / Journal of Power Sources 272 (2014) 1010e1017
without forming any sharp boundary. Eventually, an interfacial fracture might develop at the boundary within c-Si leading to rapid capacity loss. [21,22,46e48] Our DFT and KMC studies provide quantitative and mechanistic understanding on the mechanical failure of the Si anode, and these findings would provide promising future research directions to improve the silicon anode performance suitable for lithium ion battery applications. Acknowledgments This work was supported by the Industrial Strategic Technology Development Program (10041589, Development of Web based Multi-scale Simulation Platform for the Efficient Design of Energy Nano Materials) funded by the Ministry of Knowledge Economy (MKE, Korea), and by the National Research Foundation of Korea (NRF) grant funded by the Korea Government (MSIP) (No. 2012R1A3A2048841). References [1] J.M. Tarascon, M. Armand, Nature 414 (2001) 359e367. [2] P. Limthongkul, Y.I. Jang, N.J. Dudney, Y.M. Chiang, J. Power Sources 119 (2003) 604e609. [3] T.D. Hatchard, J.R. Dahn, J. Electrochem. Soc. 151 (2004) A838eA842. [4] M.N. Obrovac, L. Christensen, Electrochem. Solid St. 7 (2004) A93eA96. [5] C.K. Chan, H.L. Peng, G. Liu, K. McIlwrath, X.F. Zhang, R.A. Huggins, Y. Cui, Nat. Nanotechnol. 3 (2008) 31e35. [6] V.L. Chevrier, J.R. Dahn, J. Electrochem. Soc. 156 (2009) A454eA458. [7] V.B. Shenoy, P. Johari, Y. Qi, J. Power Sources 195 (2010) 6825e6830. [8] J. Moon, K. Cho, M. Cho, Int. J. Precis. Eng. Man. 13 (2012) 1191e1197. [9] B. Key, R. Bhattacharyya, M. Morcrette, V. Seznec, J.-M. Tarascon, C.P. Grey, J. Am. Chem. Soc. 131 (2009) 9239e9249. [10] V.L. Chevrier, J.W. Zwanziger, J.R. Dahn, J. Alloy Compd. 496 (2010) 25e36. [11] M.J. Chon, V.A. Sethuraman, A. McCormick, V. Srinivasan, P.R. Guduru, Phys. Rev. Lett. 107 (2011) 045503. [12] V.A. Sethuraman, V. Srinivasan, A.F. Bower, P.R. Guduru, J. Electrochem. Soc. 157 (2010) A1253eA1261. [13] J. Christensen, J. Newman, J. Solid State Electr. 10 (2006) 293e319. [14] S. Golmon, K. Maute, S.H. Lee, M.L. Dunn, Appl. Phys. Lett. 97 (2010). [15] A.F. Bower, P.R. Guduru, V.A. Sethuraman, J. Mech. Phys. Solids 59 (2011) 804e828. [16] Z.W. Cui, F. Gao, J.M. Qu, J. Mech. Phys. Solids 60 (2012) 1280e1295.
́
1017
[17] R. Grantab, V.B. Shenoy, J. Electrochem. Soc. 159 (2012) A584eA591. [18] K. Zhao, W.L. Wang, J. Gregoire, M. Pharr, Z. Suo, J.J. Vlassak, E. Kaxiras, Nano Lett. 11 (2011) 2962e2967. [19] K.J. Zhao, M. Pharr, S.Q. Cai, J.J. Vlassak, Z.G. Suo, J. Am. Ceram. Soc. 94 (2011) S226eS235. [20] K.J. Zhao, M. Pharr, J.J. Vlassak, Z.G. Suo, J. Appl. Phys. 109 (2011). [21] J.W. Wang, Y. He, F.F. Fan, X.H. Liu, S.M. Xia, Y. Liu, C.T. Harris, H. Li, J.Y. Huang, S.X. Mao, T. Zhu, Nano Lett. 13 (2013) 709e715. [22] M.T. McDowell, S.W. Lee, J.T. Harris, B.A. Korgel, C.M. Wang, W.D. Nix, Y. Cui, Nano Lett. 13 (2013) 758e764. [23] C.J. Wen, R.A. Huggins, J. Solid State Chem. 37 (1981) 271e278. [24] K. Yoshimura, J. Suzuki, K. Sekine, T. Takamura, J. Power Sources 174 (2007) 653e657. [25] R. Chandrasekaran, A. Magasinski, G. Yushin, T.F. Fuller, J. Electrochem. Soc. 157 (2010) A1139eA1151. [26] V.L. Chevrier, H.M. Dahn, J.R. Dahn, J. Electrochem. Soc. 158 (2011) A1207eA1213. [27] J. Rohrer, K. Albe, J. Phys. Chem. C 117 (2013) 18796e18803. [28] M. Pharr, K.J. Zhao, X.W. Wang, Z.G. Suo, J.J. Vlassak, Nano Lett. 12 (2012) 5039e5047. [29] Y.-M. Kang, S.-B. Suh, Y.-S. Kim, Inorg. Chem. 48 (2009) 11631e11635. [30] Y. Kubota, M.C.S. Escano, H. Nakanishi, H. Kasai, J. Appl. Phys. 102 (2007). [31] Y.H. Xu, G.P. Yin, P.J. Zuo, Electrochim. Acta 54 (2008) 341e345. [32] G. Kresse, J. Hafner, Phys. Rev. B 48 (1993) 13115e13118. €hnle, Phys. Rev. Lett. 77 (1996) 518e521. [33] W. Frank, U. Breier, C. Els€ asser, M. Fa [34] M.L. Tiago, S. Ismail-Beigi, S.G. Louie, Phys. Rev. B 69 (2004). [35] P.E. Blochl, Phys. Rev. B 50 (1994) 17953e17979. [36] V.L. Chevrier, J.R. Dahn, J. Electrochem. Soc. 157 (2010) A392eA398. [37] M.I. Larsson, R.F. Sabiryanov, K. Cho, B.M. Clemens, J. Appl. Phys. 94 (2003) 3470e3484. [38] P. Ramanarayanan, K.J. Cho, B.M. Clemens, J. Appl. Phys. 94 (2003) 174e185. [39] S.C. Jung, Y.K. Han, Electrochim. Acta 62 (2012) 73e76. [40] H. Kim, C.-Y. Chou, J.G. Ekerdt, G.S. Hwang, J. Phys. Chem. C 115 (2011) 2514e2521. [41] E.M. Pell, Phys. Rev. 119 (1960) 1222e1225. [42] H. Kim, K.E. Kweon, C.-Y. Chou, J.G. Ekerdt, G.S. Hwang, J. Phys. Chem. C 114 (2010) 17942e17946. [43] W. Wan, Q. Zhang, Y. Cui, E. Wang, J. Phys. Condens. Matter 22 (2010) 415501. [44] G.A. Tritsaris, K. Zhao, O.U. Okeke, E. Kaxiras, J. Phys. Chem. C 116 (2012) 22212e22216. [45] A.S. Fedorov, Z.I. Popov, A.A. Kuzubov, S.G. Ovchinnikov, Jetp Lett.þ 95 (2012) 143e147. [46] M.K.Y. Chan, C. Wolverton, J.P. Greeley, J. Am. Chem. Soc. 134 (2012) 14362e14374. [47] S.C. Jung, J.W. Choi, Y.K. Han, Nano Lett. 12 (2012) 5342e5347. [48] X.H. Liu, J.W. Wang, S. Huang, F.F. Fan, X. Huang, Y. Liu, S. Krylyuk, J. Yoo, S.A. Dayeh, A.V. Davydov, S.X. Mao, S.T. Picraux, S.L. Zhang, J. Li, T. Zhu, J.Y. Huang, Nat. Nanotechnol. 7 (2012) 749e756.