Diffusion of alkali metals in the first stage graphite intercalation compounds by vdW-DFT calculations

Abstract

Diffusion of alkali metal cations in the first stage graphite intercalation compounds (GIC) LiC₆, NaC₆, NaC₈ and KC₈ has been investigated with density functional theory (DFT) calculations using the optPBE-vdW van der Waals functional. The formation energies of alkali vacancies, interstitials and Frenkel defects were calculated and vacancies were found to be the dominating point defects. The diffusion coefficients of the alkali metals in GIC were evaluated by a hopping model of point defects where the energy barriers for vacancy diffusion were derived from transition state theory. For LiC₆, NaC₆, NaC₈ and KC₈, respectively, the diffusion coefficients were found to be 1.5 × 10⁻¹⁵, 2.8 × 10⁻¹², 7.8 × 10⁻¹³ and 2.0 × 10⁻¹⁰ m² s⁻¹ at room temperature, which is within the range of available experimental data. For LiC₆ and NaC₆ a curved vacancy migration path is the most energetically favourable, while a straight pathway was inferred for NaC₈ and KC₈. The diffusion coefficients for alkali metal vacancy diffusion in first stage GICs scales with the graphene interlayer spacing: LiC₆ ≪ NaC₈ < NaC₆ ≪ KC₈.

---

Introduction

Alkali metal graphite intercalation compounds (AM-GICs) display a rich variety of phases with different compositions, which find applications as battery electrodes and in molten salt electrolysis. The Li-GIC system has been intensively studied for decades because of its use in rechargeable lithium ion batteries,^1–8 while Na-GICs are model materials for Na-ion battery anodes.^9–13 Interactions between sodium and graphite is also a key factor in the degradation of cathode linings in aluminium electrolysis cells.¹⁴ In addition, K-GICs have been shown to be superconducting.¹⁵ Diffusion of alkali metals in graphite is imperative to the performance of batteries and the degradation of cathode linings during aluminium electrolysis. The microscopic mechanisms of diffusion in AM-GICs is however not fully understood due to the lack of reliable theoretical and experimental methods. Quasielastic Neutron Scattering (QENS)^16,17 and nuclear magnetic resonance (NMR)¹⁸ have been used to derive the diffusion coefficients of Li in LiC₆ and K in KC₈. Electrochemical methods, including potentiostatic intermittent titration technique (PITT),¹⁹ galvanostatic intermittent titration technique (GITT)²⁰ and electrochemical impedance spectroscopy (EIS),^21,22 have also been employed to determine the diffusion coefficient of Li in graphite anode materials. The diffusion coefficient of Na in graphite has been measured using a Rapoport test²³ and thermogravimetric analysis of sodium uptake in graphite.²⁴ The experimentally measured diffusion coefficients are scattered over several orders of magnitude.

Previous density functional theory (DFT) studies of alkali metal diffusion in graphite have primarily focused on Li-GICs,^25–30 while the diffusion of Na in Na-GICs and K in K-GICs has been comparatively less studied and a systematic comparison between Li, Na and K diffusion in GICs is lacking.

Here we study the diffusion of Li, Na and K in first stage graphite intercalation compounds by DFT calculations where the weak van der Waals interactions are accounted for. The optPBE-vdW density functional has recently been shown to accurately describe the structure and energetics of AM-GICs,³¹ and this functional was adopted for this study. The energetics of alkali vacancies, interstitials and Frenkel defects in first stage AM-GICs were evaluated and a microscopic diffusion mechanism in AM-GICs was proposed. The energy barriers were calculated according to climbing image Nudged Elastic Band theory³² and the diffusion coefficients were derived from transition state theory.²⁵

Transition state theory

Crystalline materials above 0 K contain point defects and these defects are generally mobile and determine the microscopic mechanism and kinetics of solid state mass transport processes. The point defect diffusion coefficient can be obtained from a hopping model using the Einstein–Smoluchowki equation³³where a is the jumping distance, Γ_T the total jumping frequency and d the dimensionality of the diffusion process. Here d is equal to 2 representing a spatially restricted 2D diffusion process. The total jumping frequency is given by eqn (2)where t is total jumping time. The factor of 1/2 is due to the success probability of the jump at the intermediate state. n_i is the number of the ith possible elementary jump. ω_i is the ith elementary jump frequency. The total jumping frequency for a curved path in MC₆ and straight path in MC₆ and MC₈ can be derived by eqn (3) and (4) respectively, where the details of elementary jumps for each path will be described further below.

