Cost-Effective Calculation of Collective Electronic Excitations in Graphite Intercalated Compounds.

Graphite/graphene intercalation compounds with good and improving electrical transport properties, optical properties, magnetic properties and even superconductivity are widely used in battery, capacitors and so on. Computational simulation helps with predicting important properties and exploring unknown functions, while it is restricted by limited computing resources and insufficient precision. Here, we present a cost-effective study on graphite/graphene intercalation compounds properties with sufficient precision. The calculation of electronic collective excitations in AA-stacking graphite based on the tight-binding model within the random phase approximation framework agrees quite well with previous experimental and calculation work, such as effects of doping level, interlayer distance, and interlayer hopping on 2D π plasmon and 3D intraband plasmon modes. This cost-effective simulation method can be extended to other intercalation compounds with unlimited intercalation species.

1. Introduction

The physical properties can be modified by intercalation significantly in graphite or graphene intercalated compounds (GICs) [1,2], such as electrical transport properties [3], optical properties [4], magnetic properties [5] and even superconductivity [6,7,8,9,10]. Superconductivity in the alkali GICs was first reported in 1965 [6], with critical temperature below 1 K. After the critical temperature rising to 11.5 K of CaC was reported [7], great interest has arisen for electronic properties of GICs [11,12,13,14,15,16].

The electronic collective excitations play a significant role in optical and dielectric properties, which have been studied both in experiments [17,18,19,20] and in calculations [21,22,23,24,25,26,27] in GICs. The low-energy plasmon mode, which exists in pure graphite at 7 eV in optical limit, was found to shift to lower energy in GICs [17,18,19,20,22,23,24]. In addition, an intraband plasmon (IP) mode appeared with energy ∼1 eV in GICs due to the doping effects [18,19,20,21,23,24]. Furthermore, acoustic plasmon (AP) may exist and could play a role in superconductivity [20,23,24].

Theoretically, the plasmon properties are studied in a tight-binding (TB) model and time-dependent density functional theory (TDDFT). Shung [21] reported that the in-plane IP mode can be well described in a layered 2D two-band TB model near the Dirac cone, where the interlayer tunneling effects are neglected and only Coulomb interaction in different layers is retained. Lin et al. [22] revealed the in-plane plasmon properties can also be depicted in this model, where the bandstructure in the whole first Brillouin zone (1BZ) has to be taken into account. Echeverry et al. [23,24] studied the dielectric properties of GICs at energies below 12 eV by the state-of-the-art first-principle calculations. Bulk plasmon (BP), plasmon, IP and AP modes are discussed in LiC, CaC, SrC, and BaC, which is in good agreement with the experimental results. However, the first-principle calculation is too time-consuming, and when the species of intercalation changes, the band structure, interlayer distance, and doping level are changed simultaneously, and it is hard to distingulish their effects on plasmons.

Computational simulation is a green tool for new materials developing and property prediction which could be used in many research studies. However, limited computing resources and insufficient precision restricted the promotion of wide usage. In the present work, in order to understand how the in-plane and out-plane electronic collective excitations are affected by different parameters in GICs, we extend the two-band TB model of AA-stacking graphite, where the interlayer tunneling effects are also taken into consideration. In this model, we study the low-energy electronic collective excitations within the random phase approximation (RPA) framework. The local field effects (LFE) are also involved. Our calculations show that doping level and interlayer hopping affect plasmons mainly by band structure effects, such as density of states (DOS) and group velocity near the Fermi level, and the forbidden effect of interband transition, while interlayer distance modifies plasmons via the long-range interlayer Coulomb correlations. In addition, stacking order can modify interlayer coupling and have a great influence on plasmon properties.

2. Method

For periodic systems, the dielectric function can be written in terms of tight-binding basis in the form where G,G’ are reciprocal-lattice vectors, q is a wave vector restricted to the first Brillouin zone (1BZ), and is the Fourier transform Coulomb potential, being the volume of the unit cell. The polarizability function is expressed as [28] where and are the eigen-energy and Fermi–Dirac occupuation for band index n and wave vector k, is the contribution of the th tight-binding basis function to the Hamiltonian eigenstate (for more details, please see Appendix A), N is the number of unit cells, is a broadening parameter, and is a phase factor. Factor 2 accounts for spin (we assume a spin-degenerate system). is defined as the charge-density wave. The index s stands for the lattice vector index L and for the indices and of the orbital. Using these relations, the dielectric matrix can be written in the form

