Sphere-in-Contact Model/Theorem for the Development of New Graphite-Based Intercalated Battery Materials

1. Introduction

Li-ion batteries are powering smartphones, laptops, electric vehicles and other electronic devices are now broadly used in many commercial applications². Their invention is credited to John B. Goodenough, M. Stanley Whittingham and Akira Yoshino in 2019 with the Noble price in chemistry “for the development of lithium-ion batteries”^3-5. The battery anode in Li-ion batteries is composed of graphite which stores the Li⁺ cations during the charging cycle of the battery⁶. The Li⁺ diffuse between the basal planes of graphite forming structures that resemble graphite intercalation compounds (GICs)⁷. The exact mechanism for the diffusion and storage of Li⁺ cations in GICs is not well understood and various proposals for this have been made⁸. In battery materials the rate of charge-discharge cycles is a critical component in the design of battery materials and setups that require less time to be charged or discharged are being developed⁹. This rate is associated to the rate of diffusion of Li⁺ cations between the graphene layers of graphite². Additionally, studies have since investigated the electrochemical properties of graphite and alternative anode materials to improve rate capabilities, thermal stability, and overall energy density; a pursuit driven by the ever-rising demand for longer-lasting, higher-capacity batteries in consumer electronics and electric vehicles^{10, 11}. We therefore in this study have used sphere-in-contact models of GICs that have an AαAα stacking sequence, to study the diffusion of Li⁺ cations for GIC having various empirical formulas. These empirical formulas where given by the placement of intercalate atoms on the basal plane of graphene in a configuration that minimizes the Coulombic repulsive energy between cations. These in general are intercalate structures that are hexagonally symmetric and place the intercalate atoms at their largest possible distances, as we have recently shown is also the arrangement that monoatomic adsorbates attain when adsorbed on close-packed surfaces of metals¹².

The sphere-in-contact model/theorem is a molecular and materials model where the atoms have their atomic radius (e.g. ionic radius, metallic and covalent). This results in spheres that are in contact for structures that have either covalent, ionic or metallic bonding and these structures show in an improved way the electron cloud volume of an atom, an ion or even a metallic atom. We have previously shown that this model can better describe the structure of carbon materials (e.g. graphene, graphite, carbon nanotube, fullerene) as it correctly shows the void space in the structure. ¹³ We have also used sphere-in-contact models to study the position of close packed nanoparticles on metal surfaces ¹⁴ and prepared miniature molecular models of heterogeneous catalysts, supported metal nanoparticles and organic molecules. ¹⁵ We have also used it in the elucidation of the cap structure of (3,3)-, (4,4)- and (5,5)-SWCNT where physical sphere-in-contact models were constructed to understand the geometry of the spherical cap. ¹⁶ Also as electromagnetic radiation interacts with the electron density cloud of atoms in sphere-in-contact models was applied to model the structure of an X-ray filter in which graphite had the rhombohedral structure. ¹⁷ Recently we have used such models to study diffusion mechanisms of adsorbates on close-packed metallic surfaces.¹² So the use of the sphere-in-contact model is useful when analytical solutions of the distances between charged surfaces of atoms is desirable. These inter-surface distances can be used to evaluate the repulsive interaction energy between adsorbates and dopants in surface and materials chemistry using the Coulomb’s law, where the partial charges are not centred at the nucleus of the atoms, but rather on the surface of a sphere that has the atomic radius of the atom. The evaluation of Coulomb terms is computationally inexpensive and a method that can show clearly why certain dopant geometries or adsorbate structures that are symmetric are preferred as they reduce the interatomic repulsive energy term. It also useful as an educational tool as it directly translates concepts that students are taught in secondary school, such as basic laws of physics (i.e. Coulomb’s law), mathematics (e.g. trigonometry, geometry) and the existence of partial charges in polar covalent bonds to structures of doped materials and ordering of adsorbates on surfaces of metals and 2D layered materials (e.g. graphite intercalation compounds).

1.1. Computational Studies of GICs as Battery Materials