The mean jumping frequency, ω_i, can be derived by two approaches. One is based on the Arrhenius equation:

where ν₀ denotes the frequency pre-factor which has been often approximated by a constant value in the order of 10¹³ s⁻¹,³⁴ ΔE_h,i is the difference between the energy at an activated state (saddle point) and the initial equilibrium state and k_B is the Boltzmann's constant. Alternatively, the mean jumping frequency can be derived by evaluating the vibrational free energy difference, ΔF_vib,i, between an activated state and an initial state, as stated in eqn (3):²⁵where the vibrational free energies are obtained through phonon calculations as described by eqn (4):²⁵where v_j is the vibration frequency of the jth normal mode. Atomic jumping in a diffusion process may consist of several steps and the total jumping frequency is obtained from the individual jumps.²⁵

Computational details

Density functional theory (DFT) calculations were done with the VASP code,^35–39 using the van der Waals (vdW) functional optPBE-vdW.^40,41 The optPBE-vdW functional has been carefully evaluated and found to be the most suitable vdW functional for describing the energetics and structural properties of alkali metal graphite intercalation compounds.³¹ The projector augmented wave (PAW) method⁴² was used with the C_h (2s, 2p), Li_sv (1s, 2s), Na_sv (2s, 2p, 3s) and K_sv (3s, 3p, 4s) potentials supplied with VASP. Electron wave functions were expanded in plane waves up to a cutoff energy of 910 eV and the SCF convergence energy was set to 10⁻⁷ eV. The high cutoff energy of 910 eV is necessary for vdW DFT calculations to converge with hard carbon potentials, and the C_h potential has previously been found to most accurately describe graphite and higher order GIC.³¹ A 2^nd order Methfessel–Paxton (MP)⁴³ smearing of σ = 0.01 eV was used for the electronic energy level occupancy. The Brillouin zone was sampled with a 7 × 7 × 9 Γ-centred k-point mesh for the 2 × 2 × 2 LiC₆ supercell and a similar k-point density was used for other GICs supercells. For the energetics of defect structures the atomic positions were relaxed until the Hellmann–Feynman forces on the ions were smaller than 10⁻³ eV Å⁻¹, with the exception of the NaC₆ supercell with a Frenkel defect at the first site and the KC₈ supercell with Frenkel defects at both sites (Fig. 1). Relaxation of these three defect structures converged to the corresponding perfect structures, hence static calculations were performed to estimate the defect formation energies. The ground state energies of Li, Na and K metals were done with a 15 × 15 × 15 Monkhorst–Pack k-point mesh and with the same convergence criteria as above. Transition state searches were done using the climbing-image nudged elastic band (cNEB)³² method as implemented in VASP. For all images along the band, atomic positions were relaxed until the Hellmann–Feynman forces on the ions converged to below 0.05 eV Å⁻¹ with fixed lattice vectors. The frozen phonon method was used to calculate lattice vibrations under the harmonic approximation, as implemented in the Phonopy code.⁴⁴ The atomic positions in the supercell were displaced by 0.01 Å in the x, y or z direction to obtain the force constants. Phonon calculations were only performed for first stage LiC₆ supercells with one Li vacancy.

Fig. 1 MC₆ and MC₈ supercells with alkali vacancy, self-interstitial and Frenkel defects for an AM-GIC seen along the c-axis. There is one possible vacancy and one possible self-interstitial site, while there are two possible Frenkel sites for each structure.

Results

Point defect in AM-GICs

Due to the strong chemical bonds within the graphite layers only point defects within the alkali metal layers in the first stage alkali-metal graphite intercalation compounds were considered: alkali vacancies, interstitials and Frenkel defects. The different configurations of these defects in an AM-GIC with MC₆ and MC₈ in-plane superstructures⁴⁵ are illustrated in Fig. 1. There is one possible vacancy and one possible self-interstitial site, while there are two possible Frenkel sites for each structure. To simulate point defect concentrations of 0.125 and 0.056, 2 × 2 × 2 and 3 × 3 × 2 supercells were used, respectively.

The defect formation energies were calculated from the change in cohesive energy E_coh (binding energy corresponding to sublimation energy) according to eqn (5)–(7) for a vacancy, self-interstitial and Frenkel defect, respectively:

