Received: 20 Aug 2018
Revised: 11 Sep 2018
Accepted: 11 Nov 2018
Published online: 12 Nov 2018
Shiqian Hu 1, 2, 3 , Zhongwei Zhang 1, 2, 3, Zhongting Wang 4 , Kaiyang Zeng 4 , Yuan Cheng 5, Jie Chen 1, 2, 3 , Email, Gang Zhang 5, Email
1 Center for Phononics and Thermal Energy Science, School of Physics Science and Engineering, and Institute for Advanced Study, Tongji University, Shanghai 200092, China
2 China–EU Joint Lab for Nanophononics, School of Physics Science and Engineering, Tongji University, Shanghai 200092, China
3 Shanghai Key Laboratory of Special Artificial Microstructure Materials and Technology, School of Physics Science and Engineering, Tongji University, Shanghai 200092, China
4 Department of Mechanical Engineering, National University of Singapore, 9 Engineering Drive 1, 117576, Singapore
5 Institute of High Performance Computing, Agency for Science, Technology and Research, 1 Fusionopolis Way, Singapore 138632, Singapore
1. Introduction
In recent years, Lithium ion (Li-ion) rechargeable batteries have served as significant power sources for most portable electronics as well as the electric vehicles because of their high energy density and rate capability. Safety is a key aspect for Li-ion batteries.1 Under ideal operation, the batteries convert the chemical energy to electrical energy with minimal heat. However, in some abusive cases, the chemical energy may directly convert into heat. For example, overcharge can generate huge heat energy within the battery. Although safety control is not determined by one criterion or parameter, thermal conductivity is still one of the most important material parameters that determine the safety of Li-ion batteries,
modules and battery packs.2 Even under normal operation, charging/discharging circuit will heat up the battery via Joule heating, and this heat is dissipated through the electrode materials in the format of thermal conduction. If the thermal conductivity of electrode material is not sufficiently high, the temperature will continue to rise. Once the temperature is over a critical value, exothermic reaction will begin and the battery itself produces heat, in which high rate reaction will cause temperature rising sharply and rapid disassembly may occur.3
Most commercial Li-ion batteries consist of transition metal oxide cathode such as LiCoO , graphitic carbon as anode and an 2 organic solvent electrolyte. In contrast to the high thermal conductivity of graphite (~thousands Wm-1K-1), room temperature thermal conductivity of LiCoO is only about 5.4 Wm-1K-1.4-6 2 Obviously, compared to graphite based anode, low thermal conductivity in LiCoO will be the bottleneck in thermal 2 management of Li-ion batteries. Moreover, during the charging cycle, lithium is moved from the cathode, while the dynamic (Li concentration dependent) thermal conductivity of LiCoO is still not 2 well understood, which is challenging in experiment. Despite playing critical role in the safety of Li-ion battery, fundamental understanding of how delithiation changes thermal conductivity and temperature dependence are still an uncharted territory.
In recent years, atomistic simulation has played a vital role in the materials design of advanced Li-ion batteries,7-12 and has successfully predicted the properties and key performance of battery systems. In particular, the predictions of a variety of properties pertinent to Li-ion batteries have been shown to be accurate, such as the energy density, ion diffusion mobility and the charge/discharge voltage profiles.13-16 However, previously, the major efforts are focused on the electronic, mechanical and morphological properties, while the thermal property of LiCoO upon delithiation is still not 2 clear, despite its crucial role in Li battery safety. Although a fully quantum-mechanical treatment (represented by first-principles calculations) of the system is highly desirable, it can only be applied to the system with a small size due to the large computational load. Molecular dynamics (MD) simulation is another powerful tool to handle many-body problems at the atomic level, and its applications have covered a wide range of material properties, such as mechanical and thermal properties.17-20 It can explore the thermal conductivity without any thermodynamic-limit assumption, and naturally includes full affecting factors in atomistic interactions. In MD simulations, the force on each atom according to the force field controls the dynamic process of the system. Therefore, tremendous efforts have been devoted to developing empirical potential with the sacrifice of exactness.
In this work, using non-equilibrium molecular dynamics simulations (NEMD), we report the Li-ion flowing modulation of the thermal conductivity in LiCoO over the range of 19.8 Wm-1K-1 to 5.2 2 Wm-1K-1 for the delithiation degree changing from zero to 40%. Based on the heat energy transmission coefficient, the delithiation induces the non-propagating heat transport mechanism is obtained. Moreover, the underlying physical mechanism of the weak temperature-dependence of the thermal conductivity in the delithiated LiCoO is uncovered by comparing the spectral energy 2 density between pristine and delithiated LiCoO . A clear dynamic 2 modulation picture of the thermal conductivity of LiCoO by 2 delithiation is obtained in this work.
1Center for Phononics and Thermal Energy Science, School of Physics Science and Engineering, and Institute for Advanced Study, Tongji University, Shanghai 200092, People's Republic of China
2China–EU Joint Lab for Nanophononics, School of Physics Science and Engineering, Tongji University, Shanghai 200092, People’s Republic of China
3Shanghai Key Laboratory of Special Artificial Microstructure Materials and Technology, School of Physics Science and Engineering, Tongji University, Shanghai 200092, People’s Republic of China
4Department of Mechanical Engineering, National University of Singapore, 9 Engineering Drive 1, 117576, Singapore
5Institute of High Performance Computing, Agency for Science, Technology and Research, 1 Fusionopolis Way, Singapore 138632, Singapore
*E-mail: jie@tongji.edu.cn (J.C.); zhangg@ihpc.a-star.edu.sg (G. Z.)
2. Computational Methods
MD simulations in this paper are performed by using LAMMPS package.21 The Buckingham potential combining with a long-range Coulombic component was used to describe the potential energy between pairs of ions in each crystalline solid. The Buckingham potential is denoted as,22