Computational methods have been applied extensively in order to investigate the structural, energetic, and transport properties of graphite intercalation compounds (GICs), offering insight into lithium staging behaviour in layered materials like graphite, interlayer interactions, and the dynamics of ion diffusion. Density functional theory (DFT) calculations that account for van der Waals interactions have demonstrated that the stability of lithium-intercalated graphite structures, such as LiC₆ and LiC₁₂ configurations, sensitively depends on intercalant type, charge density, and interlayer spacing^{18, 19}. These studies show that ionic size and the charge distribution of alkali metals influence both the spacing between graphite layers and the staging structure, which impact lithium-ion (Li⁺) mobility together. Analysis of diffusion barriers has shown that lithium tends to move more easily along flat planes of graphite rather than between them, following paths that minimize repulsive interactions consistent with geometric models²⁰. First-principles modelling of staged compounds has further helped to clarify how structural changes occur with shifts in ion concentration, emphasizing the energetic trade-offs between maintaining stage stability and transport efficiency²¹. Molecular-scale simulations of graphite-electrolyte interfaces also suggest that solvent molecules and interfacial structure influences how easily ions move into the material and how likely the graphite layers respond and separate, with size and polarity effects closely linked to the spacing and arrangement of the graphite layers^{22, 23}. These findings support the use of models that are more spatially focused such as the sphere-in-contact model/theorem, which helps visualize the space between atoms and offers an intuitive way to understand how 3D-structure affects ion movement in such layered carbon materials.

1.2. Empirical Formula of Various GIC’s Using the Sphere-in-Contact Theorem

We have explored various possibilities for GIC’s using the sphere-in-contact theorem to models these materials. These models were designed based on the assumption that the intercalates atoms will have some partial positive charge on their atoms, due to charge transfer to the graphene layers. Therefore these structures have configurations that maximise the distances of the intercalates, in order to reduce the repulsive interactions, and therefore enhance the stability of the GIC. The partial charge sign in this treatment of repulsive intercalate interactions does not effect the energetic trends found for the various GIC as a function of intercalate concentration. With the use of this approach we find some intercalate geometries that are depicted in Figure 1, which shows the position of the intercalate (red atom) on the basal plane of graphene. The second graphene layer is not shown for better clarity but we assume an AαAα stacking sequence of the graphene layers, which is known to be the stacking sequence in most GICs.¹

Having done this exploration of possible GICs we find that the empirical formula possible are MC₂, MC₆, M₃C₈, MC₃ and MC₈ although other less dense GIC are also possible. These would not be of considerable interest as the charge carrier density in such hypothetical GICs would not be high and therefore their application in ion-battery materials would not be of practical use, as such materials would not have large charge capacities. Therefore we consider here only the aforementioned empirical formulas that are shown in Figure 1.

1.3. Definition of Unit Cell in GICs

In this section we define the unit cell of the various GICs shown in Figure 1. On the basis of simple trigonometric and geometric rules we can calculate the unit cell parameters in terms of three angles (i.e. α, β, γ) and three unit cell vectors (i.e. a, b, c). The angles of the lateral unit cell are depicted in Figure 2 along with the various empirical formulas found.

In the calculation of the empirical formula we make use of the following basic geometric rules:

a)

An atom at the corner of the unit cell will have 1/8 of its mass within the unit cell if the unit cell angles are all 90°.

b)

An atom at the corner of the unit cell which has an angle φ and the other two 90°, will have φ / 720° of its mass within the unit cell.

c)

An atom at the edge of the unit cell will have ¼ of its mass within the unit cell is the unit cell angles are all 90°.

d)

An atom at the edge of the unit cell where one of the angle is φ and the other two 90°, will have φ / 360° of its mass within the unit cell.

e)

An atom at the side of the unit cell will have ½ of its mass within the unit cell.

f)

An atom within the unit cell will have an equivalent of one atom in the unit cell.

Using these rules we derive the various empirical formulas shown in Figure 2 making the assumption that the intercalates will always reside in the center of the hexagonal ring on the graphene layers, due to the higher coordination number at these positions, which reduces the charge on the intercalate and therefore the GIC would be more stable. Also the configuration of intercalates is such that there is a maximum average distance between them, which would minimise the repulsive energy between intercalates and therefore would be a lower energy configuration. This would also presumably be a preferred configuration when the ion-battery material is not in charged state.

1.4. Calculation of the Density of Various GIC’s Within the Sphere-in-Contact Theorem

The density of an arbitrary GIC using the sphere-in-contact theorem can be defined as

$$Density={\displaystyle \frac{Mass}{Volume}}={\displaystyle \frac{{(N}_{Li}\cdot m_{Li})+{(N}_C\cdot m_C)}{V_{unitcell}}}$$

(1)

The volume of any unit cell in a material can be calculated from the volume of the triclinic unit system which is given from