The separable form of the susceptibility matrix in Equation (5) enables us to calculate the inverse dielectric matrix [29] where and where the Coulomb interaction between the charge-density waves is [29]

The off-diagonal elements of the matrix describes the response of the electrons at wave vectors different from the external perturbing field and thus contain information about the inhomogeneity of the microscopic response of electrons known as the local field effect (LFE) [30]. The macroscopic dielectric function is defined as

If we neglect the LFE, it becomes . This macroscopic dielectric function is directly related to many experimental properties. For example, the optical-absorption spectrum (ABS) is given by Im. The electron energy-loss spetrum (EELS) is proportional to ). EELS is especially useful in probing the collective electronic collective excitations, known as plasmons, of bulk and low-dimensional systems.

Without loss of accuracy, we make some approximations to avoid calculating the integral of tight-binding basis function, as shown in Appendix B, to accelerate calculation greatly.

3. Results and Discussion

3.1. TB Model of AA-Stacking Graphite

In AA-stacking graphite, layers of carbon atoms locate directly on top of each other, as shown in Figure 1a. In order to obtain a reliable TB model, the geometrical optimization and electronic properties are performed by first-principle calculations (for more details, please see Appendix C). The calculated band structure in the vicinity of the Fermi level along the high-symmetry directions of 1BZ is shown in Figure 1b, which is in good agreement with the previous results [31]. The bonding and antibonding bands dominate the band structure. A hole pocket appears at K and an electron pocket appears at H.

Figure 1

Lattice and electronic structure of AA-graphite. (a) Top view of atomic structure of AA-stacking graphite. The primitive cell, lattice vectors a, a and in-plane hopping t are indicated. (b) Band structure of AA-stacking graphite, with the corresponding orbitals in cyan. The size of the symbol represents the orbital weight. (c) Comparison of the band structure for AA-stacking graphite obtained from DFT (red line) and TB (blue dashed line) calculations, respectively.

After fitting the TB parameters, we constructed TB Hamiltonian as where is the on-site energy of orbital of C atoms (0.51 eV in this work), and

The detail values of hopping parameters, schematically shown in Figure 1a, are listed in Table 1. From this TB Hamiltonian, energy eigenvalues at high-symmetry points can be calculated analytically (listed in Table 2). The band structures derived from our TB model match well with DFT calculation (Figure 1c), implying a reliable TB model.

Table 1

Hopping parameters (in eV) assigned to the simple TB Hamiltonian of AA-stacking graphite. d is the distance between the atomic sites on which the interacting orbitals are centered.

i	t_i (eV)	d (Å)
1	−3.24	1.41
2	0.36	2.44
3	−0.41	2.82
4	0.095	4.23
⊥	0.21	3.7

Table 2

The analytical energy eigenvalues at high-symmetry points.

Points	Energy Eigenvalues
Γ	ϵ_p+2t_⊥+6t_2+6t_4±3(t_1+t_3)
K	ϵ_p+2t_⊥
M	ϵ_p+2t_⊥−4t_2−2t_4±32(t_1−3t_3)
A	ϵ_p−2t_⊥+6t_2+6t_4±3(t_1+t_3)
H	ϵ_p−2t_⊥
L	ϵ_p−2t_⊥−4t_2−2t_4±32(t_1−3t_3)

3.2. Plasmon Excitations

By the method introduced in Section 2, we calculated the excitation spectra of AA-stacking graphite with and without inclusion of the LFE. K, eV, [32] (more details in Appendix B) and a dense k-mesh were used in all calculations. In Figure 2, the loss function is shown as a function of energy at variable momentum transfer q along M (a) and A (b) directions.

Figure 2

Calculated loss functions along (a) the M and (b) A directions as a function of q. The peak positions are marked by red stars. Real and imaginary part of dielectric function at (c) Å and (d) Å along the M direction.

In the M (in-plane) direction, the main feature of the loss spectra is a strong high-energy peak around 7 eV at small q. As q increases, it shows a parabolic-like positive dispersion and splits into two peaks at large momentum ( Å). This peak is attributed to the collective interband transitions from to band around M point and it is assigned as plasmon [33]. The split features of plasmon have also been reported in monolayer graphene [34]. With the help of real and imaginary parts of the dielectric function at Å (Figure 2c) and Å (Figure 2d), it is clear that the peaks splitting at large q originates from the splitted collective interband transitions. In addition, there exists a low-energy intraband plasmon (IP) mode with quite low intensity near Dirac points corresponding to low doping level of AA-stacking graphite.