where r is the distance between ions i and j, and A , ρ and C are the ij ij ij potential parameters taken from Ref 22. The time step is set as 0.5 fs. As shown in Fig. 1a, the periodic boundary conditions are used in all directions. The cut-off distance is set as 10 Å for both Coulombic and Buckingham force. Equilibrium MD simulations with the isothermal-isobaric (NPT) ensemble are first performed for 50 ps, during which the whole system has reached thermal equilibration (Fig. 1b) and the stress in the structure has been fully relaxed (Fig. 1c). Then, the NEMD simulations are performed. To establish a temperature gradient, the atoms near two ends (two shadowed regions in Fig. 1a) are coupled with Langevin heat bath22,23 at lowtemperature T , and atoms at the center are attached to high- L temperature heat bath at T . With the periodic boundary condition, H this simulation setup can effectively remove the boundary effect.25,26 The thermal conductivity κ is calculated based on Fourier’s Law,

where and J are, respectively, the temperature gradient and the heat flux. Fig. 1d shows a typical temperature profile of the LiCoO 2 at 300 K. The temperature gradient is obtained by the linear fitting to the local temperature at each side, excluding the temperature jumps closed to the heat baths. In our simulations, after the system reaching the steady state, the cumulative energy ΔE added/subtracted to the heat source/sink region is recorded for 1 ns. By applying the linear fitting to the raw data of the cumulative energy ΔE (shown in Fig. 1e), the energy change per unit time (ΔE/Δt) is obtained, which is used to calculate the heat flux J = ΔE/(2SΔt). Here S is the cross area of LiCoO defined as the width (W) multiplied by the thickness, and 2 the factor of 2 accounts for the bidirectional heat flow. The size of the simulation domain is fixed as 25 nm, 3 nm and 3 nm along x, y and z direction, respectively. The results presented here are averaged over 3 independent simulations with different initial conditions, and the error bar is obtained from the standard deviation of different runs.