$$V_{unitcell}=a·b·c·\sqrt{(1-{cos}^2\alpha -{cos}^2\beta -{cos}^2\gamma +2·cos\alpha ·cos\beta ·cos\gamma }$$

(2)

where $a$, $b$ and $c$ are the lattice vectors and $\alpha ,$ $\beta$ and $\gamma$ the angles between the lattice vectors as shown in Figure 3.

Here the primitive unit cell has rhombus base and height c. If we start from the volume of the triclinic unit cell the unit cell volume is given by,

$$V_{unitcell}=a^2·c·sin\gamma$$

(3)

where $a=2·\sqrt{3}·r_C$and $\gamma ={60}^{°}$ and $c$ is the height of the unit cell of the GIC, which is given by,

$$c=2·\sqrt{{\left(r_C+r_{Li}\right)}^2-4·{r_C}^2}$$

(4)

The trigonometric and geometric derivation of the height of the unit cell is given in section 1.5.

The density of the GIC is given by combining Eqns. 1, 3 and 4 and is given by

$$\rho \left(GIC\right)={\displaystyle \frac{{(N}_{Li}\cdot m_{Li})+{(N}_C\cdot m_C)}{{2·a}^2\cdot sin\gamma \sqrt{{\left(r_C+r_{Li}\right)}^2-4·{r_C}^2}}}$$

(5)

where $N_{Li}$ and $N_C$ is the number of lithium and carbon atoms within the unit cell of the GIC, $m_{Li}$ and $m_C$ are the atomic masses of Li and C, respectively, and $r_{Li}$ and $r_C$ are the atomic radii of the lithium ion and the covalent radius of carbon in graphite, respectively. Noteworthy, is that this equation is not limited to Li-ion GICs but can be used for any momoatomic (e.g. Li, Na, K, Rb, Cs) intercalates between graphite layers.

We have used the cationic radius of lithium to calculate the density of such GIC may depend on the charged state of the GIC. We expect that for the intercalate atoms to be mobile in the presence of an external electric field, that is parallel the graphite basal plane that these are cationic lithium species, mainly due to their smaller radius and charge compared to the atomic radius of a neutral lithium intercalates, are the mobile species. It is also conceivable that not all intercalates are in the same charged state during the charging/discharging of such materials. So the radius of the intercalate can be between the atomic radius of lithium and the ionic radius of lithium and that the charged state of the intercalate atoms mainly affects the interlayer distance, and that the carbon hexagonal framework is affected to a smaller extent, as the interatomic distances in that are mainly determined by the sigma covalent framework, that does not participate in the bonding with the intercalates. Based on MP2 calculations of the lithium cation between two benzene rings we observe little variations of the C-C bond lengths, which is not affected by the formation of the ion-π bond between lithium cations and the conjugated π-clouds in cyclic carbon rings.

1.5. Height of Unit Cell Using the Sphere-in-Contact Model/Theorem for Lithium GICs

Figure 4 shows the top view and side view of an arbitrary GIC with AαAα stacking sequence, where the intercalates are in the center of the hexagonal ring of graphite and in which the height of the unit cell is shown (c). With the use of Pythagoras’ theorem and the dimensions of the triangle shown in Figure 4, given as a function of the radius of the carbon atom and the radius of a lithium cation, we are able to calculate the height of the unit cell which is given by c. We only consider GIC models where the interlayer spacing between the graphene layers is that of LiC₆, which has been found to be a stable form of lithium graphite intercalation compound under ambient conditions of pressure and temperature. We expect this assumption to be valid if the intercalation of lithium atoms is maximum when there is one lithium atom per hexagonal ring of graphite and every layer of graphite is fully saturated with lithium atoms, but is has been previously confirmed via in situ neutron diffraction that the lattice dilation during charging of lithium ion batteries is only 4%²⁴, which corresponded to an empirical formula of LiC₆. However, in this study we explore some higher in concentration Li⁺ phases in GICs as we expect based one purely geometric considerations that such higher in density lithium phases in graphite would also result in lattice dilation of the order of 4%. This is based on a number of experimental and computational studies that suggest that such superdense phases of lithium are possible to an empirical formula of LiC₂, which is the highest density phase considered in this study. For example, XPS studies with the use of piston-cylinder high-pressure chamber indicate that lithium intercalation to LiC₂ causes insignificant changes in the chemical shifts of carbon atoms, consistent with the assumption made here that there is not a significant lattice dilation comparing LiC₆ an LiC₂ structures.²⁵ Lithium overintercalation, beyond LiC₆, has also been observed recently at ambient pressure, using static ⁷Li nuclear magnetic resonance (NMR).²⁶ A very recent study that used ball milling to prepare LiC₃ under ambient conditions found that this structure is also stable at ambient conditions.²⁷ Noteworthy, the stability of the LiC₃ phase was previously predicted on the basis of Density Functional Theory (DFT) calculations, in addition to the existence of a LiC₂ phase which was found metastable.²⁸

So the assumption of one Li⁺ per carbon ring is the maximum Li ion concentration that can be evaluated using the current approach within the sphere-in-contact theorem, and that the interlayer distance is the same for highly dense Li-ion GICs, as dense as LiC₂ would result in a possible error of less than 4%, in the volume calculation of he unit cell, of these arbitrary GICs. We note here that although we explore these high density phases of lithium in graphite, such as LiC₃ and LiC₂, such high density phases could result in considerable clustering between lithium ions into metallic lithium clusters, under reducing conditions. Such metallic clusters for lithium have been observed via infrared spectroscopy of LiC₄ and LiC₂ ²⁹ and are potential pathways for the chemical and mechanical degradation of such lithium ion GICs, in ion-battery materials, under reducing conditions (i.e. excess of electrons within these materials).

1.6. Interlayer Distance in LiC₆ and Radius of Li atoms Intercalated

There have been many crystallographic studies to accurately estimate the in interlayer spacing in LiC₆ and the stacking sequence. Crystallography suggests that the stacking sequence is AαAα with an interlayer distance of 3.70 Å, based on the stable form of this GIC at ambient temperature and pressure, in samples that where prepared by annealing or under pressure, causing lithium to intercalate between the layers.¹ This interlayer distance is also in agreement with a study that investigated it as a function of lithium intercalation obtaining the phase diagram of Li_xC₆., where x was varied between 0 and 1.³⁰ Additionally, LiC₆ prepared by electrochemical reduction of natural graphite in a LiClO₄ solution at 30 °C found a interlayer distance of 3.70 Å.³¹ This interlayer distance of 3.70 Å as well as the AαAα stacking sequence was also confirmed by neutron diffraction at various temperatures.³² This interlayer distance in comparison to that of hexagonal ABAB stacked graphite (3.354 Å)³³ is 10.3% increased due to Li intercalation. We have therefore used this distance to estimate the radius of lithium ions within the sphere-in-contact model, when the carbon atoms have their covalent radius that corresponds to half a bond length of the C-C bond in graphite, 1.422 Å.³³ With the use of Figure 4 we therefore calculated the radius of the lithium ions in the sphere-in-contact model to be 1.62 Å and that the covalent radius of carbon is 0.711 Å.

So in the next section we will perform calculations of the unit cell volume and the density in GICs using Eqn. 5 and the radius of intercalated lithium atoms, $r_{Li}=1.62\times 10^{-8}cm$ and the covalent radius for the carbon $r_C=0.711\times 10^{-8}cm$.

1.7. Density and Volume of the Unit Cell of GIC

So the volume of the unit cell of LiC₆ is,

$$V_{{LiC}_6}=36\sqrt{3}{r_C}^2\sqrt{{\left(r_C+r_{Li}\right)}^2-4·{r_C}^2}=5.822·{10}^{-23}{cm}^3$$

(6)

The GIC with the empirical formula LiC₆ has one Li atom and six carbon atom per unit cell, so the density is given by

$$\rho ={\displaystyle \frac{m_{Li}+{6·m}_C}{V_{{LiC}_6}}}$$

(7)

Where $m_{Li}$ and $m_C$ are the mass of the lithium and carbon atoms which were calculated from the atomic masses of these atoms as taken from the IUPAC periodic table and are $m_{Li}=1.1526·{10}^{-23}g$ and $m_C=1.9944·{10}^{-23}g$. So the density of LiC₆ is 2.25 g cm^-3 which is in excellent agreement with the literature value of 2.21 g cm^-3 based on XRD.¹. This suggests that Li in LiC₆ is not in a +1 charged state under ambient conditions and that the radius of Li in LiC₆ is closer to the atomic (i.e. 12 coordinated metallic radius) radius of Li (1.5 Å)³⁴ , rather than the ionic radius of Li⁺ (r_CN4 = 0.59 Å, r_CN6 = 0.76 Å, r_CN8 = 0.92 Å)³⁵ estimated from empirical analysis of interatomic distances in lithium oxides. This is an interesting result that is mainly associated with the property of lithium having a significantly larger radius than carbon which is associated to the smaller effective nuclear charge of lithium (Z_eff) compared to carbon. Furthermore, it suggests that based on the atom radius of Li in the sphere-in-contact model that the most dense phase of lithium in graphite is LiC₆ and that any further addition of Li to graphite, forming LiC₂ and LiC₃, which have been suggested as superdense phases of lithium atoms in GICs would result in a significant overlap of the electron densities of the lithium atoms, which would suggest a rather metallic 2D array of lithium atoms, rather than isolated lithium atoms. This is perhaps the reason that these superdense configurations of lithium atoms in GIC are not observed under ambient conditions and require mechanical compression or other extreme electrochemical conditions.

In Table 1 we present the calculated volume of the unit cell in the various empirical formulas studied and Table 2 we present the density along with analytical equations to calculate this based on the sphere-in-contact model of LiC₆.

1.8. Packing Fraction of GIC and Diffusion Rate of Intercalates

From the packing fraction the available space for intercalate diffusion can be estimated. The void space fraction of the atoms in a unit cell is given as 1 minus the packing fraction. Therefore, the void space in a material correlates to the space available for atom/ion (i.e. intercalate) diffusion. This correlation can be made based on the fact that the lowest barrier for intercalate diffusion within the graphene sheets will be roughly the same for any density of intercalates between the layers as intercalates are hold by weaker interactions (e.g. ion-π interactions), compared to the stronger covalent C-C bonds within the graphene layers. This barrier is expected to be lower over the C-C bonds of the graphene layers, rather than over the C positions, as the coordination is greater when the adsorbate is directly on top of the covalent bonds, connecting the carbon atoms, which suggests a lower barrier for diffusion from the center of one hexagonal carbon ring to the centre of an adjacent hexagonal ring (see Figure 2). So the parameter that limits the diffusion in the case of GICs is not the barrier for diffusion but the available void space between the layers and the path that the ions follow during this diffusion, which also determines the capacity of the batteries. This is based on physical understanding of this material which suggest that the barrier of diffusion is very high, if there is already and intercalate in the adjacent hexagonal ring. We therefore suggest that the packing fraction using the sphere-in-contact model/theorem can used in materials design.

Although there are not any reported crystallographic structures for LiC₁₀ there have been reports of GICs with empirical formula MC₈. KC₈, RbC₈, and CsC₈ are known examples of GIC with this empirical formula³⁶, which is a less dense phase of the intercalated atoms in graphite. This is perhaps associated with the atom radius of the intercalate as potassium (r_K = 2.3 Å)³⁴, rubidium (r_Rb = 2.5 Å)³⁴ and caesium (r_Cs = 2.7 Å)³⁴ which all have significantly larger radii than Li (r_Li = 1.5 Å)³⁴ and therefore are more stable at a larger separation, to prevent their transition to a 2D metallic state. This is in agreement with the interatomic distance of intercalates in the MC₈ empirical formula in which the interatomic distance between intercalates is about 4.926 Å, and consequently the nearest surface-to-surface distance between intercalates is larger compared to the interatomic distance between intercalates in the MC₆ structure which is about 4.266 Å. Therefore, the MC₈ structure in the sphere-in-contact model will have smaller repulsion between the intercalated atoms and therefore, bulkier atoms can be accommodated between the layers, at relatively larger separations. The effect of the radius of these atoms on the interlayer distance in MC₈ (d_K = 5.35 Å)³⁷, rubidium (d_Rb = 5.65 Å)³⁷ and cesium (d_Cs = 5.94 Å) has been previously reported and is consistent with the trend of their atomic radius.³⁷

Additionally, It is possible that a MC’₈ structure can exist, as based on the interatomic distances between intercalates it would have intermediate stability to that of LiC₆ and LiC₈, but in this structure the interatomic distances in one direction are different than in the perpendicular direction, so 2D anisotropy would be expected and this structure would be more stable if it relaxes to MC₈, based on the repulsive interactions between intercalates, as presented in section 1.9.

From the packing fraction the available space for intercalate diffusion can be estimated. The packing fraction is the % space occupied by the atoms in the unit cell of the material. Therefore, the packing fraction is inversely proportional to the available empty space (i.e. void space) in a material, which correlates to the space available for atom/ion (i.e. intercalate) diffusion. This correlation can be made based on evidence that the lowest barrier for intercalate diffusion within the graphene sheets will be roughly the same for any density of intercalates between the layers as intercalates are hold by weaker interactions (e.g. ion-π interactions), compared to the stronger covalent C-C bonds within the graphene layers. This barrier is expected to be lower over the C-C bonds of the graphene layer rather than over the C positions, as the coordination is greater when the adsorbate is directly on top of the covalent bonds, connecting the carbon atoms, which suggests a lower barrier for diffusion from the center of one hexagonal carbon ring to the centre of an adjacent hexagonal ring. The barrier for this diffusion has been previously estimated based on DFT calculations and it was found that for larger interlayer spacing this barrier reduces, increasing therefore the diffusion rate by 3-orders of magnitude, when the interlayer distance increases between 3.4 to 3.8 Å.³⁸ In addition to this, it is conceivable that when the intercalates are in close proximity, that the barrier for ion diffusion will become exponentially high due to overlap of their electron density clouds, increasing the repulsive interactions, considerably. So this suggests that a more important factor for the diffusion of intercalates between graphite layers is the available void space, that is inversely proportional to the atomic packing fraction ($\eta )$ given by,

$$\eta ={\displaystyle \frac{Volumeoccupiedbyatoms}{Volumeofunitcell}}$$

(8)

So we can approximate the diffusion rate parameter ($D$) is being proportional to the void space given by,

$$D\propto 1-\eta$$

(9)

Apart from this purely geometric factor the relative charge transfer rate is directly proportional to the molar intercalate density $C$ (units: mol⋅cm^-3) of the GIC which is given by the number of intercalate atoms in the unit cell divided by the volume of the unit cell and Avogadro’s number,

$$D\propto C={\displaystyle \frac{numberofintercalateatomsinunitcell}{N_{Av}{V}_{unitcell}}}$$

(10)

$$D\propto C·\left(1-\eta \right)$$

(11)

As the specific atomic radius of the intercalate atoms, for certain empirical formulas, results in considerable overlap of the electron cloud of the intercalate atoms, the void space in the sphere-in-contact model can also obtain negative values, when the intercalate density is very high (e.g. MC₂). These structures have not been considered in the trends of intercalate diffusion rate. In this analysis we only consider the structures where the intercalate atoms do not overlap based on their atomic radius, which were found to be empirical formulas MC₆, MC’₈ and MC₈, as the nearest neighbour interatomic intercalate distances were 4.266 Å, 4.524 Å and 4.524 Å, respectively. These interatomic distances are larger than two times the atomic radius of the intercalates which are 3.24 Å, 3.6 Å, 4.5 Å, 4.8 Å and 5.4 Å for Li, Na, K, Rb and Cs, respectively. Here, we find that for the other structures MC₂ and M₃C₈ and MC₃ the nearest neighbhour intercalate distance is only 2.262 Å which is less than two times the atomic radius of any of the group 1 intercalates. As such structure would result in a 2D metallic state of the intercalate atoms, and would render such a material not applicable in graphite-based battery materials. Based on this analysis we find that the trends of intercalate diffusion rate for the various group 1 elements having the MC₆ empirical formula are,

LiC₆=0.0154 > NaC₆=0.0124 > KC₆=0.0056 > RbC₆=0.0044 > CsC₆=0.0011

This trend suggests that lithium as an intercalate is the GIC that would have the largest diffusion rate of intercalate atoms, but interestingly a NaC₆ GIC would only be 20% less efficient, based on geometric considerations. So we provide here a computationally very efficient way to screen various GIC’s for their efficiency as ion battery materials using as an input only the atomic radii of atoms in the GIC and using the sphere-in-contact model. For the empirical formulas of GIC considered we provide the value of the packing fraction in Table 3 that can be used to screen for the diffusion rate of any GIC with AαAα stacking given only the radius of the intercalates and carbon atoms.

) of the various empirical formulas of GIC using the sphere-in-contact model.