ΔE_f,vac = E_coh,MCX[vac] − E_coh,MCX + E_{coh,M bulk} (8)

ΔE_f,int = E_coh,MCX[int] − E_coh,MCX − E_{coh,M bulk} (9)

ΔE_f,Fre = E_coh,MCX[Fre] − E_coh,MCX (10)

M in MC_X denotes the alkali metal and [vac], [int], [Fre] represents vacancy, self-interstitial and Frenkel defect respectively. For example, the internal energy of a vacancy defect in LiC₆ is the energy of a LiC₆ supercell with a vacancy defect minus the energy of a perfect LiC₆ supercell plus the energy of Li metal. The enthalpy of formation for a point defect can be expressed as

ΔH_f = ΔE_f + pV ≈ ΔE_f (11)

since the pV term can be neglected for condensed materials. The calculated defect formation enthalpies in LiC₆, NaC₆, NaC₈ and KC₈ are summarized in Table 1.

Table 1 The calculated defect formation enthalpies for LiC₆, NaC₆, NaC₈ and KC₈ in eV. 2 × 2 × 2 supercells were used for NaC₆, NaC₈ and KC₈

	Vacancy defect	Interstitial defect	Frenkel defect
			Site 1	Site 2
a Relaxation to these three defect structures converged to the corresponding perfect structures, hence static calculations were performed to estimate the defect formation energy.
LiC6 222	0.166	0.460	0.397	0.651
LiC6 332	0.175	0.446	0.393	0.635
NaC6	−0.068	0.946	1.273^a	1.210
NaC8	−0.057	0.882	0.332	0.857
KC8	0.523	0.561	2.698^a	1.381^a

---

Vacancies have the lowest defect formation energies for all the first stage AM-GICs and will dominate over interstitials and Frenkel defects. Vacancy diffusion is thus proposed to be the most likely diffusion mechanism in AM-GICs. The negative formation enthalpies for vacancies in NaC₆ and NaC₈ imply that vacancies will form spontaneously. This is in agreement with our previous report: all the first stage Na-GICs are thermodynamically unstable.³¹

Defect formation enthalpies were also calculated for both 2 × 2 × 2 and 3 × 3 × 2 LiC₆ supercells in order to estimate the possible error from the relatively small size of the 2 × 2 × 2 supercells. The difference in formation enthalpies between the different supercell sizes were only 0.009–0.016 eV. Increasing the size of the supercell beyond 2 × 2 × 2 does not significantly affect the defect formation enthalpies. The following cNEB calculations to study the diffusion of alkali metal in GICs were therefore performed with 2 × 2 × 2 supercells due to the excessive computational cost of using larger supercells.

Diffusion of alkali metal in the graphite intercalation compounds

Minimum energy path and diffusion barrier. The possible migration pathways for alkali vacancy diffusion in MC₆ and MC₈ structures are illustrated in Fig. 2. The migration along the crystallographic c-axis through a hexagonal carbon ring is prohibited due to the high energy barrier.^29,30 The diffusion process in AM-GICs can thus be characterised as 2D diffusion or in-plane diffusion.^27,46 As illustrated in Fig. 2b, there are two possible in-plane migration paths in the MC₆ structure: a straight path going over the top of carbon atoms and a curved path travelling over carbon bridges and the hollow of a carbon ring. The distance h in Fig. 2b defines the curvature of the curved path. Only one straight migration path is reasonable in the MC₈ structure, as shown in Fig. 2c.

Fig. 2 Possible migration paths for vacancy diffusion in (a and b) MC₆ and (c) MC₈ structures. The distance h in Fig. 1b defines the curvature of a curved path.

The energetics of the in-plane diffusion pathways for a vacancy defect in LiC₆, NaC₆, NaC₈ and KC₈ are presented in Fig. 3 where the relative energies along the proposed diffusion paths are plotted against the distance from the initial vacancy site. The saddle points (transition states) are the highest points along the path. In the MC₆ curved path and the MC₈ straight path the saddle points are the positions where alkali metals are on top of carbon bridges. In the MC₆ straight path the saddle points are the positions of alkali metal on top of carbon atoms. There are two saddle point along the paths in all the structures, E₁ and E₂ denote the two corresponding energy barriers, where E₁ is the energy difference between the saddle point and the initial state and E₂ is the energy difference between the saddle point and the intermediate state at the middle point along the diffusion pathway.

