Alkali metals inside bi-layer graphene and MoS2: Insights from first-principles calculations

Contrary to a wide-spread belief that alkali metal (AM) atoms intercalated into layered materials form singlelayer structures only, recent experiments [Nature 564 (2018) 234] showed that multi-layer configurations of lithium are possible in bi-layer graphene. Using state-of-the-art first-principles calculations, we systematically study the intercalation energetics for various AMs (Li, Na, K, Rb, Cs) in bi-layer graphene and MoS2. We demonstrate that for bi-layer graphene as host the formation energy of multi-layer structures is negative for K, Rb and Cs and only slightly positive for both Li and Na. In view of the previous experimental data on lithium, a multi-layer of Na might therefore form, while it is well-known that single-layers of Na in graphitic hosts are energetically very unfavorable. In MoS2, multi-layer structures are considerably higher in energy than the singlelayer ones, but the formation of the former can still occur, especially for the AMs with the lowest electronegativity. To rationalize the results, we assess the charge transfer from the intercalants to the host material and analyze the interplay between the ionic and covalent bonding of AM and host atoms. While our theoretical effort primarily focuses on the fundamental aspects of AM intercalation, our findings may stimulate experimental work addressing multi-layer intercalation to maximize the capacity of anode materials in AM ion batteries.


Introduction
A fossil fuel-free society is hardly possible without the development of light-weight but high-capacity rechargeable electrical power sources, such as alkali metal (AM)-ion batteries [1,2], which had extensively been studied since 1970s and entered the market in early 1990s. The significance of this scientific and technological breakthrough was reflected by the 2019 Nobel Prize in Chemistry awarded for the development of Li-ion batteries (LIBs). However, in spite of numerous applications of AM-ion batteries ranging from portable electronics to electric cars, they currently cannot match the growing demand for higher energy storage density, charging speed and cost reduction. At the same time, further improvements in their design and operation present a challenge, see Ref. [1,2] for an overview.
Specifically, modern LIBs have safety problems and narrow operating temperature ranges. Currently, graphite, which has a theoretical capacity of 372 mAhg À 1 , is normally used as the anode material [1], but other layered materials have been studied both experimentally and theoretically, including transition metal dichalcogenides (TMDs), MXenes, black phosphorus, both in bulk and nano-structured forms, see Ref. [3][4][5] for an overview. Particular attention has been paid to bi-layer graphene (BLG), as intercalation into this system is interesting not only in the context of AM storage [6][7][8][9][10], but also superconductivity [11,12]. Bi-layers can potentially provide more space for intercalants than bulk systems. Notably, a macroscopic three-dimensional bi-layer graphene foam with few defects was recently manufactured [6], and its Li-storage capacity and intercalation kinetics were systematically studied. It was concluded that in BLG Li atoms can be stored only between the graphene sheets, not on the outer surfaces. As for TMDs, and first of all, MoS 2 , the most common material in the TMD family, their layered structure and the weak van der Waals (vdW) interaction between the layers also enable the easy intercalation of Li ions without any significant increase in the volume and exhibit high Li storage capacities (up to 700 mAh/g) [13]. Understanding TMD behavior upon AM intercalation is important for controlling phase transitions from the semiconducting H to metallic T and T' phases as well [14][15][16][17][18][19][20]. Intercalation into heterostructures of these materials has also been investigated both experimentally [21] and theoretically [22,23].
All these studies on the intercalation of Li (or other AMs) into bulk and few-layer systems assumed that a single layer of the intercalant atoms is formed between the sheets of the host material. However, recent studies [24] employing in-situ transmission-electron microscopy (TEM), which is one of the most powerful techniques for getting insights into the dynamics of various processes in energy materials with atomic resolution [25], showed the unexpected formation of multi-layer close-packed Li phases between graphene sheets, hinting at the possibility of increasing the areal Li storage capacity. This gives rise to a question: is the formation of multi-layer structures between bi-layers of other 2D materials possible? Moreover, as sodium-ion [26] and potassium-ion batteries [27] have also been rapidly developed, another important issue is to understand the behavior of these (much cheaper than Li) and other AM atoms in bi-layers formed by graphene and TMDs.
Here, using density functional theory (DFT) calculations, we theoretically study the intercalation of AM atoms into bi-layer graphene and MoS 2 . We analyze the charge transfer and energetics of multi-layer configurations of various AMs between the layers of the host material and show that the formation of such structures should be possible with a few exceptions.