In the A (out-plane) direction, the main peak starts at 0.26 eV, which originates from the collective intraband transitions. As q increases, the position of this peak increases and reaches its maximum 0.84 eV at the boundary of 1BZ with Å, and then, it decreases. The intensity of the peak also increases firstly and then decreases as q increases, but it reaches its maximum at Å.

3.3. Effect of Doping Level

Intercalating different atoms in graphite will induce different doping levels, which is important to understand the behavior of electronic collective excitations. The induced electrons will occupy the band of carbon atoms, so the doping level can be regarded as an in-plane quantity to some extent. In order to study how this in-plane parameter influences the in-plane and out-plane collective excitation of electrons, we calculate the EELS of AA-stacked graphite at different doping levels. In this calculation, the doping effect is represented by ragid band approximation (RBA), which means only a rigid shift of the Fermi level, as shown in Figure 3. We calculated loss functions at different doping levels of 0.08, 0.25, 0.333, 0.42 and 0.667 e-doped per unit cell, correspoing to , , , , and cm charge carrier density addition, respectively. The quite small difference between Figure 4a,b, especially at small q, illustrates that LFE is negligible in this direction. As the doping level increases at small q, both energies and intensities of IP mode increase, while plasmon performs the inverse. When the doping level is larger than 0.42 e-doped per unit cell, the Fermi level will be lifted to the Van Hove singularity shown in Figure 3, and the plasmon disappears. Because in this case, the band at M point is occupied and the corresponding interband transition is forbidden.

Figure 3

Band structure and Density of States (DOS) at different doping levels in rigid band approximation (RBA). The Fermi level at 0.08, 0.25, 0.333, 0.42, and 0.667 e-dope per unit cell are marked by light blue, orange, green, brown and violet dashed lines, corresponding to at 0.98, 1.555, 1.777, 1.992, and 2.793 eV, respectively.

Figure 4

Calculated loss functions along the M direction (a) with LFE, (b) without LFE and along the A direction (e) with LFE, (f) without LFE at different doping levels (light blue line for 0.08 e-doped, orange line for 0.25 e-doped, green line for 0.333 e-doped, brown line for 0.42 e-doped and violet line for 0.667 e-doped per unit cell). (c,d,g,h) show positions of peaks of (a,b,e,f), respectively. The shaded background gives the loss functions of the no-doping case for comparison.

The dispersions of the in-plane modes with different electron doping levels are summarized in Figure 4c,d, corresponding to with and without inclusion of LFE, respectively. The plasmon shows a parabolic-like positive dispersion with the increase of momentum transfer. As the doping level increases, at small q, more interband transitions are forbidden and the plasmon energy shows a redshift; however, at large q, the forbidden effect is negligible and the plasmon energy shows a blueshift. As a result, a higher doping level leads to a stronger parabolic-like positive dispersion of P. The IP mode also shows a parabolic-like positive dispersion with the momentum transfer increase at low doping level. A higher doping level lifts the IP energy but reduces the dispersion, and even negetive dispersion appears at small q in 0.42 e-doped case. The negative dispersion of the in-plane IP mode has also been reported in CaC[23] and SrC [24], which is explained by band structure effects.

In layered materials, plasmon properties are anisotropic and quite different between in-plane and out-plane direction. From Figure 4, it is obvious that LFE is crucial in A direction, especially for the dispersion as shown in Figure 4g,h. When LFE is neglected, the energies of the out-plane IP mode decline to nearly 0 eV at the boundary of the second Brillouin zone (2BZ), while when LFE is taken into account, this IP mode shows a nearly flat band. The variation in energy is 0.14, 0.08, 0.09, 0.05 and 0.06 eV in 0.08, 0.25, 0.333, 0.42, and 0.667 e-doped per unit cell, respectively. This result agrees quite well with previous work in calculation for CLi [24]. As the doping level lifts, the plasmon energy increases at first, reaches maximum around 0.333 e-doped per unit cell, and then drops down. The variation of intensities as a function of doping level shares the trend with the plasmon energies. Both density of states near the Fermi surface and the group velocities determine the IP energy. As shown in Figure 3, at low doping level, the density of state increases drastically with doping level increasing, and it reaches the top at 0.42 e-doped per unit cell. So, the IP energy firstly increases with doping level increase; then, it decreases due to the drop of density of state above Van Hove singlarity. However, the average group velocity along the A direction of 0.42 e-dope is less than that of 0.333 e-dope. As a result, the IP energy of 0.42 e-dope is less than that of 0.333 e-dope.