Fig. 3 Relative energies along the in-plane diffusion pathways for vacancy defects in LiC₆, NaC₆, NaC₈ and KC₈. The points are the calculated values (images) while the lines are spline fitting curves.

The energies of the initial and final states are identical and set to 0 as a reference. E₁ and E₂ in the figures represent the energy differences between the saddle points and initial states and intermediate points, respectively. Due to symmetry, only one half of the pathways in the MC₈ structures were calculated.

The energy barriers E₁ and E₂ and the relative energies of the saddle and intermediate points for vacancy diffusion pathways in LiC₆, NaC₆, NaC₈ and KC₈ are summarized in Table 2.

Table 2 Energy barriers of a vacancy defect in LiC₆, NaC₆, NaC₈ and KC₈ structure, unit in [eV]

		LiC6	NaC6	NaC8	KC8
Straight	Saddle	0.57	0.30	0.30	0.13
	Intermediate	0.53	0.28	0.04	−0.01
	E1	0.57	0.30	0.30	0.13
	E2	0.04	0.02	0.26	0.14
Curved	Saddle	0.48	0.29	NA	NA
	Intermediate	0.23	0.21	NA	NA
	E1	0.48	0.29	NA	NA
	E2	0.25	0.08	NA	NA

---

The cNEB calculations show that the energy barrier for the first jump (E₁) is always higher than for the second jump (E₂), meaning that the first jump is the rate limiting step in the diffusion process. For diffusion in LiC₆ and NaC₆ the path with the lowest energy barriers is the minimum energy path (MEP). For LiC₆ the energy at the saddle point for the straight path is 0.09 eV higher than along the curved path. For NaC₆, this energy difference is only 0.01 eV.

The calculated energy profile for LiC₆ along the curved path in Fig. 3 is similar to the LDA (local density approximation) calculated by Toyoura et al.²⁵ As the concentration of alkali metal in a GIC increases with lower stage numbers Li–C interactions become progressively more important compared to vdW forces. This is why first stage GIC like LiC₆ can be described relatively well by conventional LDA or PBE GGA. For higher stage GICs vdW interactions will dominate progressively with increasing stage number and the use of vdW functionals is imperative.³¹ For consistency and comparison with higher stage GIC vdW functionals were also used in the present work on first stage GIC.

Jumping frequency and diffusion coefficient. The two energy barriers along the diffusion pathways imply that the diffusion process consists of two steps. As illustrated in Fig. 4, the first jump is from the initial site to the intermediate site, with 12 possible elementary jumps along a curved path, and with 6 possible jumps along a straight path. The second step is from the intermediate site to another vacancy site from which 2 possible elementary jumps exist.

Fig. 4 Possible elementary jumps for alkali metal atoms from an initial site to a vacancy site (final state) for (a) curved pathway and (b) straight pathway.

The diffusion coefficients were derived by evaluating the mean jumping frequency for each step and finally the total jumping frequency. Arrhenius plots of the diffusion coefficients of Li, Na and K in the first stage graphite intercalation compounds are shown in Fig. 5. The curved diffusion pathway was chosen for Li in LiC₆ and Na in NaC₆ since it leads to a higher diffusion coefficient compared to the straight pathway, which will be discussed further below. Experimental diffusivity data^16–24 are included for comparison.

Fig. 5 Calculated Arrhenius plots of vacancy defect diffusivity in first stage LiC₆, NaC₆, NaC₈ and KC₈. The results for LiC₆ and NaC₆ are calculated from curved diffusion pathways, while those for NaC₈ and KC₈ are from straight pathways. The points are experimental results from different methods: “”: Rapoport test,²³ “”: Na thermo-gravimetric test,²⁴ “”: QENS,^16,17 “”: NMR.¹⁸ The error bars show the data measured by electrochemical methods including PITT,¹⁹ GITT²⁰ and EIS.^21,22 All the electrochemical data were measured at room temperature. The positions of the error bars are shifted to improve the readability.