Formula	Analytical equation of density	ηLi, GIC	ηNa, GIC	ηK, GIC	ηRb, GIC	ηCs, GIC
MC2	${\displaystyle \frac{\frac{4}{3}·\pi ·{r_M}^3+\frac{8}{3}·\pi ·{r_C}^3}{{12\sqrt{3}·r_C}^2\sqrt{{\left(r_C+r_M\right)}^2-4·{r_C}^2}}}$	1.073	1.262	1.936	2.095	2.623
MC6	${\displaystyle \frac{\frac{4}{3}·\pi ·{r_M}^3+\frac{24}{3}·\pi ·{r_C}^3}{36\sqrt{3}{{·r}_C}^2\sqrt{{\left(r_C+r_M\right)}^2-4·{r_C}^2}}}$	0.461	0.513	0.717	0.767	0.936
M3C8	${\displaystyle \frac{4·\pi ·{r_M}^3+\frac{32}{3}·\pi ·{r_C}^3}{{48\sqrt{3}·r_C}^2\sqrt{{\left(r_C+r_M\right)}^2-4·{r_C}^2}}}$	0.843	0.981	1.479	2.130	1.991
MC3	${\displaystyle \frac{\frac{8}{3}·\pi ·{r_M}^3+\frac{24}{3}·\pi ·{r_C}^3}{36\sqrt{3}{·r_C}^2\sqrt{{\left(r_C+r_M\right)}^2-4·{r_C}^2}}}$	0.383	0.444	0.663	0.716	0.890
MC’8	${\displaystyle \frac{\frac{4}{3}·\pi ·{r_M}^3+\frac{32}{3}·\pi ·{r_C}^3}{{48\sqrt{3}·r_C}^2\sqrt{{\left(r_C+r_M\right)}^2-4·{r_C}^2}}}$	0.385	0.419	0.565	0.602	0.725
MC8	${\displaystyle \frac{\frac{4}{3}·\pi ·{r_M}^3+\frac{32}{3}·\pi ·{r_C}^3}{{48\sqrt{3}·r_C}^2\sqrt{{\left(r_C+r_M\right)}^2-4·{r_C}^2}}}$	0.385	0.419	0.565	0.602	0.725