Computational methods
All our calculations were carried out using the VASP software package [28,29]. The exchange-correlation functional proposed by Bj€ orkman [30] was employed to account for the vdW interaction. An energy cut-off of 600 eV was used for the primitive cell calculations, for supercells we used 400 eV. The Brillouin zones of the primitive cells of the 2D materials and bulk AMs were sampled using 12x12x1 and 12x12x12 Monkhorst-Pack grid points [31]. The equivalent or better sampling was used for supercells consisting of up to 620 atoms (see Tables S1-S16). The maximum force on each atom was set to be less than 0.01 eV for the optimized configurations. The atomic structures and the charge densities were illustrated using the VESTA package [32]. Various initial configurations of AM atoms between graphene sheets (AA and AB stacking) and MoS 2 (H, T, T 0 -phases) were created, and then the geometry was fully optimized. In order to reduce the lattice mismatch for periodic structures, the optimum sizes of the supercells and rotational angles between surfaces were defined using the Virtual NanoLab software [33].
We assessed the stability of various AM structures between the sheets of graphene and MoS 2 by calculating their formation energy E f per AM atom defined as where E½Host þAM� is the energy of the supercell containing bi-layer graphene (or MoS 2 ) and AM atoms, E½Host� is the energy of pristine bilayer (graphene or MoS 2 ), n AM is the total number of AM atoms and μ AM is their chemical potential. Eq. (1) gives the energy per atom required to take AM atoms from the bulk structure and place them between the host material. Negative values indicate that this process is energetically favorable. Charge transfer from AM atoms to graphene (MoS 2 ) was evaluated as a difference between the integrated electron densities of the composite structure and isolated 2D material and the equivalent AM structures.

