Graphite Under Compression: Shift of Layer Breathing and Shear Modes Frequencies with Interlayer Spacing

Layered materials have huge potential in various applications due to their extraordinary properties. To determine the interlayer interaction (or equivalently the layer spacing under different perturbations) is of critical importance. In this letter, we focus on one of the most prominent layered materials, graphite, and theoretically quantify the relationship between its interlayer spacing and the vibrational frequencies of its layer breathing and shear modes, which are measures of the interlayer interaction. The method used here to determine the interlayer interaction can be further applied to other layered materials.

Layered materials have huge potential in various applications due to their extraordinary properties. To determine the interlayer interaction (or equivalently the layer spacing under different perturbations) is of critical importance. In this letter, we focus on one of the most prominent layered materials, graphite, and theoretically quantify the relationship between its interlayer spacing and the vibrational frequencies of its layer breathing and shear modes, which are measures of the interlayer interaction. The method used here to determine the interlayer interaction can be further applied to other layered materials. Layered materials have unique electronic and excellent mechanical properties [1]. Many of the properties are closely related to the weak interlayer van der Waals (vdW) interaction between the 2-D molecular monolayers [1]. Vibrational modes include shear modes (CMs) and layer breathing modes (LBMs), due to relative motions of the planes. The CMs (vibrations parallel to the planes) have been experimently identified in many bulk layered materials, such as h-BN [2], MoS 2 [3] and WSe 2 [4]. The LBMs (vibrations perpendicular to the planes) are less studied, because they are optically inactive. Nevertheless, both modes are of significant importance to the understanding of various layered materials (and the full exploitation of their application potential), as direct measures of the interlayer interactions.
The interlayer interactions (and the frequencies of the CMs and LBMs) of layered materials strongly depend on the interlayer spacing. So do many of their properties. Therefore, the control of the interlayer spacing (strain) is vital to the applications -one could adjust the properties as needed by tuning the interlayer spacing. To quantify the relationship between the frequencies of the interlayer modes and the interlayer strain provides a practically convenient way to determine the strain as the phonons are usually easier to directly detect than the lattice constant in experiments. Also this relationship itself is a fundamental mechanical property for each material.
Among layered materials, we are most familiar with graphite, with just four carbon atoms in a unit cell. Graphite has been studied for over two thousand years. The relatively recent successful synthesis of graphene [5], along with the massive research that followed brought the understanding of graphite to a new level. Yet much remains behind a veil.
The CM of graphite is a Raman active E 2g mode and was measured at 42 cm −1 in 1975 [6]. Tan et al. measured the CMs of two-to eight-layer graphene and bulk graphite. These frequencies fitted well with a linear chain model [7] : where ω N is in cm −1 , N is the number of layers, c is the speed of light in cm s −1 , α CM ∼ 12.8 × 10 18 Nm −3 is the only fitting parameter, referring to the interlayer force constant for the CMs, and µ = 7.6 × 10 −27 kgÅ −2 is the mass per unit area of the single-layer graphene. For graphite, N → ∞ and ω ∞ = 1 πc α CM µ . The shift of the CM with pressure was measured by Hanfland et al. [8]. They measured the frequency at ambient pressure at 44 cm −1 and fitted their data under pressure with the following equation [8]: where δ 0 is the logarithmic pressure derivative (d lnω/dP ) P =0 and δ is the pressure derivative of d lnω/dP . They obtained δ 0 =0.110 GPa −1 and δ =0.43. The pressure-induced change of the CM frequency is intrinsically due to the change of the lattice constants, mainly a 33 . The Grüneisen parameter is commonly used to relate phonon frequencies to volume [9]. For layered materials of large anisotropy, such as graphite, the interlayer modes and intralayer modes are independently related to out-of-plane strain and in-plane strain, respectively, by a scaling parameter γ [10]: There is one scaling parameter for each mode. The scaling parameter and the Grüneisen parameter are equivalent for 3-D isotropic materials. For the CM, Hanfland et al. gave γ=1.4, for the pressure range up to 14 GPa [8].
The LBM of graphite is an optically inactive B 2g mode and was measured at 127 cm −1 by inelastic neutron scattering [11]. The LBM of multilayer Bernalstacked [12] graphene was measured by combination Raman (LO+ZO ) [13] and has not been directly probed so far at room temperature due to the low electron-phonon coupling and its symmetry. Wu et al. directly measured the LBM of twisted multilayer graphene [14] and similarly to the CM, fitted the position with the linear chain model [15]: where N is the number of layers. The shift of the LBM with pressure was measured by Alzyab et al. (the same neutron scattering experiment as mentioned above) for the pressure range up to 2 GPa and a shift rate of ω =19 cm −1 GPa −1 was observed [11]. From the shift rate, δ 0 = (ω /ω 0 )=0.15 GPa −1 were obtained. No further study on the Grüneisen parameter of the LBM of graphite or multi-layer graphene has been reported.
In this paper, we modeled graphite under hydrostatic compression and uniaxial compression along c-axis. We calculated the phonon frequencies of the CMs and LBMs in each case and obtained the Grüneisen parameters for each mode.
We used density functional theory (DFT) [16,17] as implemented in the Vienna ab initio Simulation Package (VASP) [18] with the generalized gradient approximation as parameterized by Perdew, Burke and Ernzerhof [19]. We used 900 eV plane-wave cutoff energy and the reciprocal unit cell meshed with 18×18×9 k-points. We included the vdW using the Grimme method [20] as implemented in the VASP code. We used the 2×2×2 supercell employing the finite displacement method as implemented in the PHONOPY code [21] to calculate the phonon frequencies. The calculations setup were the same as that used in our previous work [22], where more details can be found.
We started by obtaining the optimized geometry of unstrained graphite, with the in-plane bond length a 0 =1.42 A and the interlayer distance c 0 =3.20Å. The errors relative to the experimental values are 0.06% and 4.6%, respectively [23]. We expect a relatively large error in the interlayer spacing, despite the vdW add-on being used. To minimize the effect of the inaccuracy of the vdW, we focused on graphite under compressive strain, where the repulsion dominates over the vdW attraction. Furthermore, the shift of the phonon frequencies with compressive strain is mainly due to the increasing overlap of π-electrons of neighboring layers, where vdW plays only a small role.
We then modeled graphite under hydrostatic pressure, following the same method as in Ref. [22] -setting a unit cell volume and calculating the corresponding pressure.  Fig. 1, the difference between the CM (L) and CM (T) is merely to be seen and therefore we study the longitudinal mode alone as a representative for the CM. We fit the data under compressive pressure range up to 10 GPa with Eq. 2 and obtained δ 0 =0.1055 GPa −1 , δ =0.3707 for the CM and δ 0 =0.1969 GPa −1 , δ =0.3541 for the LBM. As mentioned before, the values from the experiments are δ 0 =0.110(8) GPa −1 and δ =0.43(3) for the CM (pressure range up to 14 GPa) [8], and δ 0 =0.15 GPa −1 with no error bars for the LBM (pressure range up to 2 GPa) [11]. For the CM, our results agree reasonably well with the experimental values, especially in the initial shift rate. For the LBM, we expect the calculations to be as reliable as the CM, as the mechanism of the shifts of both modes under pressure is the same (increasing interlayer interaction through overlap of π-electrons of neighboring layers). The less satisfying agreement between the calculations of the LBM and experiment is probably due to the lack of reliable experimental data.
The shift of the frequencies of interlayer modes, such as CM and LBM, should be mainly induced by interlayer strain, but we do not know if there is contribution from intra-layer strain. As our previous work showed, the contribution of interlayer strain to intra-layer modes of graphite was non-negligible [22]. To check, we quantify the relationship between the phonon frequencies and non-hydrostatic strain. We modeled graphite under compression along c-axis, by applying uniaxial strain -setting an interlayer distance and fixing the in-plane geometry, and by applying uniaxial stress -releasing the in-plane geometry. We present the data under uniaxial strain first. We plot the phonon frequencies with interlayer distance in Fig. 2 (a). We fit the data with Eq. 3 and obtain γ CM =1.6729 (2) and γ LBM =2.26 (5). The errors are from the fitting. For uniaxial stress, similarly we fit the data with Eq. 3 and obtain γ CM =1.675 (1) and γ LBM =2.27 (5) in Fig. 2 (b). The almost identical values of γ CM and γ LBM from uniaxial strain and stress indicate, that the contributions of the in-plane strain to the CM and LBM are trivial, unlike the non-negligible contributions of the out-of-plane strain to the in-plane GM [22].
We presented the data under hydrostatic pressure at beginning as they were compared to the experimental results to show the reliability of our calculations. Those data can also be used to further validate the trivial contribution of the in-plane strain to the interlayer modes by obtaining the γ under hydrostatic pressure and see if the value is close to the cases of uniaxial strain and stress. We calculated the corresponding interlayer distance to the data in Fig. 1 (the calculation input is the unit cell volume here), and plot the phonon frequencies against it in Fig. 2 (c). We obtain γ CM =1.663 (2) and γ LBM =2.29 (4). We are now confident to conclude that for the frequencies of the CM and LBM of graphite under strain, there is negligible contribution from the in-plane strain. From the fitting in Fig. 2, Eq. 3 describes the shift of the CM with strain excellently with γ CM =1.6729 (2) and the LBM reasonably well with γ LBM =2.26 (5), over the range up to 10 GPa.
For the GM of graphene, many papers used the phenomenological equation proposed by Thomsen et al. to relate the Grüneisen parameters obtained in various conditions [25]: where ω 0 is the unperturbed GM frequency and the SDP is the shear deformation potential. We use γ here to distinguish from the previous scaling parameter γ. Eq. 5 makes explicit the two-dimensional nature of the analysis and later Huang et al. derived this equation from the dynamical equation [26]: where u=(u 1 ,u 2 ) is the relative displacement of the two carbon atoms in the unit cell, and A and B are two . We follow this method to see if it describes the shifts of the CMs and LBMs of graphite better than Eq. 3. Since the in-plane strain makes negligible contribution to the two interlayer modes of graphite, the analysis becomes one-dimensional: where u is the displacement of the layers. We fit the data in Fig. 2 with this Eq. 7 and obtain C CM =3.9 (2)×10 4 , C LBM =5.8 (3)×10 5 under uniaxial strain, C CM =3.8 (2)×10 4 , C LBM =5.7 (3)×10 5 under uniaxial stress, and C CM =3.7 (2)×10 4 , C LBM =5.7 (3)×10 5 under hydrostatic pressure. The unit of C is cm −2 . More sensibly, we obtain the γ in a similar form to that in Eq. 5: where γ =C/(ω 0 (ω + ω 0 )). A caveat must be stated. For the GM, the frequency shift is fractionally small and therefore ω + ω 0 ≈ 2ω 0 , giving a linear relationship between the frequency and strain. To compare the fittings with Eq. 7 to Eq. 3, we subtract straight lines from all the data in Fig. 2 to make the data roughly on horizontal lines, and then vertically shift them for clarity. We plot the fitting curves of all the data with both Eq. 3 and Eq. 7 in Fig. 3. Here the shift rates of the frequencies of the interlayer modes fitted with Eq. 7 drop with increasing strain, opposite to the fitting with Eq. 3. We obtain γ CM = 6.1 and γ LBM = 8.9 under uniaxial strain, over the range up to 10 GPa but clearly this fitting following the phenomenological method is not as good as that with Eq. 3 for the CM. For the LBM, the fitting with Eq. 3 agrees with the data on superlinearity, while Eq. 7 shows sublinearity. Therefore, we think that γ describes the shift of interlayer modes with strain better than γ . Finally, to extend this determination method of interlayer spacing from bulk to multilayer materials, we need to relate the γ to the interlayer coupling strength α in Eq. 1. We do not consider the change in the mass per unit area of single-layer graphene µ as the in-plane contribution is negligible. So it is α = α 0 (a/a 0 ) 6γ , both for the CM and the LBMs.
In conclusion, we modeled graphite under hydrostatic pressure and uniaxial compression along c-axis. we calculated the phonon frequencies of the CM and LBM at various interlayer distance. We quantified the relationship between these two, separately by γ derived from the original definition of the Grüneisen parameter and by γ , following a phenomenological method. We found that the former method describes the relationship better and the in-plane strain makes a negligible contribution to the shifts of the interlayer modes. The γ can be further related to the interlayer coupling strength in the linear chain model and Eq. 1 and 4 can now be used to determine the interlayer strain in multilayer graphene. This strain determination method can also be applied to other layer materials in both bulk and multilayer forms.