1.9. Calculation of the Interatomic Repulsive Energy (U_net) Between Intercalates to Assess the Stability of an Arbitrary GIC

To calculate the interaction between intercalates we have used a supercell approach for the various GICs in Figure 1, shown in Figure 5. We also assumed that the interactions between intercalates are up to the 8^th nearest neighbour, which is about 8 Å, depending on the atomic radius of the intercalate. This assumption is valid as there will be significant screening of the partial charges on the surface of the atoms by the surrounding electron clouds, of nearby intercalate atoms. We have counted all 8 nearest neighbour interactions (i.e. d_i, where i = 1, 2, 3, 4, 5, 6, 7, 8) of the various configurations shown in Figure 5, which are shown in Table 1. We then used the following trigonometric and geometric derivation to calculate the repulsive interaction between partial point charges, at the nearest distance between the intercalate spheres. A similar equation was used to calculate the repulsive interaction between adsorbates on HCP and FCC, metal surfaces [Ref, preprint]. This equation is a summation of all the pairwise interactions, up to the 8^th nearest neighbour intercalates and we assume that the partial charge on the surface of the intercalate is constant, and does not depend on the interactions between the intercalate (which would generate a local electric field), but only depends on the coordination number (CN) that the intercalate has with the graphene layer (i.e. the contact points of the intercalate sphere with the carbon atoms). Since the intercalates in all GIC’s studied in this study reside at the middle of the hexagon, in each carbon ring, and the stacking of the layers is AαAα, we can calculate the internal energy of this GIC model, depending on the saturation with lithium ions, during charging and discharging of the battery material. The coordination number of the intercalate is 12, as there are 6 carbon atoms in the first layer and 6 carbon atoms in the second layer, that are in direct contact with the intercalate. As lithium is a group 1 element and it tends to loose one electron to become more stable, the charge on each lithium will be +1. But since some of this electronic charge of the ionised electron, is charge transferred to the carbon rings in direct contact with the intercalate, and then from a π-ion bond, therefore the partial charge on the intercalate will only be a fraction of +1. Perhaps, we can assume that it will be +1/12 e because the coordination number of the intercalate is 12. Here we note that the absolute value of the partial charge on the intercalate surface is not a critical parameter in the sphere-in-contact model/theorem, in order to get the trends of the interactions between intercalates, in the various GICs. So we used a partial charge of +0.08 e and +0.16 e (this is a positive partial charge) gives us an indication of the repulsive interaction energy trends between intercalates in the various GICs. This provides a trigonometric equation that quantifies the stability of the various GICs and can be used to assess the internal energy (U_net) associated to various GICs, in a charged state, by changing the partial charge on the intercalates. This equation has dependent variables the atomic radius of the carbon atom ($r_C$) and the intercalate ($r_M$) and the partial charge on the intercalate surface ($\delta$). These for Li_xC_y GICs are $r_C=0.709\times 10^{-8}cm$ and $r_{Li}=1.52\times 10^{-8}cm$ , respectively.