In general, the calculated diffusivities are in good agreement with the data acquired by Quasielastic Neutron Scattering (QENS)^16,17 and nuclear magnetic resonance (NMR).¹⁸ The data from electrochemical methods, as for the case of LiC₆, including potentiostatic intermittent titration technique (PITT),¹⁹ galvanostatic intermittent titration technique (GITT)²⁰ and electrochemical impedance spectroscopy (EIS),^21,22 are scattered over a relatively large range from 1.2 × 10⁻¹⁵ to 1.7 × 10⁻¹² m² s⁻¹. The data from electrochemical methods are scattered because the measurements cover a wide range of the state of charge (SOC): the concentration of Li in graphite varies during the measurements. The calculated diffusivities for NaC₆ and NaC₈ are also in good agreement with the measurements from Rapoport test²³ and Na thermo-gravimetric test.²⁴ However, deviations for NaC₆ and NaC₈ are also expected since these two compounds are thermodynamically unstable.³¹ The real Na concentration in graphite is expected to be considerable lower, corresponding to a higher stage number.⁴⁷

The activation energy (E_a) of a diffusion process can be determined from an Arrhenius plot according to eqn (9)

where D₀ is a pre-exponential factor. The diffusion coefficients D of Li, Na and K in graphite intercalation compounds at 25 °C and 960 °C, the activation energies (E_a) and the pre-exponential factors (D₀) are summarized in Table 3. The calculated activation energies are within the value range of the available experimental data.

Table 3 Diffusion coefficients D at 25 °C and 960 °C, activation energies (E_a) and pre-exponential factors (D₀) for LiC₆, NaC₆, NaC₈ and KC₈. Experimental data for activation energies are included for comparison

	LiC6	NaC6	NaC8	KC8
D @ 25 °C [m^2 s^−1]	1.5 × 10^−15	2.8 × 10^−12	7.8 × 10^−13	2.0 × 10^−10
D @ 960 °C [m^2 s^−1]	1.8 × 10^−9	1.0 × 10^−8	3.6 × 10^−9	1.3 × 10^−8
Activation energy [eV]	0.47	0.28	0.28	0.14
D0	2.29 × 10^−7	2.10 × 10^−7	5.24 × 10^−8	4.68 × 10^−8
Activation energy Ex [eV]	0.55 (ref. 18)	NA	NA	0.18 (ref. 17)
	1.0 (ref. 16)
	0.2 (ref. 48)

---

Discussion

The diffusion coefficients in Table 3 were derived by the Arrhenius approach according to eqn (1) and (5), where the empirical vibration pre-factor (v₀) was set to 10¹² and 10¹³, respectively,³⁴ giving diffusion coefficients within a range of one order of magnitude. The contribution from changes in the vibrational free energy, however, can also be elucidated through phonon calculations according to eqn (6) and (7). These two methods are compared in Fig. 6. It is clear that within the calculated temperature range, the results from the phonon calculations lie between the results from the Arrhenius approach with the two empirical vibration pre-factors. At room temperature, the line with v₀ equal to 10¹² is closer to the quantum statistics' result, while at elevated temperatures the line with v₀ equal to 10¹³ is closer. In general, the Arrhenius approach gives good estimates for the diffusion coefficients in the considered structures. Taking into consideration the scattered experimental data, computationally heavy phonon calculations would not provide additional insight. From the steeper slope in Fig. 6 it is clear that the activation energy inferred from quantum statistics is slightly higher than that from the Arrhenius approach; 0.50 eV and 0.48 eV, respectively.

Fig. 6 Vacancy defect diffusivity in LiC₆ structure as a function of inverse temperature. The solid line is calculated from quantum statistics; the dashed and dash-dotted lines are based on Arrhenius approach with pre-factors v₀ equal to 10¹² and 10¹³, respectively. Experimental data are also included for comparison with the same notations as in Fig. 5.

Curved vs. straight diffusion path in MC₆ structures

The diffusion coefficients of Li and Na in LiC₆ and NaC₆ obtained from straight and curved paths (Fig. 2) are listed in Table 4. The diffusion along a curved path is faster than along a straight path for both compounds, with a difference of more than one order of magnitude for Li in LiC₆ and a factor of three for Na in NaC₆. It is evident that the curved path is the minimum energy path (MEP) in both structures.

Table 4 Diffusion coefficients for vacancy diffusion at 25 °C for LiC₆ and NaC₆ structures. v₀ in eqn (2) is set to 10¹²

	LiC6	NaC6
