Ultrafast lithium diffusion in bilayer buckled graphene: A comparative study of Li and Na

TL;DR

This study reveals that curvature in bilayer buckled graphene significantly accelerates lithium diffusion, with anisotropic effects influenced by buckling direction and stacking, offering insights for advanced battery anode design.

Contribution

It provides a comparative analysis of Li and Na diffusion in buckled graphene, highlighting the impact of curvature and vibrational asymmetry on diffusion rates.

Findings

Enhanced lithium diffusion parallel to folding axis

Diffusion anisotropy depends on buckling direction and stacking

Li diffusion is faster than Na in buckled graphene

Abstract

The effect of the curvature of bilayer graphene on the interlayer diffusion of Li atoms is investigated using molecular dynamics simulations. A spectacular enhancement of the diffusion constant parallel to the folding axis is found. The ratio of the parallel to the perpendicular diffusion depends on the buckling direction and stacking type, and it increases with the degree of buckling. The strongest anisotropy is observed in the case of fixed zig-zag edges. A comparison with the interlayer diffusion of Na suggests that the strong asymmetry in the vibrational states of buckled graphene and also the smaller mass of Li are likely to contribute to the observed diffusion enhancement. This work opens a new pathway to develop highly-efficient anodes for rechargeable alkaline batteries.

Tables

Table 1. Table 1: Li diffusion coefficient at T = 300 𝑇 300 T=300  K, frequency factor, and activation energy in single-layer (1L) and AA-stacked bilayer (2L) flat graphene.

system	$D_{300}$, cm2/s	$Z$, cm2/s	$E$, eV	prior work
Li-1L	$6.6\cdot {10}^{-7}$	$4.7\cdot {10}^{-3}$	$0.23$	$E=0.28-0.33$ eV  12, 11, 13
AA-Li-2L	$6.0\cdot {10}^{-6}$	$1.6\cdot {10}^{-4}$	$0.085$	$D=4.4\cdot {10}^{-6}-7\cdot {10}^{-5}$ cm2/s 24, 25

Table 2. Table 2: Li diffusion constants at T = 300 𝑇 300 T=300  K, frequency factors, and activation energies in AA- and AB-stacked bilayer graphene. Here AC (ZZ) denotes the type of the edge that was kept fixed in simulations (the edge in the ∥ parallel-to \parallel direction). The number 90 (or 95) in the system name denotes the degree of buckling (i.e., 95 corresponds to d / L = 0.95 𝑑 𝐿 0.95 d/L=0.95 ).

system	$D_{300}^{\parallel }$, cm2/s	$Z_{\parallel }$, cm2/s	$E_{\parallel }$, eV	$D_{300}^{⟂}$, cm2/s	$Z_{⟂}$, cm2/s	$E_{⟂}$, eV
AA-Li-AC-95	$3.2\cdot {10}^{-6}$	$2.2\cdot {10}^{-4}$	$0.11$	$2.7\cdot {10}^{-6}$	$1.5\cdot {10}^{-4}$	$0.10$
AA-Li-AC-90	$6.4\cdot {10}^{-6}$	$1.4\cdot {10}^{-4}$	$0.079$	$2.1\cdot {10}^{-6}$	$0.98\cdot {10}^{-4}$	$0.10$
AA-Li-ZZ-95	$2.3\cdot {10}^{-6}$	$1.0\cdot {10}^{-4}$	$0.099$	$3.4\cdot {10}^{-6}$	$0.50\cdot {10}^{-4}$	$0.070$
AA-Li-ZZ-90	$2.0\cdot {10}^{-5}$	$2.8\cdot {10}^{-4}$	$0.068$	$5.7\cdot {10}^{-7}$	$0.83\cdot {10}^{-4}$	$0.13$
AB-Li-AC-95	$5.8\cdot {10}^{-6}$	$2.4\cdot {10}^{-4}$	$0.96$	$2.7\cdot {10}^{-6}$	$0.99\cdot {10}^{-4}$	$0.093$
AB-Li-ZZ-95	$1.7\cdot {10}^{-5}$	$2.0\cdot {10}^{-4}$	$0.063$	$6.1\cdot {10}^{-7}$	$0.48\cdot {10}^{-4}$	$0.11$

Table 3. Table S.1: Na diffusion coefficient at T = 300 𝑇 300 T=300  K, frequency factor, and activation energy in single-layer (1L) and AA-stacking bilayer (2L) flat graphene.