Fig. 1 Simulation setup, structure relaxation, temperature gradient and heat flux calculation. (a) Simulation set up, top view and side view of the LiCoO . Periodic boundary conditions are used along the 2 length (L) and the width (W) direction. The green dashed line denotes the imaginary plane which divides the system into “L” and “R” two segments. The size of the simulation domain is L= 32 nm and W= 3 nm. (b) The average temperature in the simulation system with the NPT ensemble. (c) The pressure distribution in three directions with the NPT ensemble. (d) The typical temperature profile of the LiCoO at 300 K. Linear fitting is performed to obtain 2 the temperature gradient. (e)The typical cumulative energy change in the heat source and heat sink region for the LiCoO at 300 K. Linear 2 fitting is performed to obtain the energy change per unit time (ΔE/Δt).
3. Results and Discussion
3.1 Rapid reduction in thermal conductivity upon delithiation
The delithiation process in our simulations is realized by randomly removing certain amount of Li-ions in LiCoO . The delithiation degree is defined as N /N where N denotes the number of removed R P, R Li-ions and N is the total number of Li-ions in pristine LiCoO . In P 2 our simulation, the delithiation degree is varied up to 40%, which is consistent with previous experimental work.4 As shown in Fig. 2, thermal conductivities of the delithiated LiCoO (κ ) are significantly 2 D reduced compared to that of pristine LiCoO (κ ), showing an 2 P exponential decay behavior (solid line in Fig. 2) in the delithiation (charging) process. When the delithiation degree increases to 20%, a sharp decrease in κ is observed. For instance, κ is reduced to 65% D D of κ with only 5% of delithiation, and it is only 30% of κ with 20% P P of delithiation. Further increase of delithiation degree above 20% induces less significant change in κ . This behavior is consistent with D the previous MD simulations of the defected MoS 27 and graphene 2 sheet.28

Fig. 2. Thermal conductivity of LiCoO changes with the 2 delithiation degree at 300 K. The symbols are numerical simulation results and the solid line is the exponential fitting curve. The insert shows 1/κ versus the delithiation degree, and the dashed line draws the linear dependence for comparison.
There are Li-vacancies in the delithiated LiCoO . Phonon scattering 2 by vacancies in crystals is usually treated from perturbation theory in terms of the missing mass and the missing linkages, and the scattering rate is29

where x is the density of vacancies,g(ω)is the phonon density of states, and N is the number of atoms in the crystal. If only considering the Umklapp phonon-phonon scattering
from the Matthiessen’s rule, the total phonon scattering rate
given as:

According to the phonon Boltzmann transport equation within relaxation time approximation, the thermal conductivity can be described as:

where C(ω) is the heat capacity, v(ω) is the phonon group velocity vector and τ(ω) is the effective phonon relaxation time. Thus we can have the following relationship:

where A and B are two constants independent of the density of vacancies.
The insert in Fig. 2 shows 1/κ as a function of the Li-vacancy concentration x (delithiation degree), which deviates notably from the linear relationship (dashed line in the inset) predicted by Eq. (6). This discrepancy indicates that there might exist other underlying mechanism responsible for the reduction of thermal conductivity. In the following part, we will show that delithiation induced ordered-todisordered transition plays an important role in this remarkable reduction of thermal conductivity.
3.2 Propagating and non-propagating heat energy transport
To shed light on the origin for the reduced thermal conductivity in LiCoO after delithiation, we explore the spectral 2 heat energy transmission coefficients in pristine LiCoO and the 2 delithiated LiCoO with different delithiation degrees following 2 Ref.30,31. The spectrum distribution of heat energy transmission T(ω) is calculated as31

where kB is the Boltzmann parameter and ΔT is the temperature difference between the two Langevin thermostats. Here q(ω) is the frequency dependent heat flux across the imaginary cross-section (green dashed line in Fig. 1a), which can be calculated as:30,31

where S is the cross area, t is the simulation time, and F is the inter- s ij atomic force on atom i due to atom j. Here, “L” and “R” denote the left and right segment, respectively, which are located at two sides of the imaginary cross-section.