3.4. Effect of Interlayer Distance and Hopping

Intercalating different atoms in layered materials will induce different interlayer distance and hopping, which are two important out-plane parameters to understand the behavior of electrons. In reality, they are always changed simultaneously, and it is hard to study their effect independently in experiments and first-principle calculations. In order to establish a simple picture how these two out-plane parameters influence the in-plane and out-plane electronic collective excitations, we tune them independently in our TB model and calculate the corresponding EELS.

If we fix the interlayer hopping and tune the interlayer distance d, the band structrue has no change except for rescaling of reciprocal lattice vector along -axis. So, the in-plane plasmon will be only tuned by the interlayer Coulomb correlation, while the out-plane plasmon can also be modified by rescaling effect along the -axis. For example, as shown in Figure 5, the Fermi level is fixed at 1.555 eV, and the interlayer hopping is fixed at 0.21 eV. When d increases, the interlayer Coulomb interaction becomes weaker, so the in-plane IP and plasmon energy decrease (Figure 5a–c). The intensity of IP increases as d increases, while that of plasmon decreases. Figure 5b,c show the peak positions of the loss spectra along the M direction with and without LFE, respectively. With the increase of d, both energies of IP and plasmon decrease, but with different rates obviously. When d increases from 2.44 to 4.88 Å, the plasmon energy decreases from 8.1 to 6.2 eV at small q, while the IP energy decreases only from 2.6 to 2.4 eV. As q increases, both IP and plasmon show parabolic-like positive dispersion at different interlayer distances. With the interlayer distance decreasing, dispersions of IP mode remain unchanged with and without the inclusion of LFE, and we find a weaker dispersive feature of plasmon mode with the inclusion of LFE. In addition, we plot the peak positions of the loss spectra along the A direction in Figure 5d,e, where the out-plane IP energies indeed increase as d increases due to the -axis rescaling effect.

Figure 5

(a) Calculated loss functions along the M direction with LFE at different interlayer distance (light blue line for Å, orange line for Å, green line for Å, brown line for Å, violet line for Å and brown line for Å). The shaded background gives the loss functions of Å case for comparison. (b,c) show the loss function peaks along the M direction with and without LFE, respectively. (d,e) show the loss function peaks along the A direction with and without LFE, respectively. The Fermi level is fixed at eV.

As shown in Figure 6, larger interlayer hopping results in lower Fermi level and a larger energy difference between the different points along the -axis, creating a larger group velocity along the out-plane direction. As a result, both the energy and weight of the out-plane IP mode increase dramatically with the increase of interlayer hopping , as shown in Figure 7e–h. The effect of changing interlayer hopping on the band structure is similar to tuning the doping level at different plane seperately. For example, the energy difference between the MK plane and AHL plane is , as shown in Table 2. Therefore, its effect on plasmons can be attributed to a combination of different doping level in all planes in 1BZ. Figure 7a–d shows the loss function along the M direction (in-plane). The increase of has little effect on the in-plane plasmon mode. The energy of the in-plane IP mode increases as increases, while the weight decreases significantly.

Figure 6

Band structure and density of states (DOS) at different interlayer hopping (red line for = 0.21 eV with = 1.555 eV and blue dashed line for = 0.50 eV with = 1.428 eV). The interlayer distance is fixed at 3.70 Å and the carrier density is fixed at cm.

Figure 7

Calculated loss functions along the M direction (a) with LFE, (b) without LFE and along the A direction (e) with LFE, (f) without LFE at different interlayer hopping (light blue line for eV, orange line for eV, green line for eV and brown line for eV). (c,d,g,h) show positions of peaks of (a,b,e,f), respectively. The shaded background gives the loss functions of the no doping case for comparison.

3.5. Effect of Stacking Order

In all known GICs, the stacking order of graphene layer is AAA [1], which is different from the AB or ABC stacking order in natural graphite. The stacking order can modify the interlayer coupling and electronic structure consequently. For example, in small-rotation-angle twisted bilayer graphene, the moiré superlattice alters the electronic properties significantly and has led to observations of exotic emergent electronic properties such as superconductivity and strong correlated states [35,36].