The distances between the intercalates are given by the following trigonometric equations. These are as a function of the atomic radius of the carbon atom and the atomic radius of the intercalate. So they can be applied to any GIC and are,

$$d_1=2(\sqrt{3}{r}_C-r_M)$$

(12)

$$d_2=2(3r_C-r_M)$$

(13)

$$d_3=2(2\sqrt{3}{r}_C-r_M)$$

(14)

$$d_4=\sqrt{84}{r}_C-{2r}_M)$$

(15)

$$d_5=2(3\sqrt{3}{r}_C-r_M)$$

(16)

$$d_6=2(6r_C-r_M)$$

(17)

$$d_7=\sqrt{156}{r}_C-{2r}_M)$$

(18)

$$d_8=2(4\sqrt{3}{r}_C-r_M)$$

(19)

The interaction energy can be calculated by the pairwise summation using Coulombic terms given by,

$$U_{net}\left(kJ{mol}^{-1}\right)={\displaystyle \frac{N_{Av}·k{·\delta }^2}{1000}}\sum _{i,j=1}^8{\displaystyle \frac{n_i}{d_i}}$$

(20)

where $k$ = 8.98755 x 10⁹ Nm²C^-2, e = 1.60217 x 10^-19 C, r_Li = 1.52 x 10^-10 m, r_C = 0.709 x 10^-10 m and $N_{Av}=6.022·{10}^{23}{mol}^{-1}$. So for partial charges on the intercalates of 0.08 e, δ = 0.128 x 10^-19 C this formula becomes