system	$D_{300}$, cm2/s	$Z$, cm2/s	$E$, eV	prior work
Na-1L	$6.5\cdot {10}^{-6}$	$1.5\cdot {10}^{-3}$	$0.14$	$E=0.135$ eV 9
AA-Na-2L	$2.5\cdot {10}^{-5}$	$2.0\cdot {10}^{-4}$	$0.072$

Table 4. Table S.2: Na diffusion constants at T = 300 𝑇 300 T=300  K, frequency factors, and activation energies in AA-stacking bilayer graphene. Here AC (ZZ) denotes the type of the edge that was kept fixed in simulations (the edge in the ∥ parallel-to \parallel direction). The number 90 (or 95) in the system name denotes the degree of buckling (i.e., 95 corresponds to d / L = 0.95 𝑑 𝐿 0.95 d/L=0.95 ).

system	$D_{300}^{\parallel }$, cm2/s	$Z_{\parallel }$, cm2/s	$E_{\parallel }$, eV	$D_{300}^{⟂}$, cm2/s	$Z_{⟂}$, cm2/s	$E_{⟂}$, eV
AA-Na-AC-95	$9.9\cdot {10}^{-6}$	$0.40\cdot {10}^{-4}$	$0.036$	$6.8\cdot {10}^{-6}$	$0.56\cdot {10}^{-4}$	$0.055$
AA-Na-AC-90	$9.9\cdot {10}^{-6}$	$1.7\cdot {10}^{-4}$	$0.074$	$3.0\cdot {10}^{-6}$	$2.3\cdot {10}^{-4}$	$0.11$
AA-Na-ZZ-95	$5.4\cdot {10}^{-6}$	$0.37\cdot {10}^{-4}$	$0.050$	$6.1\cdot {10}^{-6}$	$0.42\cdot {10}^{-4}$	$0.050$
AA-Na-ZZ-90	$3.1\cdot {10}^{-5}$	$1.2\cdot {10}^{-4}$	$0.034$	$1.9\cdot {10}^{-6}$	$0.61\cdot {10}^{-4}$	$0.090$

Equations

$$D_T^{\parallel }=\sum _{i=1}^N\frac{(x_i(\tau )-x_i(0){)}^2}{2N\tau },D_T^{\perp }=\sum _{i=1}^N\frac{(y_i(\tau )-y_i(0){)}^2}{2N\tau },$$

$$D_T^{\parallel (\perp )}=Z_{\parallel (\perp )}\exp \left(-\frac{E_{\parallel (\perp )}}{kT}\right),$$

$$F{T}_x(k_x,k_y)=\frac{1}{J}\sum _{j=0}^{J-1}\text{abs}\left[\sum _{n=0}^{N-1}\sum _{m=0}^{M-1}{x}_{mn}(j\cdot \Delta t)e^{2\pi i\left[k_{x}mL+k_{y}nW\right]}\right],$$

$$F{T}_x(\nu )=\frac{1}{J}\text{abs}\left[\sum _{j=0}^{J-1}x(j\cdot \Delta t)\cdot \exp (i2\pi \cdot \nu \cdot j\Delta t)\right],$$

Full text

Novosibirsk State University] Novosibirsk State University, Novosibirsk, 630090, Russia \alsoaffiliation[Nikolaev Institute] Nikolaev Institute of Inorganic Chemistry SB RAS, Novosibirsk, 630090, Russia

usc] Department of Physics and Astronomy, University of South Carolina, Columbia, SC 29208, USA

{tocentry}

Lithium-ion batteries 1, 2 have been continuously improved over the past few decades and they have now become the first choice for many consumer electronic devices, zero emission electric vehicles, and even sustainable energy technology 3, 4. A further progress in their capacity and cycling performance is directly related to the development of new electrolytes and advanced electrode materials. In particular, it has been shown that various wrapped, coated and crumpled composites 5, 6, 2 offer a superior performance compared to graphite for application as anode materials.

A fast diffusion of lithium ions through the anode material is one of the requirements to achieve high power, large capacity batteries that are able to operate at fast charge/discharge rates 7. In this paper, we demonstrate a spectacular enhancement of the interlayer diffusion of Li in the buckled bilayer graphene. Using molecular dynamics simulations, the anisotropic diffusion constants for the interlayer diffusion are found. A physical mechanism is put forward to explain the ultrafast diffusion that is observed in our simulations. Importantly, as many graphene composites contain curved graphene sheets on their surface, our findings are relevant to real experimental situations. Moreover, the anodes incorporating curved, buckled, or rolled graphene can easily be fabricated by suitable mechanical means.