Straight	3.2 × 10^−17 [m^2 s^−1]	9.5 × 10^−13 [m^2 s^−1]
Curved	1.5 × 10^−15 [m^2 s^−1]	2.8 × 10^−12 [m^2 s^−1]
h	1.20 Å	0.99 Å

---

The curved path does not pass exactly at the middle of the carbon ring, but is pushed towards the straight path due to the repulsive forces between adjacent alkali metal ions. The curvature can be quantified as h as illustrated in Fig. 2b. The lower the value of h the closer the curved path is to the straight path. The curvature h in Table 4 demonstrates that Li travels along a more curved path with h equal to 1.20 Å while Na follows a more straight path with h equal to 0.99 Å. The larger Na ions give rise to stronger repulsive forces than the smaller Li ions do. Consequently, a smaller difference in the diffusion coefficient between a curved and a straight path is found for Na in the NaC₆ structure.

Charge transfer and redistribution during vacancy migration

As demonstrated in our previous work, in GIC the valence s electron of alkali metal atoms are donated to carbon p_z orbitals where the electrons are delocalized.³¹ Alkali metal vacancies cause a redistribution of charge around the vacancy site and the charge density difference between the LiC₆ supercell with and without a Li vacancy is illustrated in Fig. 7a. The charge density difference along the central line in the graphene layer is shown in Fig. 7b. Compared to the total charge density this difference is subtle.

Fig. 7 (a) Charge density (ρ) between the LiC₆ perfect supercell and the supercell with a Li vacancy. Isosurface levels were set to be +0.01 e⁻ Bohr⁻³ (yellow) and −0.01 e⁻ Bohr⁻³ (blue). (b) Charge density (ρ) along the line in (a) through the graphene layer.

The diffusion of alkali metal ions in AM-GICs leads to a redistribution of charge around the moving ions. The redistribution of electrons close to the Fermi energy indicates how the chemical bonding is affected by the elementary jumps. The partial charge density difference between the initial state and the saddle point for Li in LiC₆ and K in KC₈ are shown in Fig. 8a and b, respectively. Only electrons with energy within the range of E_F − 2 eV to E_F were considered. The linear partial charge density differences for the graphene layer between the initial and saddle points are illustrated for LiC₆, KC₈, NaC₆ and NaC₈ in Fig. 8c–f, respectively.

Fig. 8 Partial charge density difference between a (a) LiC₆ and (b) KC₈ supercell with a vacancy at the initial site and at the saddle point. Only the electrons with the energy within the region of E_F − 2 eV to E_F were considered. Isosurface levels were set to ±0.0013 e⁻ Bohr⁻³ for LiC₆ and ±0.0003 e⁻ Bohr⁻³ for KC₈ where yellow is positive and blue is negative. Partial charge density through the graphene layer for a vacancy travelling along the line in Fig. 1(a) and (b) in (c) LiC₆, (d) KC₈, (e) NaC₆ and (f) NaC₈.

Compared to the small difference in partial charge density for a stationary vacancy in Fig. 7, a substantial partial charge density difference for an AM atom at the initial site and a saddle point is evident from Fig. 8. The electrons with energies from E_F − 2 eV to E_F can be exclusively designated to C 2p_z orbitals³¹ which form a π* band with delocalized electrons. The high mobility of the C 2p_z electrons close to E_F allows a substantial charge redistribution in the graphene layer when an AM ion migrates. The partial charge density differences in the graphene layer is smaller for the MC₈ crystal structure where the carbon coordination number is lower. The relative partial charge density differences between LiC₆ and NaC₆, and between NaC₈ and KC₈, reflect the different graphene interlayer distances.

From a structural point of view, the calculated diffusion coefficients of Li, Na and K in the first stage GICs correlate with the graphene inter-layer distance, as illustrated in Fig. 9. The largest ion, K⁺, gives the largest separation of the adjacent graphene layers and the highest diffusion coefficient. In Na-GICs, the interlayer distances are subtly larger in NaC₆ than in NaC₈ and the diffusion of Na in NaC₆ is slightly faster in NaC₆ than in NaC₈. The diffusivities in AM-GICs follow the sequence: KC₈ ≫ NaC₆ > NaC₈ ≫ LiC₆. The diffusion coefficients decrease with increasing activation energy and diffusion energy barrier (E₁), which corresponds to the first elementary jump of the diffusion process.

Fig. 9 Diffusion coefficients (at 25 °C, v₀ equal to 10¹²) for alkali metal vacancy diffusion as a function of graphene interlayer spacing.