The plasmon properties of AB and ABC stacking graphite are studied based on our TB model (for more details, please see Appendix D). The loss function of AB graphite along the M direction is almost the same as that of AA graphite. The main feature of the loss spectra is a strong plasmon peak around 7 eV at small q, and the plasmon shows positive-parabolic dispersion with an increase of q, which agrees very well with previous work [33], as shown in Figure 8a. The real and imaginary parts of dielectric functions at Å in Figure 8b and at Å in Figure 8c are similar with that of AA graphite, except that no zero point of Re[] appears at eV due to the lack of charge carriers. However, no plasmon peaks are observed in ABC graphite in the low-energy region (0–12 eV), as shown in Figure 8d. The collective character in the loss function is cofirmed by the real and imaginary parts of the dielectric function at Å in Figure 8e and Å in Figure 8f. We can see that at energy around 7 eV, Re[] crosses the zero line with a positive derivative at Å however, Im[] is still too large, implying a strong Landau damping effect. As a result, plasmon will have a very large linewidth, and no peak appears in the loss function. Along the out-plane A direction, no prominent plasmon peaks appear in the low-energy region (0–12 eV) for both AB and ABC stacking graphite due to the low level of charge carrier concentration.

Figure 8

Calculated loss functions along the M direction of AB (a) and ABC (d) graphite with LFE. The peak positions are marked by red stars. Real and imaginary part of dielectric function at (b) Å and (c) Å for AB graphite and (e) Å and (f) Å for ABC graphite.

As the doping level increases, the density of state and the group velocity near the Fermi surface change rapidly according to the band structure, and more interband transitions from band to band are forbidden at small q. As a result, with the increase of doping level, the energy blueshifts for in-plane IP mode, blueshifts at first and then redshifts for out-plane IP mode, and blueshifts at small q but redshifts at large q for in-plane plasmon mode. The variation in the intensity of plasmon peaks is the same as plasmon energies. When the doping level is so large that the band at the M point is occupied, the in-plane plasmon disappears as a result of a forbidden interband transition from the band to band.

With the interlayer distance increasing, the band structrue has no change except for rescaling of the reciprocal lattice vector along the -axis when the interlayer hopping remains unchanged, and the interlayer Coulomb correlation becomes weaker. As a consequence, the energy redshifts slightly for the in-plane IP mode, redshifts remarkably for the in-plane plasmon mode, and blueshifts for the out-plane IP mode. The variation in the strength of plasmon peaks is the same as plasmon energies for the in-plane plasmon and out-plane IP modes, but it is opposite for the in-plane IP mode. The effect of changing interlayer hopping on the band structure is similar to tuning the doping level at different planes separately and group velocity along the out-plane direction. We find that the interlayer hopping has nearly no effect on in-plane plasmon mode, but it affects the IP mode dramatically. As the interlayer hopping increases, the energies of both the in-plane and out-plane IP mode blueshift, and the strength of plasmon peak increases for the out-plane mode but decreases for the in-plane mode. AB-stacking graphite shows few differences with AA-graphite, while no plasmon peaks appear in ABC-stacking graphite because of the strong Landau damping effect.

A comparison of the excitation spectra obtained with and without inclusion of the local-field effects demonstrates that in layered AA stacked graphite, the LFE have a significant impact on the dielectric properties, especially along the out-plane direction, where they flatten the out-plane IP dispersion at least in the second Brillouin zone (2BZ).

4. Conclusions

Based on the computing resource and precision limitation of current simulation on layered materials, we have constructed a TB model of AA-stacking graphite to mimic GICs and investigated the plasmons properties within the RPA framework, with high efficiency and sufficient precision. The local-field effects are involved in our model; the effects of doping level, interlayer distance, interlayer hopping, and stacking order on 2D plasmon and 3D intraband plasmon modes are studied independently; and the corresponding evolutions of plasmons are presented clearly. Our results are in very good agreement with the previous experimental and calculation work. This tight-binding calculation does not need a self-consistent process and the basis set contains only several orbitals, so it is very computationally efficient. At the same time, the precision is higher than many methods and satisfies most of the prediction requirements. Our method is easy to extend to other intercalation materials with unlimited intercalation species and has great significance for understanding how plasmons are tuned in layered materials.

Table A1

TB parameters for AB stacking and ABC stacking graphite. All values are in eV. The TB Hamiltonian is valid in the whole three-dimensional Brillouin zone of graphite.

Stacking Order	AB	ABC
t_1	−3.04	−3.042
t_2	0.254	0.267
t_3	−0.312	−0.313
t_4	0.06	0.066
t_5	−0.022	−0.018
γ_0	0.177	0.177
γ_1	0.098	0.048
γ_2	0.045	−0.042
γ_3	−0.035	0.084
γ_4		−0.024
E_0	0.317	0.333
Δ	0.004	0