There are many theoretical and experimental studies dedicated to the thermodynamic properties of lithium in graphitic materials. In particular, it was shown that lithium atoms have a fairly high adhesion energy $\sim 1.3$ eV8, 9 ($\sim 0.7-1$ eV for Na 10, 9), and that this energy is increasing in the vicinity of defects 8, 11. This leads to the directional asymmetry of the diffusion barrier. In pristine graphene, this barrier is $\sim 0.3$ eV for Li 12, 11, 13 and $\sim 0.13$ eV for Na 9. However, when an atom displaces towards a defect (such as a vacancy, divacancy or Stone-Wales defect), the barrier energy may decrease by up to $\sim 60\%$ for Li 8 and $\sim 80\%$ for Na 14, 15. A similar effect occurs when an atom moves towards an edge 10, 12.

Moreover, it has been shown that the diffusion barrier for alkaline atoms absorbed on a monolayer graphene depends on its curvature 9. Density functional theory (DFT) calculations have demonstrated that the energy barrier decreases/increases if the atom is absorbed from the concave/convex side. In a wide range, for any type of buckling, the energy barrier for Na atoms is about $0.1$ eV lower than that for Li atoms 9. Ref. 9 also reports a diffusion anisotropy with a variation of the energy barrier of $0.08$ eV for Li and $0.03$ eV for Na.

However, it is more difficult to make a similar analysis for the diffusion in the bilayer graphene because of a larger number of degrees of freedom. The effect of the interlayer distance on Li and Na intercalation and diffusion was evaluated in Refs. 16 and 13, respectively. Recently, Zhong et. al. 17 performed a series of ab initio calculations of lithium diffusion in AA and AB stacked graphene. According to these calculations, in the case of AA stacking, the diffusion is isotropic, and the energy barrier is $0.34$ eV. In the case of AB stacking, there are two distinct possibilities for lithium migration with the barrier heights of $0.07$ and $0.25$ eV. However, these previous results neglect the impact of the vibrational states on diffusion, which may be of significant importance 18.

Molecular dynamics (MD) simulations were carried out with NAMD2 software package 19 (NAMD was developed by the Theoretical and Computational Biophysics Group in the Beckman Institute for Advanced Science and Technology at the University of Illinois at Urbana-Champaign). To describe the interaction between carbon atoms, the CHARMM-like potential (including standard 2-body spring bond, 3-body angular bond, including the Urey-Bradley term, 4-body torsion angle and Lennard-Jones potential energy terms) that was previously optimized for graphene 20 was used. Knowing that the lithium-carbon interaction can be quite accurately described by the linear superposition of power-law functions of the type $1/r^{n}$ 21, the Lennard-Jones potential was employed to approximate Li-C and Na-C interactions. The Li-C bonding parameters were optimized to fit the equilibrium distance and energy from ab-initio calculations 21. A similar procedure was carried out to fit the ab-initio Na-C distance 22 and energy 9. The Van-der-Waals interactions were gradually cut off, starting at $10$Å from the atom until reaching zero $12$Å away.

MD simulations were performed with 1 fs time step. The Langevin dynamics with a damping parameter of $5\textnormal{ps}^{-1}$ was used for the temperature control. In all our simulations, the energy was first minimized in 2000 steps. Next, the system dynamics was simulated for $\tau=0.5$ns. Finally, the energy was minimized in 2000 steps.

Fig. 1 shows a typical simulation setup. Here, an initially optimized bilayer graphene sheet of $L=93.62$Å length (along $x$-axis) and $W=86.82$ Å width (along $y$-axis) is buckled around $x$-axis. In each simulation, a Li or Na atom was placed between the layers of AA- or AB-stacked bilayer graphene in a random initial position. The buckling was carried out by fixing the edge atoms at a distance $d<L$ (in case of fixing armchair edge (AC)) or $d<W$ (in case of fixing zig-zag edge (ZZ)). The role of the temperature was analyzed by performing simulations at three different temperatures (with a step of $50$K) for each graphene geometry. At each temperature, the statistics were collected in 250 runs. To correlate our results with data available in the literature, the alkaline diffusion coefficients were also found for a monolayer (1L) and bilayer (2L) graphene without any fixed atoms.

