Article Cite This: J. Phys. Chem. C 2018, 122, 1447−1455
pubs.acs.org/JPCC
Tunable Anisotropic Thermal Conductivity and Elastic Properties in Intercalated Graphite via Lithium Ions Zhiyong Wei,*,† Fan Yang,‡ Kedong Bi,† Juekuan Yang,† and Yunfei Chen† †
Jiangsu Key Laboratory for Design & Manufacture of Micro/Nano Biomedical Instruments and Department of Mechanical Engineering, Southeast University, Nanjing 210096, People’s Republic of China ‡ The Molecular Foundry, Lawrence Berkeley National Laboratory, Berkeley, California 94720, United States S Supporting Information *
ABSTRACT: The electrochemical intercalation of metal ions into layered materials is a useful strategy to reversibly tune thermal-transport properties, but the fundamental mechanisms are not well understood. In this study, we systematically investigated the effects of lithium-ion concentrations on the anisotropic thermal conductivity of intercalated graphite by molecular-dynamics simulations and a continuum-mechanics method. It was found that the in-plane thermal conductivity rapidly decreased to 34.3% of that of graphite in the same direction when the lithium-ion concentration increases to 0.1667. However, the cross-plane thermal conductivity first decreased to 23.7% of that of the graphite. Then, it surprisingly recovered to a thermal conductivity even higher than that of graphite when the lithium-ion concentration increased further. These two different trends and thermal-conductivity anisotropies were explained by extracting the phonon lifetimes and elastic constants of intercalated graphites with various lithium-ion concentrations. At lithium-ion concentrations lower than 0.05, the reduction of both the in-plane and cross-plane thermal conductivities in the intercalated graphite was attributed to the increased phonon scattering owing to the interactions between the lattices and ions. However, at lithium concentrations higher than 0.05, the thermal transport of the intercalated graphite was mainly influenced by anisotropic elastic constants. The rapidly increasing cross-plane elastic constants increased the cross-plane thermal conductivity, while simultaneously weakening the phononfocusing effects along the in-plane direction, which resulted in the opposite tendencies of the in-plane and cross-plane thermal conductivities with increasing lithium-ion concentrations. This study provides important guidance in the active regulation of the anisotropic thermal conductivities of layered materials for thermal management in energy storage and conversion. thermoelectric generators,15 thermal cloakings,16 and energy storage.17 Note that most layered materials have large interlayer spaces (such as 3.35 Å for graphite), which are usually larger than the van der Waals radii of alkali metals.18 When the guest ions are driven into the interlayer space, the intercalation can induce some additional interactions between lattices and ions, which may influence both the kin and kc and thus the thermal anisotropy.19 When the guest ions are removed from a layered material, its thermal-transport property can be regained.20 Therefore, the intercalation of layered materials may provide an important route for reversibly tuning thermal-transport properties. Electrochemical intercalation is a common method to drive metal ions into and extract them from layered materials and is very similar to the charge or discharge processes of chemical batteries. By electrochemically intercalating lithium ions (Li ions) into layered MoS2, Zhu et al.19 measured both the inplane and cross-plane thermal conductivity of intercalated
1. INTRODUCTION Two-dimensional van der Waals materials, such as graphite, molybdenum disulfide, and black phosphorus, have been of much interest in recent years because of their unique lattice and electron properties.1−3 In particular, the thermal-transport properties of these layered materials have been extensively investigated.4−6 Previous investigations have shown that the number of layers can change the in-plane thermal conductivity of suspended graphene7,8 and encased graphene.9 However, when the number of layers is larger than five, the thermal conductivity of suspended graphene stays constant.7,8 Similarly, substrate coupling only influences the thermal transport of very few layers graphene.10,11 Aside from the number of layers, stress also influences the thermal conductivity. By applying stress, Chen and co-workers12 have shown that the cross-plane thermal conductivity of graphite can be controlled. However, few studies focus on tailoring the thermal anisotropy of bulk layered materials over a wide range, such as from unity to 30 or so.13,14 Thermal anisotropy (R) is generally defined as the ratio of the in-plane thermal conductivity (kin) to the cross-plane thermal conductivity (kc). Tunable thermal anisotropy has special advantages for various applications, such as in transverse © 2018 American Chemical Society
Received: September 30, 2017 Revised: December 7, 2017 Published: January 4, 2018 1447
DOI: 10.1021/acs.jpcc.7b09717 J. Phys. Chem. C 2018, 122, 1447−1455
Article
The Journal of Physical Chemistry C Table 1. Structure Parameters of Various Intercalated Graphites sample
graphite
LiC20
LiC18
LiC16
LiC12
LiC8
LiC6
x a (Å) c (Å) ρ (kg/m3)
0 2.42 6.78 2314
0.0500 2.42 7.48 2169
0.0556 2.42 7.49 2169
0.0625 2.42 7.53 2168
0.0833 2.42 7.59 2178
0.1250 2.41 7.48 2265
0.1667 2.41 7.27 2393
Figure 1. Preparation of the equilibrium atomic structure for intercalated graphite by annealing at 3000 K. (a) Initial atomic structure before annealing and (b) equilibrium structure after annealing.
2. MODEL AND METHOD OF CALCULATION 2.1. Sample Preparation. In order to investigate the effects of Li-ion concentrations on the thermal transport of intercalated graphite, the atomic structure of intercalated graphite with various Li-ion concentrations should be built. Previous studies have shown that combinations of lithium and graphite can form regular crystals,28 such as LiC6, LiC12, and LiC18. However, during the practical electrochemical intercalation process, the Li-ion concentration in intercalated graphite may change continuously, and the atomic structure strongly depends on the context of the thermodynamics.29 In order to prepare a reasonable atomic model, we first chose intercalated graphites of seven different Li-ion concentrations. Each Li-ion concentration, x, which is defined as the ratio of the number of atoms of Li to that of C, is shown in Table 1, with the special case of x = 0 for bulk graphite. Next, we built the atomic structure of bulk graphite with an initial layer distance of 3.35 Å and a carbon−carbon (C−C) distance of 1.42 Å by a homemade code according to the ordered atomic position in Cartesian coordinates. For the intercalated graphite, we placed the lithium atoms uniformly in the interlayer space of the graphite. We used the same interatomic potentials for both the sample preparation and thermal-conductivity calculations, as shown in next section. Because an annealing temperature of ∼3000 K is in the range of the optimal temperature for singlecrystal graphite,30 we then performed MD simulations at a temperature of ∼3000 K and a pressure of ∼0 under an isothermal−isobaric ensemble. A high annealing temperature is favorable because it increases the atomic kinetic energy, makes it easier to overcome the energy barrier, and causes equilibrium positions to be reached more quickly. After the simulation system was annealed to 300 K with a rate of about 10 K/ps, we finally obtained the equilibrium atomic structures of intercalated graphite, which were used to calculate the thermal conductivities and elastic constants below. Figure 1 compares the atomic structure of LiC6 before and after annealing. It is obvious that the Li ions in the interlayer space of graphite are rearranged after annealing, whereas the
MoS2 using time-domain thermoreflectance (TDTR). It was found that the thermal anisotropy reached a maximum at an appropriate Li-ion concentration. This result is counterintuitive, because increasing the disorder of the structure generally leads to low thermal anisotropy.21−23 Zhu et al. postulate that lithiation-induced structural and compositional disorder tends to reduce the phonon mean free path more along the crossplane direction than along the in-plane direction. Recently, Qian et al.24 investigated the thermal transport of intercalated graphite from equilibrium-molecular-dynamics (MD) simulations and found that the intercalation of Li ions has an anisotropic effect on tuning the kin and kc. Kang et al.20 also used the TDTR method to investigate intercalated black phosphorus and found that both kin and kc are reduced with increasing Li-ion concentrations. Their experimental results agree well with the proposed theoretical model at low Li-ion concentrations.20 However, the model’s predictions deviate from the experimental results at high Li-ion concentrations at all fitting-parameter options. As we know, high elastic anisotropy in most layered materials such as graphite is closely correlated with phonon focusing,25−27 which may significantly influence the thermal conductivity of an intercalated composite, although this is neglected in the above studies. In order to explore how elastic anisotropy and phonon focusing influence the anisotropic thermal transport of an intercalated composite, the lattice thermal conductivity of intercalated graphite with Li ions was systematically investigated in this study by continuously adjusting the Li-ion concentration. MD simulations show that kin decreases with increasing Li-ion concentrations, while the cross-plane thermal conductivity, kc, first decreases and then increases. The opposite tendencies of kin and kc at high degrees of lithiation make the active regulation of the thermal anisotropy over a wide range of carrier concentrations possible. The phonon lifetimes and elastic constants of intercalated graphites with various Li-ion concentrations are calculated, and the Li-ion-concentrationdependent anisotropic thermal-transport properties are well explained. 1448
DOI: 10.1021/acs.jpcc.7b09717 J. Phys. Chem. C 2018, 122, 1447−1455
Article
The Journal of Physical Chemistry C
Figure 2. In-plane (a) and cross-plane (b) thermal conductivity of bulk graphite by an equilibrium molecular-dynamics simulation. The bold red curves are the averages of six independent simulations. The directions of heat transfer are shown in the insets.
Although the ensemble average of the heat current is 0 in an equilibrium state, the instantaneous heat current is not 0 and can be obtained by the statistical-mechanics theory.37 At last, the thermal conductivity along the α direction (α = x, y, or z) was calculated from the Green−Kubo formula:36,38,39
crystal structure is still ordered. On the basis of the AB-stacked graphite, variations in the equivalent lattice constants a and c were also carefully checked for the intercalated graphites in Table 1. After annealing and thermal relaxation, it was found that c increased about 6.7% when the Li-ion concentration increased from x = 0 (bulk graphite) to x = 0.1667 (LiC6), which agreed well with previously reported value of about 10.5%.24 However, because the intralayer covalent bonds among the C−C are very strong, the average in-plane C−C distance changed by only about 0.4%. 2.2. Molecular-Dynamics Simulation. The MD simulations were performed using the LAMMPS package.31 The AIREBO potential 32 was used to describe the C−C interactions, including the covalent bonds and nonbond van der Waals interactions. When Li ions were introduced into the interlayer space, the carbon atoms were negatively charged in order to neutralize the positive Li ions. Thus, we used additional electrostatic-interaction potentials and LennardJones (L-J) potentials to describe the interactions between Li and C. For the electrostatic interactions, we set the charge for every Li ion to +0.822 ec,33 where ec is the charge of an electron. The charge value for the Li ions was obtained using Bader charge analysis on the basis of the first-principles calculations of Zhang et al.33 It was assumed that the transferred electrons were uniformly shared by all the carbon atoms in the graphite. Thus, the charge for each carbon atom could be evaluated to keep the neutrality of the simulation system. For the L-J interactions, the well depths and equilibrium distances between the Li and Li were 1.090 85 meV and 2.1835 Å,34 respectively. The corresponding values between the Li and C were 2.255 018 meV and 2.766 75 Å,34 respectively. In this study, the velocityverlet integral algorithm was used to solve Newton’s equations of motion, with a time step of 0.5 fs. After the equilibrium atomic structure for an intercalated graphite with a specified Li-ion concentration was obtained, an equilibrium-MD simulation35,36 was performed to calculate the thermal conductivity. This equilibrium method was based on linear-response theory, which did not require the application of nonequilibrium heat currents across the super models. The simulation sizes along the three principal directions were about 24.6 × 25.6 × 26.8 Å. During the MD simulation, the system was first run for 250 ps to reach equilibrium at 300 K. Then, the temperature adjustment was closed, and the simulation system was further run at a constant number-of-particles, volume, and energy (NVE) ensemble for 40 ns, during which the heatcurrent components along the three directions were recorded.
kα =
1 VkBT 2
∫0
τm
⟨Jα (τ )Jα (0)⟩dτ
(1)
where V was the volume of the simulation system, kB was the Boltzmann constant, T was the system temperature, τm was the integration time, and Jα was the heat-current component. The angle brackets denoted autocorrelation, so that the bracketed term was the heat-current autocorrelation function (HCACF). For every atomic structure, several independent simulations were performed by changing the seed of the initial velocity to verify the simulation finally converged. To reduce the simulation uncertainty, a total of six simulations were performed and averaged to obtain the present thermal conductivity.
3. RESULTS AND DISCUSSION 3.1. Calculation of Anisotropic Thermal Conductivity. To obtain both the in-plane and cross-plane thermal conductivity of graphite, their convergences are checked first, as shown in Figure 2. It is found that the in-plane thermal conductivity of bulk graphite, kin(x = 0), converges very quickly to reach about 381.2 W m−1 K−1. This result is much lower than the experimental result (about 2000 W m−1 K−1),40 although it is nearly the same as that of the equilibrium-MD simulation on bulk graphite (between 300 and 400 W m−1 K−1) by Oliveira and Greaney41 with the same interaction potential. Because an AIREBO potential can accurately predict the inplane velocity of graphite,32 this may imply that it significantly underestimates the phonon lifetime of bulk graphite because of the strong nonharmonic interactions.32 However, Figure 2b shows that the obtained cross-plane thermal conductivity of bulk graphite, kc(x = 0), is about 8.0 W m−1 K−1, which is very close to the experimental result (6.8 W m−1 K−1).42 The effects of size on the thermal conductivity of graphite are also checked. The present simulation size (24.6 × 25.6 × 26.8 Å) has been extended by eight times (2 × 2 × 2). The calculated in-plane and cross-plane thermal conductivities with the larger sizes (371.7 and 7.6 W m−1 K−1) are nearly the same as the present values (381.2 and 8.0 W m−1 K−1). Therefore, the current simulation box is large enough to avoid the size effects. 1449
DOI: 10.1021/acs.jpcc.7b09717 J. Phys. Chem. C 2018, 122, 1447−1455
Article
The Journal of Physical Chemistry C
3.2. Phonon-Lifetime Analysis. It is generally accepted that lattice thermal conductivity is related to two major factors: the phonon lifetime and the phonon-group velocity. In order to elaborate the effects of Li-ion concentrations on the thermal conductivity of intercalated graphite, we will discuss the phonon lifetime in this subsection and the group velocity effect in the next one. Although spectral-energy density (SED) is a common method to qualify phonon lifetimes,48 our previous investigation49 of single-layer graphene shows that this method is difficult to apply here because of the large number of atomic positions and velocity data. For distinguishing the acoustic and optical modes and considering that the decay of HCACF from the above equilibrium-MD simulations is directly related to the phonon scattering and has been verified as an effective method to extract phonon lifetimes in the literature,38,50,51 we evaluate the average phonon lifetime of intercalated graphite by fitting the HCACF with doubleexponential-decay functions Y(t), as suggested by Che et al.:51
We know that the present AIREBO potential underestimates the in-plane thermal conductivity. We still use this force field for several reasons. First of all, we use the conventional method43 to normalize the thermal conductivity to the bulk intrinsic thermal conductivity (at x = 0) in both directions. This can significantly reduce the effects of underestimation. In addition, few force fields can accurately predict both the inplane and cross-plane thermal conductivities of graphite at the same time.44,45 Among all of them, the AIREBO potential is still the most popular potential in the literature.46,47 Finally, our main objective in this study is to explore the variation of anisotropic thermal conductivities in graphite with Li-ion concentrations and not to obtain the most accurate thermal conductivity. The relative in-plane and cross-plane thermal conductivities of the intercalated graphites with various Li-ion concentrations are shown in Figure 3a,b. It is found that the kin/kin(x = 0)
Y (t ) = A a ·e−t / τa + Ao ·e−t / τo
(2)
where Ai and τi (i = a or o) are the fitting parameters, and the subscripts a and o denote slow acoustic modes and fast optical modes, respectively. Before double-exponential fitting, it is necessary to carry out low-pass filtering for the raw HCACF, if a high-frequency component exists within it, to reduce the numerical noise. A high-frequency component generally results from the optical phonons around the Γ point of the first Brillouin zone (FBZ).52 Although it does not influence the calculated thermal conductivity,53 the high-frequency component makes the HCACF oscillate rapidly, as it is shown in the left inset in Figure 4. The normalized HCACF after low-pass filtering is also
Figure 3. Relative in-plane (a) and cross-plane (b) thermal conductivities of the intercalated graphites with various Li-ion concentrations. (c) Corresponding thermal anisotropy for these samples. The solid and empty blue circles are the modeling results of Qian et al.24 from a different interaction potential.
continuously decreases with increasing Li-ion concentrations, whereas the kc/kc(x = 0) first decreases and then increases sharply with increasing the Li-ion concentrations. Our calculations show that the in-plane thermal conductivity of the intercalated graphite can be as low as 34.3% of that of bulk graphite. Additionally, the variation of the cross-plane values is in the range of 23.7 to 104.7% when the Li-ion concentration is adjusted. We also show the modeling results of Qian et al.24 in Figure 3a,b. Note that the exact thermal conductivities of their results are not the same as those in this study because of the different interaction potentials used in the MD simulations, but the tendency of the thermal conductivity is similar to that of the two calculations with Li-ion concentrations. The thermal anisotropy, R, of intercalated graphite, which is defined as R = kin/kc, is shown in Figure 3c. It is found that R can be actively regulated from ∼15 to ∼100. This thermal-conductivity anisotropy agrees with Zhu and colleagues’ experimental observations on both the trend and values (∼52 to ∼110) of intercalated MoS2 with Li ions.19 At low concentrations (x), the increased thermal-conductivity anisotropy with respect to x is counterintuitive. When the concentration increases further, the anisotropy decreases, which is normal. Zhu et al.19 postulate a possible mechanism through a combination of phonon focusing and different length scales of disorder between the in planes and cross planes. In our work, by analyzing the thermal conductivities in both directions, as explained in the next few subsections, we conclude that the mechanism for thermal anisotropy is phonon focusing.
Figure 4. Phonon lifetimes of intercalated graphites with various Liion concentrations. The solid red squares are obtained by doubleexponential-fitting independent heat-current autocorrelation functions (HCACF) after filtering is performed. The normalized HCACF before and after the filtering of the high-frequency oscillations are shown in the insets.
shown in this inset. After the filtering, the oscillation is removed, and the HCACF becomes monotonic. Figure 4 shows the fitted phonon lifetimes for various Li-ion concentrations. By comparing the results between bulk graphite and the other intercalated graphites, it is found that the introduction of Li ions significantly reduces the average phonon lifetime. Thus, in the low-concentration regime in Figure 3a,b, the reduction of kin and kc is attributed to the reduction of the phonon lifetime caused by the lattice−ion interactions. However, the reduced 1450
DOI: 10.1021/acs.jpcc.7b09717 J. Phys. Chem. C 2018, 122, 1447−1455
Article
The Journal of Physical Chemistry C lifetime alone cannot explain the opposite trends of kin and kc in the high-concentration regime in Figure 3a,b, especially considering the phonon lifetime is not obviously dependent on the concentration, x, in this regime. Thus, the group-velocity effect is discussed in the following subsection. 3.3. Calculation of the Elastic Constant. In order to explain the opposite tendencies of kin and kc at high Li-ion concentrations, we should check the variation in the phonongroup velocities, which are another important factor that influences thermal transport, according to phonon kinetic theory.54 In crystals, phonon-group velocities are strongly related to the elastic constants of the materials. Thus, the elastic constants are the key physical quantities needed to identify variation in the phonon-group velocities that can be calculated from MD simulations, as demonstrated by our previous work.17 During the MD simulation, it is assumed that the whole atomic structure of the intercalated graphite is still in a hexagonal structure. For this structure, there are a total of six independent elastic constants, Cij. For each atomic structure prepared by high-temperature annealing, the isothermal− isobaric ensemble is first run to readjust the simulation system to 300 K and no pressure. Then, one specified strain is controlled over a narrow range while the other five strains are kept at 0. Meanwhile, the corresponding stress is calculated from the MD simulation.55 Lastly, each component of the elastic-constant tensor is obtained by linearly fitting the relations between the applied strain and the corresponding stress, according to Hooke’s law.56 Figure 5 shows the fitting
plane direction as C33/ρ and C44 /ρ ,57,58 respectively, where ρ is the density. Thus, the increased kc in Figure 3b at high Liion concentrations can be attributed to the increased crossplane elastic constants.59−61 This indicates that the phonongroup velocity indeed plays a role in the thermal transport of intercalated graphite. Interestingly, the tensile elastic constant C11 and shear elastic constant C12 and C66 first decrease and then increase slightly with increasing Li-ion concentrations, as shown in Figure 6a,b,f. C11 and C66 are related to the longitudinal and transverse sound velocities along the in-plane direction. Moreover, the in-plane elastic constant has a minimum at a concentration around x = 0.0625. In the present study, our MD simulation shows the Li-ion-concentrationdependent elastic constants and sound velocities. A better method to explain this dependence is the use of the electronic structure from the first-principles calculation, because the intercalation of Li ions can induce charge transfers. However, that is beyond the scope of this work. Our main objective is to explain the variation in anisotropic thermal conductivities with different Li-ion concentrations. Thus, we left the first-principles calculation of the electronic structure for a future work. 3.4. Phonon-Focusing Effects. It should be noted that the thermal conductivity is a collective property of all the phonons in the FBZ. In anisotropic materials or structures, such as graphite, the principal thermal conductivity (kin or kc) is determined not only by the phonons in the high-symmetry directions (the in-plane or cross-plane directions) but also by the phonons in all the other low-symmetry directions of the FBZ.25,62 This means that not only does the phonon-velocity amplitude affect the thermal conductivity, but the velocity direction does as well. This velocity-direction-related effect is called the phonon-focusing effect.27 For instance, our previous investigation25 has shown that even if only the interaction strength along a specified direction is changed, the thermal conductivity in its orthogonal direction is also changed, even though the interaction strength along the orthogonal direction stays constant. From the calculated elastic constants in Figure 6, it is found that from LiC16 to LiC6, the C33 and C44 of the crossplane direction increase by factors of ∼20 and ∼30, respectively. However, the C11 and C66 of the in-plane direction change by only about 10.6 and 7.9%, respectively. The large variation in the cross-plane properties (C33 and C44) should play an important role in the in-plane thermal conductivity, kin, because of phonon-focusing effects. The acoustic iso-energy surface is equivalent to the phonondispersion relation under a continuum approximation. The difference when comparing it to phonon dispersion is that it contains the phonon property along all directions, but the phonon dispersion usually only shows in some highly symmetrical directions. Thus, the acoustic iso-energy surface is more suitable to the study of how a low-symmetry-direction phonon would affect the high-symmetry-direction thermal conductivity, such as in the current case of graphite. In order to elaborate the phonon-focusing effects, the acoustic iso-energy surfaces of the intercalated graphites with various Li-ion concentrations must be calculated, because they provide direct demonstrations of how the cross-plane elastic constants influence the in-plane thermal conductivity. The method of calculating the iso-energy surfaces is in the Supporting Information. Figure 7 shows the projection of the iso-energy surface of intercalated graphite on the x−o−z and x−o−y planes. The specified frequency in all the panels of Figure 7 is
Figure 5. Typical fitting process of elastic-constant components by controlling the tensile strain εxx (a) and shear strain εxz (b). The unfitted, scattered points in (b) change slightly and irregularly with varied εxz, implying that the corresponding elastic-constant components are nearly 0.
process of some typical elastic-constant components by applying the tensile strain εxx and shear strain εxz. It is found that the three principal stresses, σxx, σyy, and σzz, increase linearly with increasing εxx in Figure 5a. Similarly, the shear stress σxz also increases linearly with increasing strain εxz. However, the shear stresses σxy and σyz change slightly and irregularly, which may imply that their corresponding elastic components are nearly 0 and agrees with the above hexagonalstructure assumption. Figure 6 shows the variation of the six elastic-constant components with different Li-ion concentrations in intercalated graphites. It is found that the tensile elastic constant C33 and shear elastic constants C13 and C44 rapidly increase with increasing Li-ion concentrations. C33 and C44 are related to the longitudinal and transverse sound velocities along the cross1451
DOI: 10.1021/acs.jpcc.7b09717 J. Phys. Chem. C 2018, 122, 1447−1455
Article
The Journal of Physical Chemistry C
Figure 6. Six extracted, independent elastic constants of intercalated graphites with various Li-ion concentrations. (a,b,f) In-plane elastic constants C11, C12, and C66, respectively. (c−e) Cross-plane elastic constants C13, C33, and C44, respectively.
Figure 7. Iso-energy surface projections of bulk graphite (a,e), LiC16 (b,f), LiC12 (c,g), and LiC6 (d,h) on the x−o−z plane (first row) and x−o−y plane (the second row) in the reciprocal space. The units for the x, y, and z axes are all 109 m−1. In all panels, the blue, green, and red curves represent the longitudinal-acoustic (LA), transverse-acoustic (TA), and flexible-acoustic (ZA) branches, respectively. The phonon-group-velocity vectors are always normal to the iso-energy surfaces and are represented by black arrows, which are gradually shifted along the cross-plane direction from (b) to (d).
Different from the projections of the iso-energy surface on the x−o−y plane, the projections on x−o−z show complex, highly anisotropic shapes and include regions with nearly no curvature, which are labeled by black arrows in Figure 7a−d. The phonon-velocity vectors in such regions have larger components along the in-plane directions. Therefore, Figure 7a−d shows that phonon focusing exists along the in-plane direction. With increasing Li-ion concentrations from LiC16 to LiC6, it is clear that the extent of each iso-energy surface in the cross-plane direction shrinks. For example, the TA branch on the x−o−z plane is roughly an ellipsoid revolving with long axis and parallel to the cross-plane direction. Increasing the Li-ion concentration reduces the ellipsoid’s z-direction diameter while leaving the x-direction diameter essentially constant. Simultaneously, increasing the Li-ion concentration gradually shifts the phonon-group-velocity vector along the cross-plane direction, as clearly demonstrated in Figure 7b−d, because the direction of the phonon-group velocity is always normal to the iso-energy surface. This weakens the phonon focusing along the in-plane direction and can reduce the in-plane thermal conductivity.25
0.1 THz. Generally, a small wavevector on the iso-energy surface means a higher sound speed under the long wavelength limit. Because the sound speed of the longitudinal mode is higher than that of the transverse mode, the blue solid lines in all the panels of Figure 7 are the longitudinal-acoustic (LA) branches. The green and red solid lines are the transverseacoustic (TA) branches and flexible-acoustic (ZA) branches, respectively. The iso-energy surface of bulk graphite in Figure 7a is very similar to that from our previous results calculated by the lattice-dynamics method,25 which validates this calculation. For the iso-energy surface, as shown on the x−o−y plane, we found that all three acoustic-phonon branches are nearly circular. This means that the acoustic properties of intercalated graphite are isotropic in the in-plane direction. Additionally, the contour profiles of LA and TA only change slightly. This is because C11 and C66, which are directly related to the in-plane sound speed,57 also change very slightly, as shown in Figure 6a,f. Because C44 rapidly decreases with decreasing Li-ion concentrations, the contour profile of ZA on the x−o−y plane also changes in Figure 7e−h. 1452
DOI: 10.1021/acs.jpcc.7b09717 J. Phys. Chem. C 2018, 122, 1447−1455
Article
The Journal of Physical Chemistry C
trations, which could result in a decreased kin. This study theoretically uncovers the phonon-transport mechanism of intercalated graphite and also provides essential guidance for tuning anisotropic thermal conductivity in layered materials, through which potential applications in thermoelectrics, energy storage, and thermal devices may be found.
Therefore, at high Li-ion concentrations, the in-plane thermal conductivity decreases with increasing Li-ion concentrations. We know that Qian et al.24 used different potentials to investigate the effects of Li-ion concentrations on the thermal conductivities of intercalated graphites. The similarity between their results and ours is the Li-ion-concentration-dependentthermal-conductivity trend, as shown in Figure 3. The major difference is the explanation of the origin of the thermalconductivity variation. For the reduced in-plane thermal conductivity, Qian et al. attributed it to the decreased phonon velocity and phonon lifetime due to the phonon hybridization between the Li-ion vibrational modes and acoustic modes of graphite. However, in our work, the in-plane-thermalconductivity reduction is caused by phonon focusing. As shown in Figure 7, the energy-flux direction is normal to the iso-energy surface. When the Li-ion concentration increases, the iso-energy surface significantly changes from a flat ellipsoid shape in Figure 7a to a circle shape in Figure 7d. This means the phonon-focusing effect decreases. Thus, the velocity components increase in the cross-plane direction and decrease in the in-plane direction. This leads to the increase in the thermal conductivity in the cross-plane direction (kc) and the decrease in the in-plane direction (kin). We consider the increased elastic constants in the cross-plane direction to be the main reason for the increased cross-plane thermal conductivity. Qian et al.24 attributed it to the increased phonon-group velocity, which is essentially equivalent to our conclusion.
■
ASSOCIATED CONTENT
S Supporting Information *
The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcc.7b09717. Method for the calculation of the acoustic iso-energy surfaces of intercalated graphite with various Li-ion concentrations (PDF)
■
AUTHOR INFORMATION
Corresponding Author
*E-mail:
[email protected]. ORCID
Zhiyong Wei: 0000-0001-8686-451X Yunfei Chen: 0000-0002-8682-868X Notes
The authors declare no competing financial interest.
■
ACKNOWLEDGMENTS The authors acknowledge financial support from the National Key R&D Program of China (2017YFB0406000) and the Natural Science Foundation of Jiangsu Province (BK20160072).
4. CONCLUSIONS In summary, we investigated the anisotropic thermal conductivity in intercalated graphite by continuously adjusting the Li-ion concentration. The molecular dynamics simulations show that kin gradually decreases with increasing Li-ion concentrations, whereas kc first decreases and then rapidly increases with increasing Li-ion concentrations. Thus, we obtained wide thermal anisotropy (from ∼15 to ∼100) in intercalated graphite by adjusting the Li-ion concentration, a method which can also be applied to other layered materials. In order to uncover the fundamental physics that dominate the anisotropic thermal-transport properties, we first extracted the phonon lifetime by fitting HCACF with a double-exponential function. The calculations show that the introduction of Li ions into intercalated graphite significantly reduces the phonon lifetime because of the interactions between the crystals and ions. This reduction explains the reductions in kin and kc in the low-concentration regime as compared to the kin and kc in bulk graphite. To explain the thermal-conductivity enhancement in the high-concentration regime in the cross-plane direction, we then extracted the elastic constants of the intercalated graphites with various Li-ion concentrations by linearly fitting the relations between the applied strains and corresponding stresses. We found that the cross-plane tensile elastic constant C33 and shear elastic constant C44 rapidly increased with increasing Li-ion concentrations, which indicated that the phonon-group velocity played an important role at high Li-ion concentrations and could explain the increase in kc. However, the slight variation in the in-plane elastic constants meant that adjusting the Li-ion concentration has only limited effects on the in-plane-group velocity. By further calculating the acoustic iso-energy surfaces of the intercalated graphite at various Li-ion concentrations, we clearly demonstrated that the phonon focusing along the inplane direction was weakened with increasing Li-ion concen-
■
REFERENCES
(1) Novoselov, K. S.; Geim, A. K.; Morozov, S. V.; Jiang, D.; Zhang, Y.; Dubonos, S. V.; Grigorieva, I. V.; Firsov, A. A. Electric field effect in atomically thin carbon films. Science 2004, 306, 666−669. (2) Wang, Q. H.; Kalantar-Zadeh, K.; Kis, A.; Coleman, J. N.; Strano, M. S. Electronics and optoelectronics of two-dimensional transition metal dichalcogenides. Nat. Nanotechnol. 2012, 7, 699−712. (3) Qiao, J.; Kong, X.; Hu, Z.-X.; Yang, F.; Ji, W. High-mobility transport anisotropy and linear dichroism in few-layer black phosphorus. Nat. Commun. 2014, 5, 4475. (4) Balandin, A. A.; Ghosh, S.; Bao, W.; Calizo, I.; Teweldebrhan, D.; Miao, F.; Lau, C. N. Superior thermal conductivity of single-layer graphene. Nano Lett. 2008, 8, 902−907. (5) Yan, R.; Simpson, J. R.; Bertolazzi, S.; Brivio, J.; Watson, M.; Wu, X.; Kis, A.; Luo, T.; Hight Walker, A. R.; Xing, H. G. Thermal conductivity of monolayer molybdenum disulfide obtained from temperature-dependent raman spectroscopy. ACS Nano 2014, 8, 986− 993. (6) Lee, S.; Yang, F.; Suh, J.; Yang, S.; Lee, Y.; Li, G.; Choe, H. S.; Suslu, A.; Chen, Y.; Ko, C.; et al. Anisotropic in-plane thermal conductivity of black phosphorus nanoribbons at temperatures higher than 100 K. Nat. Commun. 2015, 6, 8573. (7) Ghosh, S.; Bao, W.; Nika, D. L.; Subrina, S.; Pokatilov, E. P.; Lau, C. N.; Balandin, A. A. Dimensional crossover of thermal transport in few-layer graphene. Nat. Mater. 2010, 9, 555−558. (8) Wei, Z.; Ni, Z.; Bi, K.; Chen, M.; Chen, Y. In-plane lattice thermal conductivities of multilayer graphene films. Carbon 2011, 49, 2653− 2658. (9) Jang, W.; Chen, Z.; Bao, W.; Lau, C. N.; Dames, C. Thicknessdependent thermal conductivity of encased graphene and ultrathin graphite. Nano Lett. 2010, 10, 3909−3913. (10) Chen, J.; Zhang, G.; Li, B. Substrate coupling suppresses size dependence of thermal conductivity in supported graphene. Nanoscale 2013, 5, 532−536.
1453
DOI: 10.1021/acs.jpcc.7b09717 J. Phys. Chem. C 2018, 122, 1447−1455
Article
The Journal of Physical Chemistry C (11) Wang, Z.; Xie, R.; Bui, C. T.; Liu, D.; Ni, X.; Li, B.; Thong, J. T. L. Thermal transport in suspended and supported few-layer graphene. Nano Lett. 2011, 11, 113−118. (12) Chen, J.; Walther, J. H.; Koumoutsakos, P. Strain engineering of kapitza resistance in few-layer graphene. Nano Lett. 2014, 14, 819− 825. (13) Wei, Z.; Wehmeyer, G.; Dames, C.; Chen, Y. Geometric tuning of thermal conductivity in three-dimensional anisotropic phononic crystals. Nanoscale 2016, 8, 16612−16620. (14) Song, N.; Jiao, D.; Cui, S.; Hou, X.; Ding, P.; Shi, L. Highly anisotropic thermal conductivity of layer-by-layer assembled nanofibrillated cellulose/graphene nanosheets hybrid films for thermal management. ACS Appl. Mater. Interfaces 2017, 9, 2924−2932. (15) Teichert, S.; Bochmann, A.; Reimann, T.; Schulz, T.; Dressler, C.; Toepfer, J. An oxide-based thermoelectric generator: Transversal thermoelectric strip-device. AIP Adv. 2015, 5, 077105. (16) Vemuri, K. P.; Bandaru, P. R. Geometrical considerations in the control and manipulation of conductive heat flux in multilayered thermal metamaterials. Appl. Phys. Lett. 2013, 103, 133111. (17) Wei, Z.; Yang, F.; Bi, K.; Yang, J.; Chen, Y. Thermal transport properties of all-sp(2) three-dimensional graphene: Anisotropy, size and pressure effects. Carbon 2017, 113, 212−218. (18) Bondi, A. Van der Waals volumes and radii. J. Phys. Chem. 1964, 68, 441−451. (19) Zhu, G.; Liu, J.; Zheng, Q.; Zhang, R.; Li, D.; Banerjee, D.; Cahill, D. G. Tuning thermal conductivity in molybdenum disulfide by electrochemical intercalation. Nat. Commun. 2016, 7, 13211. (20) Kang, J. S.; Ke, M.; Hu, Y. Ionic intercalation in twodimensional van der Waals materials: In situ characterization and electrochemical control of the anisotropic thermal conductivity of black phosphorus. Nano Lett. 2017, 17, 1431−1438. (21) Renteria, J. D.; Ramirez, S.; Malekpour, H.; Alonso, B.; Centeno, A.; Zurutuza, A.; Cocemasov, A. I.; Nika, D. L.; Balandin, A. A. Strongly anisotropic thermal conductivity of free-standing reduced graphene oxide films annealed at high temperature. Adv. Funct. Mater. 2015, 25, 4664−4672. (22) Luckyanova, M. N.; Johnson, J. A.; Maznev, A. A.; Garg, J.; Jandl, A.; Bulsara, M. T.; Fitzgerald, E. A.; Nelson, K. A.; Chen, G. Anisotropy of the thermal conductivity in GaAs/AlAs superlattices. Nano Lett. 2013, 13, 3973−3977. (23) Wang, Y.; Gu, C.; Ruan, X. Optimization of the random multilayer structure to break the random-alloy limit of thermal conductivity. Appl. Phys. Lett. 2015, 106, 073104. (24) Qian, X.; Gu, X.; Dresselhaus, M. S.; Yang, R. Anisotropic tuning of graphite thermal conductivity by lithium intercalation. J. Phys. Chem. Lett. 2016, 7, 4744−4750. (25) Wei, Z.; Chen, Y.; Dames, C. Negative correlation between inplane bonding strength and cross-plane thermal conductivity in a model layered material. Appl. Phys. Lett. 2013, 102, 011901. (26) Chen, Z.; Wei, Z.; Chen, Y.; Dames, C. Anisotropic Debye model for the thermal boundary conductance. Phys. Rev. B: Condens. Matter Mater. Phys. 2013, 87, 125426. (27) Wolfe, J. P. Imaging Phonons; Cambridge University Press: Cambridge, 1998. (28) Sacci, R. L.; Gill, L. W.; Hagaman, E. W.; Dudney, N. J. Operando NMR and XRD study of chemically synthesized LiCx oxidation in a dry room environment. J. Power Sources 2015, 287, 253−260. (29) Reynier, Y. F.; Yazami, R.; Fultz, B. Thermodynamics of lithium intercalation into graphites and disordered carbons. J. Electrochem. Soc. 2004, 151, A422−A426. (30) Knibbs, R. H.; Mason, I. B. Thermal expansion of pyrolytic and its variation due to non-alignment of crystallites. Nature 1964, 203, 58−60. (31) Plimpton, S. Fast parallel algorithms for short rang molecular dynamics. J. Comput. Phys. 1995, 117, 1−19. (32) Stuart, S. J.; Tutein, A. B.; Harrison, J. A. A reactive potential for hydrocarbons with intermolecular interactions. J. Chem. Phys. 2000, 112, 6472−6486.
(33) Zhang, C.; Lan, J.; Jiang, H.; Li, Y.-L. Stable and metastable structures in compressed LiC6: dimensional diversity. J. Phys. Chem. C 2016, 120, 10137−10145. (34) Rappe, A. K.; Casewit, C. J.; Colwell, K. S.; Goddard, W. A.; Skiff, W. M. UFF, A full periodic-table force-field for molecular mechanics and molecular dynamics simulations. J. Am. Chem. Soc. 1992, 114, 10024−10035. (35) McGaughey, A. J. H.; Kaviany, M. Quantitative validation of the Boltzmann transport equation phonon thermal conductivity model under the single-mode relaxation time approximation. Phys. Rev. B: Condens. Matter Mater. Phys. 2004, 69, 094303. (36) Volz, S. G.; Chen, G. Molecular-dynamics simulation of thermal conductivity of silicon crystals. Phys. Rev. B: Condens. Matter Mater. Phys. 2000, 61, 2651−2656. (37) Guajardo-Cuellar, A.; Go, D. B.; Sen, M. Evaluation of heat current formulations for equilibrium molecular dynamics calculations of thermal conductivity. J. Chem. Phys. 2010, 132, 104111. (38) Li, X.; Maute, K.; Dunn, M. L.; Yang, R. Strain effects on the thermal conductivity of nanostructures. Phys. Rev. B: Condens. Matter Mater. Phys. 2010, 81, 245318. (39) Schelling, P. K.; Phillpot, S. R.; Keblinski, P. Comparison of atomic-level simulation methods for computing thermal conductivity. Phys. Rev. B: Condens. Matter Mater. Phys. 2002, 65, 144306. (40) Klemens, P. G.; Pedraza, D. F. Thermal conductivity of graphite in the basal-plane. Carbon 1994, 32, 735−741. (41) Oliveira, L. d. S.; Greaney, P. A. Thermal resistance from irradiation defects in graphite. Comput. Mater. Sci. 2015, 103, 68−76. (42) Taylor, R. The thermal conductivity of pyrolytic graphite. Philos. Mag. 1966, 13, 157−166. (43) Mu, X.; Wu, X.; Zhang, T.; Go, D. B.; Luo, T. Thermal transport in graphene oxide−from ballistic extreme to amorphous limit. Sci. Rep. 2015, 4, 3909−3909. (44) Evans, W. J.; Hu, L.; Keblinski, P. Thermal conductivity of graphene ribbons from equilibrium molecular dynamics: Effect of ribbon width, edge roughness, and hydrogen termination. Appl. Phys. Lett. 2010, 96, 203112. (45) Ni, Y.; Chalopin, Y.; Volz, S. Significant thickness dependence of the thermal resistance between few-layer graphenes. Appl. Phys. Lett. 2013, 103, 061906. (46) Ong, Z.-Y.; Pop, E. Molecular dynamics simulation of thermal boundary conductance between carbon nanotubes and SiO2. Phys. Rev. B: Condens. Matter Mater. Phys. 2010, 81, 155408. (47) Pei, Q.-X.; Sha, Z.-D.; Zhang, Y.-W. A theoretical analysis of the thermal conductivity of hydrogenated graphene. Carbon 2011, 49, 4752−4759. (48) Thomas, J. A.; Turney, J. E.; Iutzi, R. M.; Amon, C. H.; McGaughey, A. J. H. Predicting phonon dispersion relations and lifetimes from the spectral energy density. Phys. Rev. B: Condens. Matter Mater. Phys. 2010, 81, 081411. (49) Wei, Z.; Yang, J.; Bi, K.; Chen, Y. Mode dependent lattice thermal conductivity of single layer graphene. J. Appl. Phys. 2014, 116, 153503. (50) Chen, J.; Zhang, G.; Li, B. How to improve the accuracy of equilibrium molecular dynamics for computation of thermal conductivity? Phys. Lett. A 2010, 374, 2392−2396. (51) Che, J. W.; Cagin, T.; Deng, W. Q.; Goddard, W. A. Thermal conductivity of diamond and related materials from molecular dynamics simulations. J. Chem. Phys. 2000, 113, 6888−6900. (52) Landry, E. S.; Hussein, M. I.; McGaughey, A. J. H. Complex superlattice unit cell designs for reduced thermal conductivity. Phys. Rev. B: Condens. Matter Mater. Phys. 2008, 77, 184302. (53) Qiu, B.; Ruan, X. Molecular dynamics simulations of lattice thermal conductivity of bismuth telluride using two-body interatomic potentials. Phys. Rev. B: Condens. Matter Mater. Phys. 2009, 80, 165203. (54) Chen, G. Nanoscale Energy Transport and Conversion; Oxford University Press: New York, 2005. (55) Thompson, A. P.; Plimpton, S. J.; Mattson, W. General formulation of pressure and stress tensor for arbitrary many-body 1454
DOI: 10.1021/acs.jpcc.7b09717 J. Phys. Chem. C 2018, 122, 1447−1455
Article
The Journal of Physical Chemistry C interaction potentials under periodic boundary conditions. J. Chem. Phys. 2009, 131, 154107. (56) Nye, J. F. Physical Properties of Crystals: Their Representation by Tensors and Matrices; Oxford University Press: Oxford, 1957. (57) Kelly, B. Physics of Graphite; Applied Science: London, 1981. (58) Prasher, R. Thermal boundary resistance and thermal conductivity of multiwalled carbon nanotubes. Phys. Rev. B: Condens. Matter Mater. Phys. 2008, 77, 075424. (59) Toberer, E. S.; Baranowski, L. L.; Dames, C. Advances in thermal conductivity. Annu. Rev. Mater. Res. 2012, 42, 179−209. (60) Slack, G. A. Nonmetallic crystals with high thermal conductivity. J. Phys. Chem. Solids 1973, 34, 321−335. (61) Jia, T.; Chen, G.; Zhang, Y. Lattice thermal conductivity evaluated using elastic properties. Phys. Rev. B: Condens. Matter Mater. Phys. 2017, 95, 155206. (62) Chen, Z.; Dames, C. An anisotropic model for the minimum thermal conductivity. Appl. Phys. Lett. 2015, 107, 193104.
1455
DOI: 10.1021/acs.jpcc.7b09717 J. Phys. Chem. C 2018, 122, 1447−1455