Results and discussion
Following the previous work [24] on multi-layer structures of Li atoms between graphene sheets, we studied first the intercalation of AMs into BLG, as sketched in Fig. 1 (a). We considered both the AB and AA stacking in BLG, as our calculations, in agreement with the previous results [22,34], showed that the AA stacking is only slightly (about 10 meV per atom) higher in energy, so that both configurations are generally possible. AM atoms with various densities were placed between the graphene sheets as single-or multi-layered structures, as schematically illustrated in Fig. 1(a), and then the geometry of the system was optimized. Examples of the optimized configurations are shown in Figs. S1, S2, S5, S6.
We found that the arrangements of AM atoms (M) in single-layers were the same as in the AMC 6 phase in the bulk host materials, which formed a commensurate ð ffi ffi ffi 3 p � ffi ffi ffi 3 p ÞR30 o superstructure (Fig. S1). We did not consider the AMC 8 phase, as it has a lower concentration of AM atoms, and overall our goal was to compare the energetics of multi-layer and single-layer AM structures. As for the geometry optimization of the multi-layered AM configurations, the initial atom positions were cut from the bulk fcc, hcp, bcc metals along the low-energy directions (e.g. [111] for the fcc structures) and placed between the graphene layers. The orientation of the metal crystal lattice with regard to graphene and sizes of the supercells were chosen to minimize the strain in the system, which normally did not exceed 1%.
The results of our calculations are presented in Fig. S3 and a summary is shown in Fig. 1(b). As the fcc, hcp, bcc structures of AMs in BLG have very close energies (for the same number of layers), Fig. S3, and because many bulk AM are bcc crystals, we presented here only the results for the bcc arrangement of atoms in the multi-layers. Note that E f for Li multi-layer configurations are slightly higher than those reported in Ref. [24], as they are for the bcc, not fcc/hcp arrangements of AM atoms.
It is evident that for all AMs, except for Na, E f is negative for a single layer of AMs in AA-stacked BLG. The trend is the same in AB-stacked BLG, but the values of E f are slightly higher due to the asymmetric positions of C atoms in the graphene planes with respect to AM atoms. Positive values of E f for Na in the C 6 NaC 6 ð ffi ffi ffi 3 p � ffi ffi ffi 3 p ÞR30 o phase and negative values for K have been also previously reported for graphite [1,9,35,36]. Formation of multi-layer AM structures is less energetically favorable than single-layers (except for Na): the energy increases with the number of layers, but remains negative. E f for two/three layers of Na atoms is positive, but it is comparable to that of Li, though, and as the latter have been experimentally observed [24], the multilayer Na structures may be experimentally realized in BLG. For K, Rb, and Cs, E f for multilayer structures increases with the number of layers, making the intercalation energetically less favorable, but theoretically possible, as the values are still negative. We stress that when the number of AM atom layers increases, E f should approach zero value, corresponding to an infinitely thick slab of AM atoms inserted between two graphene sheets. We note also that contrary to single-layer AM structures, multi-layers and graphene are incommensurate, as the spacing between AM atoms is governed by not only the interaction of AMs with the host material, but also between each other, as in the bulk crystal. Thus there are always AM atoms on top of different sites in graphene, but the 'net' effect is that the energies are lower for the A-A stacking in BLG.
For MoS 2 we considered both H and T/T 0 phases, as illustrated in Figs. S5 and S6. We found that the H phase had lower energy in most cases, as evident from Tables S7-S16 and Fig. S7. Although the T 0 -phase is expected to be energetically preferable at high concentrations of Li atoms, the concentration for single-layer structures we considered was always lower than the critical value (for Li, about 0.4 AM per formula unit [19,20]), while for multi-layer structures charge transfer into MoS 2 was in most cases not sufficient to make the T/T 0 phases energetically favorable. We also note that a direct comparison of the energetics of the H and T phases was not always straightforward due to different sizes and thus different number of atoms in the supercells. Because our main goal was to look at the differences in the behavior of the system upon formation of multi-layer AM structures, we present results for the H-phase only, Fig. 1(c).
The behavior of AMs between MoS 2 sheets is qualitatively the same for all AMs, except for Na: E f is smaller for a single layer of AMs than in multi-layers, and bonding becomes stronger from Na to Cs. Note a stronger bonding of AMs to MoS 2 as compared to graphene. For single layers of AMs in both graphene and MoS 2 bi-layers, our results are very close to those obtained for bulk graphite and MoS 2 serving as the hosts [36,37].
To get further insight into the energetics of the intercalants between the layers of the host material, E f can be represented [37] as where E coh is the cohesive energy of the bulk metal (energy required to split the bulk crystal into separate atoms, E coh > 0), E vdW is the loss in the vdW energy when the distance between the sheets of the host material is changed upon intercalation ðE vdW > 0Þ, and E b (E b > 0) is the binding energy of an AM atom to the host material (adsorption energy).
The contributions to E f are listed in Table S17. Similar to intercalation into bulk host materials, the behavior of E f for all AMs, except for Li, is defined by the difference between E coh and E b , as E vdW is much smaller. E b depends on the interplay between the charge transfer related to the electronegativity of AMs and covalent bonding between AM atoms and the host material. As shown previously for bulk hosts [36,37], the covalent bonding between BLG and Li is very strong affecting charge transfer [38] and the ionic contribution.
To assess charge transfer in the bi-layer systems upon intercalation and to compare it to that in the corresponding bulk host materials, we calculated the difference Δρ between electron densities of the combined and separated systems as where ρ½Host þAM� is the total electron density of the host material and embedded AM atoms, and ρ½Host� and ρ½AM� are the electron densities of the isolated host structures and AM layers for the same positions of the atoms, respectively.
For BLG and a single layer of AM atoms, a cross-section of Δρ is presented in the top panel of Fig. 2, and the corresponding values for MoS 2 bi-layer are shown in the bottom panels. It is evident that strong covalent bonds between C and Li atoms have been formed, and although the electron density is depleted between the AMs, a substantial part of it is still localized near the atoms in the covalent bonds. For MoS 2 bi-layer as a host, the situation is different: the covalent interaction is roughly the same for all the AMs. It is also evident that charge transfer from AMs to the host materials should increase from Na to Cs in case of MoS 2 , in agreement with the trend in the electro-negativity of these metals. However, for graphene the situation is more complicated, as discussed below.
We also visualized Δρ for double and triple layers of AMs between graphene and MoS 2 sheets. The results for double layers of AMs are presented in Fig. 3 and for triple layers in Fig. S8. It is evident that the interaction between Li and the host material for multi-layers is also rather strong, it even gives rise to a redistribution of the electron density between Li atoms, making the bond between them more 'covalent'. The overall charge transfer from AMs to the host material increases from Na to Cs as well, as expected, similar to the case of a single layer of AMs. It is also clear that the AM atoms facing the host material mostly donate the electrons, while the inner atoms (e.g., in the triple layer) preserve their charges.
Quantifying the charge transfer in these systems is not straightforward, and the results may depend on the methods used to evaluate it, as discussed at length previously [38]. The Bader analysis [39,40] based on finding the extrema in the electron density and splitting the space to atomic volumes accordingly is normally used. However, it accounts for the geometrical charge transfer which occurs in every system upon adding atoms, even when there is no physical charge transfer between the host and the added atoms [38].
Thus, we evaluated charge transfer using two approaches: using the Bader method and by averaging Δρ within the planes parallel to the sheets of the host material, and plotting it as a function of z-coordinate (perpendicular to the planes) as schematically shown in the inset in Fig. 4(a). Charge difference for a single layer of AMs between graphene sheets is shown in Fig. 4(a). It is evident that the averaged electron density is depleted within the AM plane (compare it to the cross-section of Δρ shown in Fig. 2), and is increased considerably in the area between graphene sheets and the AM metals, especially for Li and Na. Note that the graphene planes are at different locations due to different atomic radii of the AM atoms. Then we evaluated the total charge transfer Δq z by integrating over the range of z where Δρ is negative, as illustrated for Li in Fig. 4(a). We stress that the choice of the integration range is not unique, but if the same criterion is used for all AMs, this approach makes it possible to analyze the trends.
The charge transfer Δq z for single-and multi-layers of AMs between graphene sheets is shown in Fig. 4(b). We also present similar results (for the AMC 8 phase) for graphite taken from Ref. [38]. For a single layer of Li and Na atoms, when the covalent interaction is strong and atomic radii are small, the results for bi-layer graphene and the bulk system are very close, while the difference for AMs with larger atomic radii can be explained by the different geometry: there is no layer of AMs on the other side of the graphene plane and the system has more space to adjust its geometry. For the double layer of AM atoms, charge transfer increases from Li to Cs, and the same is true for triple layers.
The results obtained using the Bader method, Fig. 4(c), give qualitatively (and even quantitatively for all AMs, except for Li) the same picture for double and triple layers of AM atoms. However, for single layers our findings are somewhat counterintuitive, although the results for Li in BLG (0.84e per Li atom) agree well with the previous [41] Bader results (0.88e). They indicate that charge transfer decreases from Li to Rb, contrary to what one can expect from lower electronegativity of heavier AMs. The reason for that is the strong covalent bonding between graphene and light AM atoms, which gives rise to a redistribution of the electron density and its build-up in the covalent bonds between carbon and AM atoms. This is evident from 4(a): charge build-up is clearly observable, and even a part of the electron density from graphene sheets went into the new bonds. Charge transfer (per AM atom) decreases with the number of the AM layers, as only the top layers are depleted of electrons [7] and the inner layers mostly retain their metallic character.
In order to get further insight into the nature of chemical bonding between graphene and AMs, we also calculated the electron localization function (ELF) [42,43] for single layers of AMs. The ELF makes it possible to assess the distribution of the electron density in chemical bonds. The ELF for Li in graphene is shown in Fig. S9. It is evident that the picture is dominated by strong carbon-carbon bonds in graphene. Thus, in order to improve visualisation, we subtracted the values of the ELF obtained for isolated graphene sheets, following the approach previously used for other systems with different types of bonding [44]. The results are presented in Fig. S10. It is clear that the bonds between AM and carbon atoms, when moving from Li to Cs, become less and less covalent, confirming the picture obtained from the analysis of the electron density.
Having analyzed charge transfer in BLG, we moved on to MoS 2 . Similar to graphene, we carried out averaging in the planes parallel to the sheets. The charge difference for a single layer of AMs between MoS 2 sheets is shown in Fig. 4(d). The averaged electron density is also depleted close to the AM plane, but its build-up in the covalent bonds with regard to the depletion area is smaller and rather uniform for all AMs. As a result of this, the total charge transfer Δq z integrated over the range of z where Δρ is negative, increases from Li to Cs, Fig. 4 (e). The same is true for AM double and triple layers, and the charge transfer per AM atom is smaller. The results obtained using the Bader analysis in Fig. 4 (f) also indicate that Δq increases from Li to Cs, and gets smaller when the number of AM layers increases.
Using the calculated geometries for multi-layers of AMs in graphene and MoS 2 bi-layers, we estimated the theoretically achievable capacity of such systems, which is presented in Table S18. For triple-layered Li and Na structures in BLG it proved to be equal to 828 mAhg À 1 and 440 mAhg À 1 respectively. This indicates that the capacity of BLG with intercalated multi-layered Na structures can exceed that for graphite with lithium (372 mAhg À 1 ). In the case of lithium between doublelayered MoS 2 the capacity reaches 251 mAhg À 1 for three layers of Li and could further be increased by adding more Li layers. The capacity of three-layered sodium in bi-layered-MoS 2 (157 mAhg À 1 ) is rather low, though. It should be pointed out that in general the capacity can be increased by making vertical vdW graphene-MoS 2 heterostructures [23]. Such heterostructures are expected to be beneficial for multi-layer AM storage due to the weight reduction of the host by using the lighter graphene and decrease of the formation energy due to the dichalcogenide layer. There have been only a few reports (see Ref. [45] for an overview) on vdW heterostructures used for energy storage due to many challenges in their fabrication processes.