$$U_{net}\left(kJ{mol}^{-1}\right)=88.675·{10}^{-11}\left(kJ{mol}^{-1}\right)\sum _{i,j=1}^8{\displaystyle \frac{n_i}{d_i}}$$

(21)

and for partial charges of the intercalates of 0.16 e, δ = 0.256 x 10^-19 C

$$U_{net}\left(kJ{mol}^{-1}\right)=3.547·{10}^{-9}\left(kJ{mol}^{-1}\right)\sum _{i,j=1}^8{\displaystyle \frac{n_i}{d_i}}$$

(22)

To establish an initial estimate of the pairwise sphere–sphere interactions, 2 researchers performed independent manual counts (i.e. using the pen and paper approach) using the structural image from Figure 5. Distances were approximated by comparing the spacing between spheres relative to known lattice features. Both sets of manual counts were recorded in a separate table, allowing us to compare consistency across independent observers and identifying the counts that had to be repeated, to achieve agreement between the 2 observers.

The centre of each sphere was identified and the distance of adjacent spheres was estimated based on the length of the C-C bond, the vertical and horizontal width of the benzene ring and whether the distance vector was on one of the crystallographic lattice vector or offset by one hexagonal ring, of the graphite layer. So with the use of these lattice features we were able to get very accurate counts of the pairwise interactions in combination with the use of a large supercell, to ensure that the repulsive energy of interaction between intercalates is as accurate as needed to be able to estimate the relative trend of stability in the these GIC’s, and its nearest neighbours were determined based on perceived based on nearest neighbour perceived surface-to-surface proximity.