Real carbon materials like graphite or amorphous carbon contain grain boundaries, surfaces and other imperfections. Although the energy barriers for diffusion can differ from crystalline bulk in amorphous carbon and in the vicinity of structural defects, the crystalline GIC studied in this work is the ideal model system for real carbon materials used in batteries and electrodes.

Conclusion

Alkali metal vacancies were shown to be the dominating type of point defects in first stage alkali metal intercalation compounds GIC from van der Waals density functional theory calculations. Diffusion of alkali metals in first stage GIC is thus inferred to occur through the diffusion of vacancies. In LiC₆ and NaC₆ a curved migration path parallel to carbon bridges in adjacent graphene layers was found to be the most energetically favourable migration pathway. A straight path crossing the carbon bridges in adjacent graphene layers is the most favourable migration pathway in NaC₈ and KC₈. The calculated diffusion coefficients are within the range of experimentally found values. The diffusivities of alkali metals in GIC were shown to be strongly correlated with the graphene inter-layer spacing, and scale as follows: KC₈ ≫ NaC₆ > NaC₈ ≫ LiC₆.

Acknowledgements

The present work was carried out in the project Durable Materials in Primary Aluminium Production (DuraMat), financed by the research Council of Norway, Hydro Primary Aluminum Technology, Sør-Norge Aluminium (Søral), and Elkem Carbon. Permission to publish the results is gratefully acknowledged. Computational resources were provided by NOTUR (The Norwegian Metacentre for Computational Science) through the projects nn2962k, nn9268k and nn9264k.

Notes and references