Conclusions
By using first-principles calculations we have studied the energetics of various AMs (Li, Na, K, Rb, Cs) in bi-layer graphene and MoS 2 . We have demonstrated that the formation energy of multi-layer AM structures between graphene sheets is negative for K, Rb and Cs. Both Li and Na multi-layer structures have positive formation energy, but the values are small. In MoS 2 , AM multi-layer structures are considerably higher in energy than the single-layer ones. Multilayer structures can still occur, as intercalation is energetically favorable, especially for the AMs with the lowest electro-negativity. To rationalize the results, we assessed the charge transfer from the intercalants to the host material and analyzed the interplay between the ionic and covalent bonding of alkali metal and host atoms. Our results, in conjunction with the recent experimental discovery of multi-layer lithium in bi-layer graphene [24], raise hope that even for Na a multi-layer structure in bilayer graphene can form, although the single-layer configuration is energetically strongly unfavorable for this metal. This result is intriguing in the context of the recent progress in manufacturing macroscopic amounts of bi-layer graphene [6] and the low cost/abundance of Na, key prerequisites to envisage the implementation of Na-ion batteries based on nano-structured 2D materials. Thus, while our theoretical effort primarily focuses on the fundamental aspects of AM intercalation, our findings may stimulate experimental work addressing multi-layer intercalation to maximize the capacity of anode materials in AM ion batteries.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Dr. Zakhar I. Popov received his PhD in Physics from Kirensky Institute of Physics, Siberian Branch of the Russian Academy of Sciences, Krasnoyarsk, Russia in 2013. Currently, he is a senior researcher at the Emanuel Institute of Biochemical Physics, Russian Academy of Sciences, Moscow, Russia and a visiting researcher at NUST "MISiS", Moscow, Russia. His scientific interests lie in the areas of computational materials science, electronic structure calculations, two-dimensional materials, chemical reactions, and catalysis. Throughout his career, he has co-authored more than 60 peer-refereed papers.
Dr. Jurgen Smet received his PhD in Electrical Engineering and Computer Science from MIT in 1994. Currently he is the head of the Solid State Nanophysics Group at the Max Planck Institute for Solid State Research in Stuttgart, Germany. His scientific interests include Coulomb interaction, correlation and spin phenomena in two dimensional electron systems hosted by conventional semiconductor and van der Waals heterostructures as well as ionic transport in 2D materials. He has co-authored about 200 peer-refereed papers. Awards include the Walter Schottky Prize (DPG) and the Gerhard Hess (DFG) and Nano-Futur (BMBF) young investigator awards.
Dr. Arkady V. Krasheninnikov received his PhD in Physics from Moscow State Engineering Physics Institute in 1995. Currently he is a Group Leader at the Institute of Ion Beam Physics and Materials Research, Helmholtz-Zentrum Dresden-Rossendorf, Germany, and a Visiting Professor at Aalto University, Finland. His scientific interests lie in the areas of computational materials science, electronic structure calculations, twodimensional materials, and irradiation effects in solids. Throughout his career, he has co-authored more than 200 peerrefereed papers. His awards and recognitions include HZDR Research Award, Highly Cited Researcher (Physics) from Clarivate Analytics (ex. Web of Science); American Physical Society "Outstanding Referee".