Fig. 3 Heat energy transmission coefficient of the pristine LiCoO 2 and delithiated LiCoO . The heat energy transmission coefficient of 2 pristine LiCoO (solid line) and the delithiated LiCoO with 2 2 delithiation degrees as 5% (dashed line) and 10% (dotted line) respectively
As shown in Fig. 3, we find that after the delithiation process, the heat energy transmission coefficients in the delithiated LiCoO 2 (dashed and dotted line) are significantly suppressed compared to that in pristine LiCoO (solid line) for nearly the entire frequency 2 range. For the delithiated LiCoO , a further decrease is observed 2 when the delithiation degree increases from 10% to 40%. In the phonon-dominated heat transport, low frequency phonons are less sensitive to vacancy32,33 scattering compared to high frequency phonons, within the frame of perturbation theories by Ratsifaritana and Klemens.29 Unexpectedly, in LiCoO , the reduction in heat 2 energy transmission in low frequency regime is even more remarkable than that in the high frequency regime.
In crystals, phonons are the primary carriers of heat energy, which can be described by the kinetic theory such as Boltzmann transport equation.34 In contrast, in disordered solids, phonons are illdefined due to the lack of crystalline periodicity, which leads to the loss of definition of the wave vector.35-39 In 1990s, Allen et al.40 classified the heat carriers in disordered solids as propagons, diffusions, and locons, depending on the degree of localization of each carrier. Simply speaking, propagons are phonon-like carriers in a manner analogous to uncoupled particles, while diffusions carry heat in a random-walk manner. These two are delocalized and the primary heat carriers in disordered solids because locons are localized and nonpropagating.35 In delithiated LiCoO , the 2 delithiation disrupts the lattice periodicity, and introduces degree of disorder. Therefore the delithiation of LiCoO can give rise to 2 diffusions associated with the phonon-vacancy scattering, and suppress the heat energy transmission in both high and low frequency regimes. As low frequency vibrational modes usually carry more heat energy than high frequency modes, because of the higher density of modes and higher speed of sound, this results in the unexpectedly large reduction in thermal conductivity of LiCoO in 2 delithiation process.
3.3 Temperature dependence in thermal conductivity and anharmonic effect
So far, we have investigated the delithiation effect on the thermal transport of the LiCoO at 300 K and the underlying mechanism is 2 revealed by analyzing the phonon transmission coefficient. To further study the temperature effects on the thermal transport properties of LiCoO and delithiated LiCoO , we compute the 2 2 temperature dependence of κ and κ with different delithiation P D degrees.
Fig. 4 shows the temperature dependence of the thermal conductivity in pristine LiCoO and delithiated LiCoO . Obviously, 2 2 all samples exhibit a temperature dependence that does not follow the typical expectation from Debye-Callaway model (T-1), which suggests the presence of a disordering event. The results show that the κ and κ are decreased at elevated temperature and follow T -α P D behavior due to the intrinsic Umklapp phonon scattering.41And the fitting values of α are 0.53, 0.49, 0.22 and 0.02 corresponding to the pristine and delithiated LiCoO with percentage of 5%, 10%, and 2 40%, respectively. In other words, comparing with the pristine LiCoO , the thermal conductivity of the delithiated LiCoO is less 2 2 sensitive to the variation of temperature, especially when the delithiation degree is high.
For the pristine LiCoO , the thermal conductivity follows T-0.53 2 by fitting, which shows a deviation from the Debye-Callaway model (T-1).42 The possible reason is given as follows. In bulk system, when the phonon-phonon scattering dominates the thermal transport, the

Fig. 4 The temperature effect on the thermal transport properties of LiCoO and delithiation LiCoO . Log-log plot for the temperature 2 2 dependence of thermal conductivity of LiCoO and delithiated 2 LiCoO with different delithiation degrees, showing a T -α behavior. 2 The symbols are simulation results and the solid lines are the linear fitting lines. The error bars are standard deviation of 3 independent simulations with different initial conditions.
temperature dependence of the thermal conductivity will follow the Debye-Callaway model,42 that is, the exponent α equals to 1. However, in our simulations, we fixed the system length L as 32 nm (quite short compared to the phonon mean free path), the boundary scattering should play a significant role in the thermal transport. Our previous work28 has been verified that the value of α is sensitive to the total length of the simulations domain. And as the length L increases (the boundary effects subside), the value of α will be increased. Hence, for the pristine case, the deviation in α from 1 is due to the boundary scattering causing by the very small system length setting in our simulations.
Anharmonic phonon scattering plays dominated role in temperature dependent thermal conductivity. To understand the underlying physical mechanism for the distinct temperaturedependence of pristine and delithiated LiCoO , we have calculated 2 the spectral energy density (SED) from atom trajectories, The spectral energy density can be calculated as43,44