Anisotropic diffusion coefficients at the temperature $T$, $D_{T}^{\parallel}$ and $D_{T}^{\perp}$, were determined using the mean square displacement (MSD). Here, $\parallel$ is the direction about which the buckling takes place ($x$ in Fig. 1) and $\perp$ is the complementary direction ($y$ in Fig. 1). In the case of Fig. 1,

[TABLE]

where $x_{i}(\tau)$ and $y_{i}(\tau)$ are the coordinates of Li or Na atom at the final moment of time $\tau$ in the $i$-th run, $x_{i}(0)$ and $y_{i}(0)$ are the initial coordinates, and $N=250$ is the number of runs.

The temperature dependence of diffusion was fitted by the Arrhenius equation

[TABLE]

where $Z_{\parallel(\perp)}$ is the frequency factor 23 in the $\parallel$ (or $\perp$) direction, $E_{\parallel(\perp)}$ is the corresponding activation energy, and $k$ is the Boltzmann constant. In all simulated cases, MD results were fitted by Eq. 2 with a RMS error less then $1\%$.

To correlate our results with data available in the literature, the diffusion in flat graphene was simulated. The results of these calculations for Li are presented in the table 1, while for Na - in the Supporting Information (SI) table S.1 and Fig. S.1. Our general observations are that the sodium diffusion occurs more easily than lithium diffusion, and that the interlayer diffusion is more efficient than the surface diffusion. Fig. S.1 demonstrate an isotropic diffusion with MSD increasing with the temperature, which was exactly the expected behavior. Overall, the numerical diffusion parameters listed in tables 1 and S.1 are in good agreement with the experimental and numerical literature data. An order of magnitude deviation from the ultrafast experimental value $D_{uf}=7\cdot 10^{-5}$ cm2/s could be related to various factors in the experimental setup 25 not taken into account in our simple simulation. The authors of Ref. 25 also report lower values for $D$ (e.g., in the first lithiation cycle $D=5\cdot 10^{-6}$ cm2/s), and associate its further growth with the establishment of a Li intercalation pathway. Meanwhile, because the diffusion constant in buckled graphene can be of the same order in magnitude as $D_{uf}$, an (effective) buckling is a possible explanation for the observed ultrafast diffusion constant 25.

MD simulations for the diffusion in buckled structures were performed using the same procedure as for the flat ones. Some general observations are that the diffusion barrier for Na is slightly lower than that for lithium, the preferable diffusion direction is along the fixed edge (the $\parallel$ direction), and the greater the curvature of graphene the more pronounced is the diffusion anisotropy. For both Li (Fig. 2) and Na (SI Fig. S.2), the anisotropy of diffusion in AA-stacked graphene is stronger for the case of fixed zig-zag edges. For Li, a similar result was found for the case of AB stacking (see the last two lines in the table 2).

Our most striking finding is the spectacular increase in the $\parallel$ diffusion constant of Li atoms in buckled structures. Here, we refer to the results for AA-Li-ZZ-90 and AB-Li-ZZ-95 structures (table 2) showing about 30-fold increase in the diffusion compared to the flat single-layer graphene or 3-fold increase compared to the flat bilayer graphene (table 1). Contrastingly, a high diffusivity of Na in the flat bilayer graphene (SI table S.1) decreases with buckling, and is close to the “flat” value only in the AA-Na-ZZ-90 system (SI table S.2).

Naively, one can expect an increase in the frequency factor ($Z{{}_{\parallel(\perp)}}$) and diffusion energy barrier ($E_{{\parallel(\perp)}}$) with buckling (due to the higher rigidity of buckled structures). In reality, however, this expectation was not met. The numerical constants presented in tables 2 and S.2 show a mixed behavior. Finally, the energy barrier in the flat undistorted bilayer graphene (with all fixed carbon atoms) was identified by scanning the coordinate between two nearest local minima. It was found that $E=2.5$ eV for Li and $E=16$ eV for Na at the interlayer distance of $h=3.4$ Å. However, this barrier does not correspond to any real processes because in this calculation all of the carbon atoms were kept fixed.