With this approach we have listed these counts in Table 4 along with the corresponding repulsive interaction energy for two partial charge values. The partial charges used on the intercalate where 0.08 e (i.e. U_{net, δ = +0.08e} ) and 0.16 e (U_{net, δ = +0.16e}). We also attempted counting pairwise interactions using an image analysis software that extracted the centroid coordinates of the sphere in the images in Figure 5 and the sorted the distance between centroid cartesian coordinate with a computer script. But these counts were not found to be correct in most cases and therefore are not listed here.

Based on Table 4 we can observe that the relative stability trends of the various GIC do not depend on the amount of partial charge on the intercalates. We also observe that LiC₆, LiC’₈ and LiC₈ are considerably more stable by almost one order to magnitude in kJ mol^-1 than the higher intercalate density GIC’s.

So the trends in relative internal energy of the various GIC found is

LiC’₈ < LiC₈ < LiC₆ < LiC₃ < Li₃C₈ < LiC₂

Here we not the stability trend calculated using the sphere-in-contact model is very reliable as we found that the number of pairwise interactions considered by only considering a two-dimensions grid of the partial charges, is as accurate as considering that the partial charge on the intercalates is smeared over the surface of a sphere. In general we found that the repulsive energy term converges once the number of partial charges in a 2 dimensional circle is larger than 8. We also found that considering a circle to distribute these partial charges, rather than a sphere causes only overestimation of the repulsive energy between two interaction circles with equally spaced partial charges around the periphery by 9.2%.