where k is the wave vector, p is the phonon polarization branch,is the atomic mass, N is the total number of unit cells, is the velocity of the component of the jth atom in the lth unit cell which is dumped from MD simulations, is the equilibrium position of the atom, and is the total integration time. Here, for both pristine and delithiated LiCoO , the spectral energy density is computed along x- 2 axis direction. j m
Fig. 5 shows the low frequency part of phonon dispersion curves (< 5THz) in the pristine and delithiated LiCoO at 300 K. The 2 shading on the plot indicates the magnitude of the spectral energy density for each phonon mode.43 And the magnitude increases when the color varies from blue to red. In the pristine LiCoO (Fig. 5a), the 2 longitudinal acoustic (LA) branch and two transverse acoustic (TA ,

Fig. 5. The spectral energy density of pristine LiCoO and delithiation LiCoO . Spectral energy density of the pristine LiCoO (a) and the 2 2 2 delithiated LiCoO with delithiation degree of 5% (b) and 10% (c), respectively.
TA ) branches are very clear. And a slight broadening of SED curve 2 is observed, which originates from the temperature-dependent anharmonic phonon-phonon interactions. 43 In contrast, a more pronounced broadening of SED curve is observed in the delithiated LiCoO (Fig. 5b-c), and the magnitude of SED is substantially 2 reduced compared to that in the pristine LiCoO . For each phonon 2 mode, the broadening range of SED, which is related to the phonondefect scattering intensity, directly reflects the phonon relaxation time.44 The larger broadening in SED indicates a smaller phonon relaxation time. Therefore, the broadening and reduced magnitude of SED curve provide a direct evidence that the phonon relaxation times in delithiated LiCoO are reduced compared to the pristine 2 case due to the strong phonon-defect scatterings.45,46 Moreover, there exists a further broadening and magnitude reduction, revealing a further reduction of the relaxation time, when the delithiation degree changes from 5% to 10%.
In previous works,47,48 McGaughey et al. has demonstrated that the thermal conductivity contributed by short mean free path (MFP) phonons is temperature independent, while the long MFP phonons dominate the temperature dependence behavior. In the delithiated LiCoO structure, the long MFP phonons are suppressed due to the 2 strong phonon defect scattering induced by the delithiation. Thus, the relative contribution from the short-MFP phonons to κ increases D compared with that in the pristine case, which leads to κ present a D less sensible dependence on temperature. When the delithiation degree increases, the long MFP phonons are further suppressed, and thus the temperature dependence becomes more insensitive correspondingly.
4. Summary
In summary, we have investigated the effect of delithiation on the thermal transport properties of LiCoO via molecular dynamics simulations. With the increase of delithiation degree, the thermal conductivity of delithiated LiCoO decreases monotonically, and 2 becomes less sensitive to temperature. The calculation results of the spectral heat energy transmission coefficient and the spectral energy distribution together reveal that the delithiation induces the nonpropagating heat transport, giving rise to the substantially decreased thermal conductivity in delithiated LiCoO . Our study provides 2 valuable guidance for the efficient thermal management in the Li-ion battery.
Acknowledgements
This project is supported in part by the grants from the National Key Research and Development Program of China (Grant No. 2017YFB0406000), the National Natural Science Foundation of China (Grant Nos. 51506153 and 11334007), and Shanghai Science and Technology Committee (Grant Nos. 18JC1410900 and 17ZR1448000). J. C. acknowledges support from the National Youth 1000 Talents Program of China. S. H. acknowledges support from the International Exchange Program for Graduate Students at Tongji University (No. 2017XKJC-010).
Reference