It is well-known that the mechanical properties of graphene are anisotropic  26. The buckling by itself is another source of anisotropy. To reveal the asymmetry of lattice vibrations in the systems under consideration, Fig. 3 exhibits the time-averaged Fourier transform of carbon displacements

[TABLE]

where $J$ is the number of time steps, $M$ and $N$ define the grid size (each grid point corresponds to an atom), $x_{mn}(j\cdot\Delta t)$ is the atomic $x$-coordinate at the grid point $(m,n)$ at time $j\cdot\Delta t$, $\Delta t$ is the time step, $k_{x}=\tilde{m}/(LM)$, $k_{y}=\tilde{n}/(WN)$), $\tilde{m}$ and $\tilde{n}$ are integers, and $\text{abs}[...]$ denotes the absolute value. $FT_{y}$ and $FT_{z}$ were plotted using similar expressions.

Fig. 3 reveals high intensity swing modes at $\vec{k}\approx(0,0)$ in each graph. In the case of the free bilayer graphene, two symmetric lines of high amplitude oscillations (defined by $k_{x}=0$ and $k_{y}=0$) are also visible. In buckled structures, however, the vibrations in the $\perp$ direction are suppressed. In particular, Fig. 3 shows high-intensity vibrations of only about $k_{x}\approx 0$ when $x-$edges are in the $\parallel$ direction (AC case), and about $k_{y}\approx 0$ when the $y-$edges are in the $\parallel$ direction (ZZ case). This implies a strong correlation between the symmetry of high intensity vibrations and diffusion anisotropy. Given that the energy barrier calculated at all carbon atoms fixed is much higher than that from MD simulations at finite temperatures ($2.5$ eV versus 0.085 eV for Li at 300 K in free graphene), the physical displacement of carbon atoms is of crucial importance for the ultrafast diffusion discovered in our MD simulations.

In the buckled structures, while $k_{\perp}$ modes are suppressed, the oscillations with high-$k$ values ($|k_{x}|\approx 0.5$, $|k_{y}|\approx 0.4$) are enhanced. The wavelengths of these high-$k$ values modes are comparable to few lattice periods. Because such vibrations do strongly modify the order in the $x$ direction, and this effect is most noticeable for the case of fixed ZZ edges, one can associate the decrease in the energy barrier as a function of the curvature in ZZ systems with enhancement of such high-$k_{\parallel}$ vibrations.

Being lighter than Na, Li atoms must be affected more by lattice vibrations. To verify this statement, we calculate the frequency dependence of the oscillation density

[TABLE]

where $\nu_{n}=n/T_{0}$ is the frequency ($n$ is an integer), $x(j\cdot\Delta t)$ is the $x$-coordinate of an atom at time $j\cdot\Delta t$. In our work, the time-averaging was performed with the time step $\Delta t$. $FT_{y}$ and $FT_{z}$ were found using similar expressions. We note that the oscillation frequencies calculated based on Eq. (4) are in a good agreement with known DFT results 27.

Fig. 4 shows that the frequencies of lithium are higher than these of sodium. Interestingly, the most intense Li peak matches the peak of graphene. This frequency match could be associated either with the small mass of lithium or with a parametric resonance through which the lattice vibrations pump the lithium atoms. For some (not completely clear) reasons, the highest amplitude of Li oscillations is observed in $y$-direction, while the peak of oscillations along $x$-axis has a higher frequency. Whatever the reason (the resonance or small mass), the same prevailing frequencies of lithium and graphene oscillations provide in-phase oscillation condition that might contribute to the diffusion barrier decrease.

Our MD simulations have demonstrated a strong anisotropy of Li diffusion in the buckled bilayer graphene. Our main finding is the ultrafast diffusion of Li atoms in buckled structures with fixed zig-zag edges. In particular, the molecular dynamics simulations presented in this work have shown about 30-fold increase in the Li diffusion in buckled graphene compared to single-layer graphene (or about 3-fold increase compared to the free bilayer graphene). This effect increases with the degree of buckling and seems to be independent on the stacking type (although only AA- and AB-stacking were considered). Qualitatively, our findings could be explained by peculiarities of lattice vibrations in buckled structures possibly assisted by a parametric resonance.

Our work may facilitate the development of more efficient Li-ion batteries, and could be useful in various other applications and studies related to the diffusion of alkaline atoms.

{acknowledgement}

This work has been partially supported by the Russian Science Foundation grant No. 19-73-10068.

Supporting Information