1. R. Yazami and P. Touzain, J. Power Sources, 1983, 9, 365–371.
2. N. Takami, A. Satoh, M. Hara and I. Ohsaki, J. Electrochem. Soc., 1995, 142, 371–379.
3. A. Mabuchi, K. Tokumitsu, H. Fujimoto and T. Kasuh, J. Electrochem. Soc., 1995, 142, 1041–1046.
4. J. B. Goodenough, J. Solid State Electrochem., 2012, 16, 2019–2029.
5. J. S. Filhol, C. Combelles, R. Yazami and M. L. Doublet, J. Phys. Chem. C, 2008, 112, 3982–3988.
6. J. Sangster, J. Phase Equilib. Diffus., 2007, 28, 561–570.
7. M. Winter, J. O. Besenhard, M. E. Spahr and P. Novak, Adv. Mater., 1998, 10, 725–763.
8. M. Endo, C. Kim, K. Nishimura, T. Fujino and K. Miyashita, Carbon, 2000, 38, 183–197.
9. Y. L. Cao, L. F. Xiao, W. Wang, D. W. Choi, Z. M. Nie, J. G. Yu, L. V. Saraf, Z. G. Yang and J. Liu, Adv. Mater., 2011, 23, 3155–3160.
10. R. Alcantara, P. Lavela, G. F. Ortiz and J. L. Tirado, Electrochem. Solid-State Lett., 2005, 8, A222–A225.
11. S. W. Kim, D. H. Seo, X. H. Ma, G. Ceder and K. Kang, Adv. Energy Mater., 2012, 2, 710–721.
12. D. A. Stevens and J. R. Dahn, J. Electrochem. Soc., 2000, 147, 1271–1273.
13. M. D. Slater, D. Kim, E. Lee and C. S. Johnson, Adv. Funct. Mater., 2013, 23, 947–958.
14. K. Tschöpe, A. Støre, A. Solheim, E. Skybakmoen, T. Grande and A. Ratvik, JOM, 2013, 65, 1403–1410.
15. R. A. Jishi and M. S. Dresselhaus, Phys. Rev. B: Condens. Matter Mater. Phys., 1992, 45, 12465–12469.
16. A. Magerl, H. Zabel and I. S. Anderson, Phys. Rev. Lett., 1985, 55, 222–225.
17. H. Zabel, A. Magerl, J. J. Rush and M. E. Misenheimer, Phys. Rev. B: Condens. Matter Mater. Phys., 1989, 40, 7616–7632.
18. J. Langer, V. Epp, P. Heitjans, F. A. Mautner and M. Wilkening, Phys. Rev. B: Condens. Matter Mater. Phys., 2013, 88, 094304.
19. M. D. Levi and D. Aurbach, J. Phys. Chem. B, 1997, 101, 4641–4647.
20. H. Yang, H. J. Bang and J. Prakash, J. Electrochem. Soc., 2004, 151, A1247–A1250.
21. N. Takami, A. Satoh, M. Hara and T. Ohsaki, J. Electrochem. Soc., 1995, 142, 371–379.
22. P. Yu, B. N. Popov, J. A. Ritter and R. E. White, J. Electrochem. Soc., 1999, 146, 8–14.
23. A. P. Ratvik, A. Støre, A. Solheim and T. Foosnæs, Light Met., 2008, 2014, 973–978.
24. Z. Wang, A. P. Ratvik, E. Skybakmoen and T. Grande, Light Met., 2014, 2014, 1239–1244.
25. K. Toyoura, Y. Koyama, A. Kuwabara, F. Oba and I. Tanaka, Phys. Rev. B: Condens. Matter Mater. Phys., 2008, 78, 214303.
26. K. Toyoura, Y. Koyama, A. Kuwabara and I. Tanaka, J. Phys. Chem. C, 2010, 114, 2375–2379.
27. K. Persson, V. A. Sethuraman, L. J. Hardwick, Y. Hinuma, Y. S. Meng, A. van der Ven, V. Srinivasan, R. Kostecki and G. Ceder, J. Phys. Chem. Lett., 2010, 1, 1176–1180.
28. K. Persson, Y. Hinuma, Y. S. Meng, A. Van der Ven and G. Ceder, Phys. Rev. B: Condens. Matter Mater. Phys., 2010, 82, 125416.
29. V. Meunier, J. Kephart, C. Roland and J. Bernholc, Phys. Rev. Lett., 2002, 88, 075506.
30. S. Thinius, M. M. Islam, P. Heitjans and T. Bredow, J. Phys. Chem. C, 2014, 118, 2273–2280.
31. Z. Wang, S. M. Selbach and T. Grande, RSC Adv., 2014, 4, 4069–4079.
32. G. Henkelman, B. P. Uberuaga and H. Jónsson, J. Chem. Phys., 2000, 113, 9901–9904.
33. H. Mehrer, Diffusion in Solids, Springer, 2006.
34. A. Van der Ven, G. Ceder, M. Asta and P. D. Tepesch, Phys. Rev. B: Condens. Matter Mater. Phys., 2001, 64, 184307.
35. P. E. Blöchl, Phys. Rev. B: Condens. Matter Mater. Phys., 1994, 50, 17953–17979.
36. G. Kresse and J. Hafner, Phys. Rev. B: Condens. Matter Mater. Phys., 1993, 47, 558–561.
37. G. Kresse and J. Hafner, Phys. Rev. B: Condens. Matter Mater. Phys., 1994, 49, 14251–14269.
38. G. Kresse and J. Furthmüller, Comput. Mater. Sci., 1996, 6, 15–50.
39. G. Kresse and J. Furthmüller, Phys. Rev. B: Condens. Matter Mater. Phys., 1996, 54, 11169–11186.
40. M. Dion, H. Rydberg, E. Schröder, D. C. Langreth and B. I. Lundqvist, Phys. Rev. Lett., 2004, 92, 246401.
41. G. Román-Pérez and J. M. Soler, Phys. Rev. Lett., 2009, 103, 096102.
42. G. Kresse and D. Joubert, Phys. Rev. B: Condens. Matter Mater. Phys., 1999, 59, 18.
43. M. Methfessel and A. T. Paxton, Phys. Rev. B: Condens. Matter Mater. Phys., 1989, 40, 3616–3621.
44. A. Togo, F. Oba and I. Tanaka, Phys. Rev. B: Condens. Matter Mater. Phys., 2008, 78, 134106.
45. M. S. Dresselhaus and G. Dresselhaus, Adv. Phys., 2002, 51, 1–186.
46. Graphite Intercalation Compounds I, ed. H. Zabel and S. A. Solin, Springer-Verlag Berlin Heidelberg, Germany, 1990.
47. R. C. Asher and S. A. Wilson, Nature, 1958, 181, 409–410.
48. H. Estrade, J. Conard, P. Lauginie, P. Heitjans, F. Fujara, W. Buttler, G. Kiese, H. Ackermann and D. Guérard, Physica B+C, 1980, 99, 531–535.

---