# Temperature Dependence of Phonon Energies and Lifetimes in Single-and Few-layered Graphene

Markos Poulos, Konstantinos Papagelis, Emmanuel N Koukaras, George Kalosakas, Giorgia Fugallo, Konstantinos Termentzidis

► To cite this version:

Markos Poulos, Konstantinos Papagelis, Emmanuel N Koukaras, George Kalosakas, Giorgia Fugallo, et al.. Temperature Dependence of Phonon Energies and Lifetimes in Single-and Few-layered Graphene. Physical Review B, 2024, 110 (7), pp.075434. [arXiv:0477752v1](https://arxiv.org/abs/040777752v1)

HAL Id: hal-04777752

https://hal.science/hal-04777752v1

Submitted on 12 Nov 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Temperature Dependence of Phonon Energies and Lifetimes in Single- and Few-layered Graphene

Markos Poulos,¹, * Konstantinos Papagelis,²,³ Emmanuel N. Koukaras,³,⁴

George Kalosakas,⁵ Giorgia Fugallo,⁶ and Konstantinos Termentzidis⁷

¹INSA Lyon, CNRS, CETHIL, UMR5008, 69621 Villeurbanne, France

²Aristotle University of Thessaloniki, School of Physics, 54124 Thessaloniki, Greece

³Institute of Chemical Engineering Sciences, FORTH/ICE-HT, 26504 Patras, Greece

⁴Aristotle University of Thessaloniki, School of Chemistry, 54124 Thessaloniki, Greece

⁵University of Patras, Department of Materials Science, GR-26504 Rio, Greece

⁶CNRS, Université de Nantes, LTEN, UMR6607, 44306 Nantes, France

⁷CNRS, INSA Lyon, CETHIL, UMR5008, 69621 Villeurbanne, France

In this work, we have studied the phonon properties of multi-layered graphene with the use of Molecular Dynamics (MD) simulations and the k-space Autocorrelation Sequence (k-VACS) method. We calculate the phonon dispersion curves, densities of states and lifetimes $\tau$ of few-layered graphene of 1-5 layers and graphite. $\Gamma$-point phonon energies and lifetimes are investigated for different temperatures ranging from 80 K to 1000 K. The study focuses on the impact of the interlayer interaction and temperature on the energies and lifetimes of the $\Gamma$-point phonons, as well as the type of interlayer potential used. For the later we used the Kolmogorov-Crespi (KC) and the Lennard-Jones (LJ) potentials. We have found that the number of layers $N$ has little effect on the intra-layer (ZO and G) mode energies and greater effect on the inter-layer (Layer Shearing and Layer Breathing) modes, while $\tau$ is generally affected by $N$ for all modes, except for the Layer Shear mode. The trend of $N$ on the lifetimes was also found to be independent of the type of potential used. For the Raman-active G phonon, our calculations show that the lifetime increase with $N$ and that this increase is directly connected to the strength of the interlayer coupling and how this is modelled.

Keywords: Few-layered graphene; Phonons; Phonon lifetimes; Molecular Dynamics; kVACS; Temperature-Dependence; Interlayer Interaction

## I. INTRODUCTION:

Graphene since its discovery has attracted widespread attention due to its exceptional properties [1], especially for its excellent electronic conductivity [2], its performance in optoelectronic devices [3], its mechanical strength [4] and its high thermal conductivity [5]. This increased interest arised from the fact that high-quality monolayer graphene (1L) can be produced by means of mechanical exfoliation of Bulk Graphite (BG), as per the famous scotch-tape experiment [6]. There has also been recent interest in few-layered graphene (FLG), as it also exhibits interesting properties, ranging from high electrical [7] and thermal conductivity [8–10], to enhanced lubricity [11, 12] which is makes it an ideal filler material in composites for thermal management applications [13], while most recently, unconventional superconductivity was discovered in twisted bilayer graphene [14]. All the above properties are strongly phonon-related, as electrical conductivity is hypothesized to be phonon-limited [15], while the thermal conductivity of graphene is almost exclusively due to phonons [16, 17]. It is therefore of interest to map the change of the phonon properties of graphene as a function of the number of layers $N$, as well as of temperature.

In monolayer graphene (1LG), besides the 3 acoustic modes (longitudinal LA, transverse TA and out-of-

plane ZA) there exist 2 optical modes, the out-of-plane ZO and the in-plane LO and TO modes which at the $\Gamma$-point become the doubly-degenerate $E_{2g}$ mode, also known as the G-mode due to the Raman G-band that it is responsible for. The presence of multiple layers modifies the phonon structure of graphene, as new modes appear: existing intralayer optical modes split into Davydov multiplets, and new interlayer (IL) modes appear due to the folding of the Brillouin zone edge acoustic phonons [18, 19]. These emergent interlayer modes are categorised as layer Shear and layer Breathing modes, where the layers move as units in the in-plane and the out-of plane direction, respectively, with a phase difference between each other. The Layer Breathing mode (LBM) can also be found noted as ZO', while the Shear mode is also denoted as C-mode, again due to the Raman C-band signal.

Davydov splitting, also known as factor group splitting, refers to the removal of degeneracy of internal degrees of freedom due to weak external interactions, in materials whose unit cell is composed of replicated units, such as molecular crystals and layered Van der Waals (VdW) materials [20, 21]. The resulting modes usually have weakly split energies ($\Delta\omega < 5 \text{ cm}^{-1}$) and not far from the unperturbed system [22], although in cases of strong interlayer interaction such as in layered 2D antimonene, large splits ($\sim 70 \text{ cm}^{-1}$) have been reported [23]. As a result, in N-layered graphene there exist (N-1) Shear, (N-1) Breathing, N ZO and N G phonons. For a more detailed analysis of the phonon modes of few-layered graphene, we refer the interested

<footnote>* Corresponding Author: markos.poulos@insa-lyon.fr</footnote>
readers to the work of Saha et al. [24].

The dependence on temperature *T* of the phonon energies ω(*T*) and lifetimes τ(*T*) mainly arises due to anharmonicity in the interatomic interactions. As has been long established by Maraduddin and Fein (1962), by evaluating the phonon self-energy Σ(ω) of an interacting phonon system [25], one can obtain from its real Δ(*T*) and imaginary part Γ(*T*), respectively, the anharmonic phonon frequency shift and finite linewidth in the energies ω₀ of the harmonic non-interacting phonons. The T-dependence then mainly arises due to the dependence of Γ and Δ on the phonon populations n(ω, *T*). Lattice thermal expansion also contributes to Σ(ω) by the implicit inclusion of higher order anharmonic terms in the renormalization of the harmonic energy [26] leading to a temperature-dependent ω₀(*T*) since at higher temperatures these terms can no more be considered as small perturbations [27]. However, if one can determine the lattice thermal expansion independently (e.g. experimentally), the shift at a given temperature can be obtained computing the harmonic frequency at the lattice parameters corresponding to that temperature.

The phonon energies of layered graphene and graphite have been experimentally probed quite extensively in the literature, mostly concerning the zone-center Raman-active G phonon (ωG ≈ 1580 cm⁻¹). For graphite, the phonon dispersion curves by means of In-elastic X-ray Scattering (IXS) have been reported [28], and for 1L graphene studies on the T-dependence of ωG(T) by means of conventional Raman spectroscopy are also present in the literature, with excellent agreement between them [29–31]. Studies of the effect of the number of layers *N* on the phonon energies of FLG have also been conducted [32–36]

For the experimental study of the phonon lifetimes of 1L graphene, conventional Raman spectroscopy is usually performed on suspended samples in order to exclude the effect of the substrate [31, 37]. The lifetime of a Raman-active mode can be estimated by the linewidth γ in cm⁻¹ as τ = (2γπc)⁻¹ [38, 39]. The major contributions to the linewidth are due to el-ph and ph-ph interactions γG = γG<sup>ph-ph</sup> + γG<sup>el-ph</sup>. Most notably for the G-mode, the major contribution to γG is due to the electron-phonon part [31, 40, 41]. For the G-mode, the ph-ph contribution can be also probed separately by means of Time-Resolved Incoherent Anti-Stokes Raman Spectroscopy (TRIARS). It appears commonly accepted that for graphite τG<sup>ph-ph</sup> ranges between 2-2.5 ps at room temperature [40, 42]. For 1L, however, reported values are dispersed, from ~1-1.5 ps [38, 39] to 2.5 ps [41] for supported samples, and as high as 4.9 ps for samples supported on colloidal solutions [43]. The growth method and the sample quality have a big impact in the measured value. An overview of the different measured values reported in the literature can be found in Katsiaounis et al [38].

Phonon energies and lifetimes are usually studied theoretically by means of *ab initio* methods based on Perturbation Theory. Initially calcula-

tions in the literature included only the electron-phonon (e-ph) and three-phonon (3ph) contributions to the temperature-dependence of τ<sub>G</sub>(T) in single-layer graphene [44]. These calculations showed that the e-ph contribution overwhelmingly dominates the temperature dependence, causing γ<sub>G</sub> to decrease with temperature. This finding starkly contrasts with all the experimental works cited above. The discrepancy was only recently resolved. It has been reported that the 4ph contribution is actually much more significant than 3ph, although the e-ph contribution still remains dominant [31]. The same authors reported that phonon renormalization due to lattice thermal relaxation was crucial to avoid overestimating the 4-ph scattering rates [31].

An alternative way to calculate phonon properties without the above limitations is to employ the k-space Velocity Autocorrelation Sequence method within classical Molecular Dynamics (MD) simulations, which we have previously employed to obtain the temperature dependence of the Γ-point phonon energies in SLG [45]. With this method we can obtain phonon energies and lifetimes for all *k*-points in the Brillouin Zone allowed by the Periodic Boundary Conditions [45, 46]. For each *k*-point the Spectral Energy Density is obtained, which is then fitted by Lorentzian functions. The positions and widths of the peaks correspond to the phonon energies and linewidths respectively. The Phonon Density of State (PDOS) can be also extracted within the same scheme and from the same data sets. The two main advantages of this method is that it inherently takes into account by default all orders of anharmonicity and lattice expansion, while it is also applicable to both periodic and non-periodic systems. The method has been described extensively in a previous works and readers interested in the details are referred to it [45]. The further details and a theoretical proof of the method are also provided in Appendices B and C.

kVACS has been extensively used in the past to study vibrational properties in condensed matter applications, from liquids [47] to molecules [48] and the temperature dependence of phonon modes in Si [49], to surface modes in Cu [50, 51]. In more recent studies, it has been applied to obtain the phonon spectra of carbon nanotubes [52], monolayer graphene [53] and nanoscale Si [54]. Recently it was also used to investigate the contribution of the phonon vibrational structure in the mechanism of thermal rectification of certain nanostructures, as in the example of asymmetric two-phase Si-nanowires [55], and in MoSe₂-WSe₂ [56] and graphene-hBN Van der Waals lateral heterojunctions [57].

In layered materials, adjacent layers interact with each other through weak Van der Waals (VdW) forces arising from an attractive spontaneous dipole-dipole dispersion interaction [58] and a repulsive interaction arising from the Pauli repulsion of orbital overlap [59]. The attractive part is long-ranged, while the repulsive part is short-ranged, and each term is dominant in their respective distance range [59]. In MD simulations this in-
teraction has been modelled quite extensively with the classical Lennard-Jones (LJ) potential, however, due to its being a two-body radial potential, it has been found to be too smooth to describe corrugation and variations in the stacking alignment of layered graphitic materials [60]. For this reason, a new type of interlayer potential for graphitic materials was proposed by Kolmogorov and Crespi (herby denoted as KC), which is a registry-dependent potential that takes specifically into account the directionality of the repulsive overlap interaction between π-orbitals of adjacent layers [60]. This is the main interlayer potential employed in this work. Other registry-dependent potentials have also been developed ever since, such as the ILP potential which includes Coulomb interactions for polar layered materials like hBN [61], the DRIP potential which includes dihedral interactions for twisted bilayer graphene [62] and the Lebedeva potential which includes a term dependent on the vertical distance of the $p_z$ orbitals of adjacent atoms, to better describe the compressibility of graphite [63].

The rest of the article is organised as follows: First, a synopsis of the computational methods and interatomic potentials that we used in this paper is made, followed by the discussion of the results. There, we present the calculated phonon dispersion curves and Density of States for Bulk Graphite, and the temperature dependence of the energies and lifetimes of Γ-point phonons 1-2-3-5 layered graphene and graphite, calculated with the KC potential and compared with all available experimental data. The calculations of energies and lifetimes are repeated using the LJ potential and the trends are compared to the results with KC. Finally the effect of the interlayer coupling strength of the LJ potential on the lifetime of the G phonon is investigated for 2-layered graphene and graphite.

## II. METHODOLOGY:

In this work we have performed MD simulations on single- and few-layered pristine graphene (FLG) and graphite (BG), using the kVACS method. The phonon properties of the above systems were studied under temperatures ranging from 80 K to 1000 K in order to obtain the temperature dependence of the phonon energies and lifetimes.

The interatomic C-C interactions within the same layer were modelled with the reparametrised Tersoff potential by Broido et al [66]. It shall be pointed out that in a previous work the LCBOP and the AIREBO potentials were shown to give a very poor temperature dependence of the energy of the G phonon in 1L graphene [45]. A Machine-Learning based Gaussian Approximation Potential for graphene has also been reported to give very promising results and reproduced very well the DFT dispersion curves and thermal expansion coefficients [67]. However the increased computational cost of such potentials was prohibitive for the

purposes of our work, which mainly focus on probing the interlayer interactions, and which are independent on the intralayer potential used. For the above reasons, the revised Tersoff potential was deemed the best solution. For the interlayer interactions we used the reparametrised full Kolmogorov-Crespi (KC) potential with a taper function as described in [68] and whose form is given in Appendix A. In order to investigate the impact of the nature of the interlayer potential and the effect of the coupling strength, we repeated part of the calculations using the classical Lennard-Jones (LJ) potential for comparison, given by

$$V_{LJ}(r) = 4\epsilon \left[ \left(\frac{\sigma}{r}\right)^{12} - \left(\frac{\sigma}{r}\right)^{6} \right] \quad (1)$$

For the KC potential, a global cutoff of $r_c = 12$ Å was chosen. For the LJ potential, the parameters $\epsilon= 4.6$ meV and $\sigma = 3.276$ Å used were taken from Lindsay et al [18]. In all of our calculations, all layers can interact with all the others, which mostly happens via the long-range part of the VdW interaction, which "passes through" layers. As it has been underlined in previous works, in few-layered graphene the cleavage energy (CE) (the energy required to split a bulk material in half along a basal plane) and the interlayer Binding Energy (BE) (the energy required to separate a bulk material into individual non-interacting layers) is different, which is a sign of interlayer interactions longer than between adjacent layers [69]. Experimental measurements of CE and the BE of graphite support this claim, as their difference was found to be ~15% [70].

The FLG systems studied were 1L, 2L, 3L, 5L graphene and BG, all in Bernal (-AB-) stacking. The computational cell used consisted of 40x40 primitive cells with Periodic Boundary Conditions (PBC), having sides of 1 nm. Each system was first structurally and thermally relaxed for every temperature step: an initial conjugate gradient (CG) static energy minimization followed by a 125 ps NVE initial equilibration, a 250 ps NPT structural relaxation and a final 125 ps NVT thermal equilibration.

During the production step to obtain the kVACS spectra, 10 independent equilibrium NVE runs with randomized initial velocities were performed. This statistical sampling was done in order to obtain good ensemble averages and to ensure ergodicity. For the simulations using KC as the interlayer potential, each run consisted of 655.360 timesteps or 327.68 ps, for a total of 3.28 ns of production time per studied case. The sampling rate of the velocity data was every 10 timesteps or 5 fs, giving an $\omega$-resolution of 1.1 cm$^{-1}$. Comparative convergence tests with higher resolution were performed for the 2L case and are presented in Section 2 of the Supplemental Material [71]. No significant deviation of the calculated $\omega$ and $\tau$ with resolution were found. For the case of the LJ potential, convergence tests showed that this resolution gave poor results. For this reason, for all the LJ results presented in this paper the per-run simulation time was 655.36 ps, double that of KC, and the sampling rate
<figure><img src="image_2.png" /><figcaption>FIG. 1. (Left) The calculated phonon dispersion curves of Bulk Graphite (full symbols) at $T = 300$ K, together with data from IXS measurements (open symbols) [28, 64]. The magnitude of the kVACS symbols is a measure of the calculated lifetime. Straight lines connecting points are b-splines serving as visual guides. (Right) The calculated phonon DOS for Bulk Graphite at $T = 300$ K (straight lines), along with data from INS measurements [65] (open symbols). All PDOS are normalized to unit area.</figcaption></figure>

was every 15 steps, giving an $\omega$-resolution of 0.5 cm$^{-1}$. The number and frequency of sampled data points was eventually the result of a compromise between computational efficiency and sufficient $\omega$-resolution. For all simulations, the timestep used was 0.5 fs. All MD simulations were performed using the LAMMPS code [72].

The phonon properties were calculated from the kVACS method and the spectra obtained were fitted using simple Lorentzian functions. The Brillouin Zone was sampled along the path $\Gamma$-K-M-$\Gamma$. The Lattice Dynamics (LD) simulations were also performed using LAMMPS's dynamical matrix command to calculate and output the force constant matrix of the systems under study. The dynamical matrices for specific k-points were calculated by performing a spatial Fourier transform over the equilibrium positions and finally the phonon energies were obtained by the diagonalization of the dynamical matrix.

Finally, an important point should be raised, concerning the resolution of the k-points examined. As we are using finite computational cells and PBC, not any arbitrary k- can be sampled, but only those that are linear combinations of $\frac{1}{N_1}\mathbf{b}_1$ and $\frac{1}{N_2}\mathbf{b}_2$, where $\mathbf{b}_1$ and $\mathbf{b}_2$ are the primitive inverse lattice vectors of the respective unit cell and $N_1, N_2$ the number of unit cells along each lattice vector, respectively. Here, $N_1 = N_2 = 40$. The procedure to obtain the correct k-point resolution has been extensively documented in a previous study [45] and for further details we refer the reader to it.

### III. RESULTS AND DISCUSSION:

For all multi-layered systems, as has been mentioned in the Introduction, all optical modes (Shear, LBM, ZO and G) appear in multiplets, therefore a choice of representation must be made. In the rest of this work, the highest-energy ZO Davydov multiplet is chosen to represent the said mode, which corresponds to a vibration where the crystallographically equivalent C atoms of each layer vibrate in phase. For the G mode, the lowest-energy multiplet is chosen, which also corresponds to the in-phase vibration of the equivalent C atoms. This guarantees that the G-mode chosen will be Raman active, independent of the number of layers [24]. Concerning the interlayer modes, both for the LBM and the C-mode branches for each N the highest-lying branch corresponds to an anti-phase layer vibration, while the lowest-energy branches correspond to phonons where mostly the outer layers vibrate [36]. For reasons of consistency with the experimental studies, we will represent the shear modes by the highest-energy branch (anti-phase), while for the breathing mode we choose the lowest-lying one, where layers oscillate like an accordion. The reason for this choice is that these modes represent an actual layer shearing/breathing vibrational pattern, while also being always Raman-active and giving the highest intensity [34]. An important remark should be made about the LBM of graphite, which doesn't follow the aforementioned vibrational pattern, as all layers vibrate out-of-phase. A visual illustration of the vi-
ibrational patterns of all the modes mentioned above is provided in Figs. SM-1 and SM-2 of the Supplemental Material [71].

From a first comparison of the results, it turns out that the number of layers *N* affects differently the intralayer (acoustic, ZO and G) and the interlayer (C and LBM) modes across the entire Brillouin zone and throughout the entire temperature range. Specifically for the intralayer modes, the calculated energies increase slightly going from from 1L to FLG and BG in the order of 1-5 cm⁻¹, therefore no visible difference was seen between the equivalent phonon dispersion curves for different *N* and different T, especially at the energy scale covering all the mode energies (0-1700 cm⁻¹). For this reason, the dispersion curves of graphite were thus chosen to represent the intralayer curves of FLG for all *N*, which we present in Fig. 1, along with information about the calculated lifetimes, in comparison with data from Inelastic X-Ray Scattering (IXS) measurements [28, 64]. We also provide the calculated Density of States (PDOS) of graphite, compared with IXS measurements from a different author [65], all at T=300 K. The interlayer modes are affected in a more intricate way and further details will be provided later in the text. A comparison of all phonon branches for all *N* is also provided in Fig. SM-8 of the Supplemental Material [71].

We directly observe from Fig. 1 that the acoustic modes are described with very good accuracy by the Tersoff-2010 potential with respect to experimental data, as was the case in the initial work of potential reparametrization by Broido et al [66], where the dispersion curves were obtained by simple LD. Concerning the calculated lifetimes, they fall within the range of ~1-20 ps and two principal remarks are to be made. First of all, the out-of-plane intralayer modes (ZA and ZO) have noticeably much larger lifetimes than the in-plane ones (TA, LA, LO, TO), ranging between ~5-20 ps and increasing close to the Γ-point. Interestingly, the ZA and ZO mode lifetimes have close values. This is in stark contrast with the in-plane modes, where the acoustic ones globally have larger lifetimes (5-10 ps) than optical ones (1-5 ps).

The energies of the LO/TO branches are overestimated by ~10%, while the ZO branch is significantly overestimated at about 45%. Also, the description of the LO and TO branches is satisfactory on the Γ-M path (armchair), but significantly overestimates about 27% the splitting of the LO branch at the ΓK (zig-zag) and KM (parallel to the C-C bond axis) paths. Also, the TO branch has higher energies than the LO branch, which has also been observed in [45] and can be attributed to the Tersoff potential. It should be further noted that the discrepancy in the energies the optical modes shall not be attributed to the kVACS method, since even in the initial work, the dispersion curves for the optical modes, obtained by Lattice Dynamics, exhibited the same level of discrepancy with experiment [66]. Besides, as it will be later shown, the relative difference between LD and

kVACS energies of the ZO and G phonons is of the order of 1%. The original authors attributed the failure of the simultaneous fitting of both optical and acoustic modes to the short range of the potential [66].

The PDOS calculated with kVACS is compared with Inelastic Neutron Scattering (INS) measurements for graphite [65], both presented in Fig. 1 (right). Both calculated and experimental PDOS are normalized to unit area. We attribute the differences with experimental data to the optical modes, since from the comparison of the Dispersion Curves in Fig. 1 (left), the acoustic modes are perfectly matching the experimental data. Also, in our calculations no account of the out-of-plane direction (Γ-A) of the Brillouin Zone of graphite was taken, those modes are therefore missing from the PDOS.

<figure><img src="image_3.png" /><figcaption>FIG. 2. The calculated dependence of the FLG (a) Shear and (b) Breathing mode at T = 300 K on the number of layers *N* (full symbols), along with the fitting of the energies as given by eqs. (3) and (4) (straight line). Experimental measurements by Tan et al. [34] and Lui et al. [36] respectively are also shown for comparison (hollow symbols).</figcaption></figure>

It is also of interest to investigate the effect of *N* on the energy of the interlayer phonon modes. For this, in Figs 2a) and 2b) we plot the energy of the true Shear and Breathing modes at T=300 K, as a function of the *N*. As it can be clearly seen, the two modes have opposite trends, as the C-mode hardens with increasing *N*, while the LBM softens, as is also attested by the exper-
imental measurements. This trend can be explained by a simple linear chain model, where layers are modelled by single particles and the interlayer mode frequencies are given by [73]

$$ \omega_{\nu}^{(i,N)} = \sqrt{\frac{2k_{\nu}}{\mu} \left\{ 1 - \cos\left(\frac{(i-1)\pi}{N}\right) \right\}} \quad (2) $$

where N is the number of layers present, $\nu = 1, 2, 3$ represents the interlayer mode branch ($\nu = 1, 2$ the doubly-degenerate C-mode and $\nu = 3$ the LBM), of N modes each, indexed by $i = 1, 2, 3...N$ ($i = 1$ refers to rigid translation). Also, k is a measure of the force constant per unit area that corresponds to each vibration (shear force for the C-mode and tensile force along the c-axis for the LBM) and $\mu$ is the mass per unit area of a single layer, taken equal to $\mu = 7.6 \times 10^{-27}$ kg·Å$^{-2}$. The model has been successful to capture the trends in experimental data for the Shear [36] and Breathing modes [34], by simply setting $i = N$ and $i = 2$ in eq. (2) respectively. After some manipulation, these energies, measured in cm$^{-1}$ become

$$ \omega_C(N) = \frac{1}{\pi c} \sqrt{\frac{\alpha}{\mu}} \cos\left(\frac{\pi}{2N}\right) \quad (3) $$

$$ \omega_{LBM}(N) = \frac{1}{\pi c} \sqrt{\frac{b}{\mu}} \sin\left(\frac{\pi}{2N}\right) \quad (4) $$

In eqs. (3) and (4), $\alpha$ and $b$ are different symbols used to represent $k_\nu$ for the shear and tensile force constants respectively, and $c = 2.997 \times 10^{10}$ cm·s$^{-1}$ is the speed of light. It should also be noted that in eqs. (3) and (4) the pre-factors of the sinusoids stand for the C-mode and LBM energies of graphite, respectively. Again for the LBM of graphite, this energy does not derive from $\lim_{N \to \infty} \omega_{LBM}(N)$ of eq. (4), since the former refers to an anti-phase vibration and not the “accordion-like” motion of the true LBM (see Fig. SM-1 of the Supplementary Material [71]), which for $N \to \infty$ tends to zero.

From the comparison with the model and the experimental data, we can see that the kVACS simulations follow the model in an excellent way. Namely, for graphite $\omega_C^{calc} = 37.8$ cm$^{-1}$ and $\omega_C^{exp} = 43.4$ cm$^{-1}$, while the fitting of the C-mode energies using eq. (3) gives a shear force of $\alpha_{calc} = 0.96 \times 10^{19}$ N·m$^{-2}$ for the interlayer interaction, which is very close to the value of $\alpha_{exp} = 1.27 \times 10^{19}$ N·m$^{-2}$ obtained from the experimental fitting [34]. Accordingly, the trend of the LBM energies follows very closely the model of eq. (4), as it predicts an energy of $\omega_{LBM}^{calc} = 108.8$ cm$^{-1}$ for graphite which is very close to the experimental value $\omega_{LBM}^{exp} = 114.6$ cm$^{-1}$. Again by using eq. (4), we obtain $b_{calc} = 7.97 \times 10^{19}$ N·m$^{-2}$, while $b_{exp} = 8.84 \times 10^{19}$ N·m$^{-2}$. Interestingly, if we calculate the same interlayer

phonon energies for graphite using Lattice Dynamics with the same potential parameters, which completely ignores thermal motion and anharmonicity, we obtain $\omega_C^{LD} = 42$ cm$^{-1}$ and $\omega_{LBM}^{LD} = 123$ cm$^{-1}$, which give $\alpha_{LD} = 1.19 \times 10^{19}$ N·m$^{-2}$ and $b_{LD} = 10.18 \times 10^{19}$ N·m$^{-2}$.

From the above discussion we can extract two useful conclusions. First, even though the KC potential is a complex multi-parameter potential, and despite the inclusion of thermal motion and anharmonicity within kVACS, the behaviour of the interlayer phonon energies eventually follows a simple linear chain model, with a single adjustable parameter, the interlayer force constants $k_\nu$. Secondly, because of the latter argument, a significant part of the discrepancy with experiment can be tracked down to the parametrization of the KC potential. This is also supported by the LD calculations, which, especially for the case of the LBM, reveal that thermal motion and anharmonicity can significantly reduce the effective interlayer force constants and phonon energies. This is more or less expected, as the parametrization of the KC potential used in this work was done via DFT, a zero-temperature method [68].

In the following, we obtained the temperature dependence of phonon energies and lifetimes for the Shear, Breathing, ZO and G modes of FLG. For temperatures higher than the mode-specific Debye temperature $\Theta_D = \hbar\omega/k_B$, where $\omega$ is the harmonic energy of the phonon in question, $\omega(T)$ is linear or quadratic (classical limit), depending on whether the dominant anharmonic process is the 3- or 4-phonon [77]. As MD obeys classical statistics, though, this behaviour is expected to be present even at low temperatures, due to the Equipartition of energy being valid at all temperatures in Boltzmann statistics, but only at temperatures $T \gg \Theta_D$ for Bose statistics. Therefore the slopes $\chi_T = \frac{d\omega}{dT}$ are a characteristic measure of phonon anharmonicity. Another possibly interesting metric that can quantify the effect of the thermal motion on the phonon softening could be the relative difference between the LD energy value and the lowest-temperature kVACS value, that is $\delta\omega_{LD} = (\omega_{LD} - \omega_{T=80K})/\omega_{T=80K}$. The slopes $\chi_T$ for all $\Gamma$-point phonon modes and numbers of layers are presented in Table I.

<table><caption>TABLE I. The energy slopes with temperature $\chi_T$ for each $\Gamma$-point phonon mode and for various numbers of layers calculated via kVACS. Experimental slopes are taken from [75, 76]. All values are given in cm$^{-1}$K$^{-1}$.</caption><thead><tr><th># of Layers</th><th>Shear</th><th>LBM<br>(×10<sup>-3</sup> cm<sup>-1</sup>K<sup>-1</sup>)</th><th>ZO</th><th>G</th></tr></thead><tbody><tr><td>1</td><td>-</td><td>-</td><td>-36</td><td>-51</td></tr><tr><td>2</td><td>-6.2</td><td>-9.7</td><td>-35</td><td>-47</td></tr><tr><td>3</td><td>-5.9</td><td>-6.3</td><td>-34</td><td>-48</td></tr><tr><td>5</td><td>-5.7</td><td>-3.6</td><td>-35</td><td>-49</td></tr><tr><td>BG</td><td>-4.4</td><td>-5.6</td><td>-35</td><td>-50</td></tr><tr><td>Exper.</td><td>-</td><td>-</td><td>-</td><td>-21</td></tr></tbody></table>
<figure><img src="image_4.png" /><figcaption>FIG. 3. The calculated temperature dependence of the true FLG (a) Shear (b) Breathing (c) ZO and (d) G phonon mode energies with the number of layers *N*. The correspondence of symbols with *N* is as follows. 2L: green squares, 3L: red diamonds, 5L: blue stars, BG: black circles. The hollow symbols represent experimental data by Cong et al [74] Yang et al [75] and Liu et al [76]. The straight lines with the same color correspond to the least-square fittings of the associated data. The dashed lines represent the LD energies for each *N*.</figcaption></figure>

In Figs. 3a) and 3b) we present the calculated T-dependence of the interlayer Shear and Breathing mode energies for FLG, along with linear fittings and the corresponding energies obtained from Lattice Dynamics. Experimental measurements for the 2LG Shear mode by Cong et al. [74] are also presented. Even down to T=80 K, the calculated $\omega(T)$ appears linear, which of course contrasts with the real behaviour, as quantum phenomena are absent in MD and more prevalent at low T. Also, with increasing *N*, the slopes have distinctly decreasing values for the interlayer modes. For the Shear mode, they go from with $-6.2 \cdot 10^{-3} \text{ cm}^{-1} \text{ T}^{-1}$ for bilayer, down to $-4.4 \cdot 10^{-3} \text{ cm}^{-1} \text{ T}^{-1}$ for graphite. The same trend is observed for the LBM, where the slopes are ranging from $-9.7 \cdot 10^{-3} \text{ cm}^{-1} \text{ T}^{-1}$ for bilayer, down to $-3.6 \cdot 10^{-3} \text{ cm}^{-1} \text{ T}^{-1}$ for pentalayer. Again for the LBM of graphite, neither the energies nor the slope follow the above trend, as the energies are higher than those of bilayer, while the slope is between the values of 3ML and 5ML, at $-5.6 \cdot 10^{-3} \text{ cm}^{-1} \text{ T}^{-1}$. For both interlayer modes, $\delta\omega_{LD} \approx 5-10 \%$. To our knowledge, no experimental work on the temperature dependence of the energies of the LBM of FLG has been published so far.

A completely different picture is observed for the intralayer modes, where both the energies and the slopes are almost insensitive to *N*. In Figs. 3c) and 3d) we present the temperature dependence of the intralayer

ZO and G mode phonon energies for FLG's with different *N*. The energies are slightly upshifted by ~1-5 cm⁻¹ going from 1L to graphite. Experimental measurements on the G mode also confirm these findings [32–34]. The LD energy is not affected by *N* either. The slopes $\chi_T$ are also practically constant with *N* and one order of magnitude larger than for the interlayer modes (see Table I). Interestingly it is found that for both intralayer modes $\delta\omega_{LD} \approx 0.3 \%$, which is an order of magnitude lower than $\delta\omega_{LD}$ for the interlayer modes, although on the contrary, the slope $\chi_T$ of the intralayer modes is an order of magnitude higher than $\chi_T$ for the interlayer modes (see Table I).

By comparing with the experimental measurements, we notice again that the Tersoff-2010 potential significantly overestimates the G-mode energies by a factor of ~5.5 %, and the slope $\chi_T$ by a factor of ~2.5 [75, 76]. Furthermore, the experimental measurements of Yang et al in multilayered graphene demonstrated that going from monolayer to FLG's, $\omega_G$ is slightly reduced [75], in contrast to our findings. However, all the measurements in their work were performed on 1ML supported on a Si-substrate. Whether the discrepancy is due to the limitations of the Tersoff potential or to the presence of the substrate in the experimental measurements cannot be determined by this comparison.

Finally the phonon lifetimes due to anharmonicity $\tau^{ph-ph}$ were determined from the linewidth $\Gamma^{ph-ph}$ of
<figure><img src="image_5.png" /><figcaption>FIG. 4. The calculated temperature dependence of $\tau^{ph-ph}$ of the FLG (a) Shear, (b) Breathing, (c) ZO and (d) G modes with the number of layers $N$. The correspondence of symbols with $N$ is as follows. 2L: green squares, 3L: red diamonds, 5L: blue stars, BG: black circles. For the G mode, the experimental data from Chatzakis et al [40] are also presented (hollow circles). Dashed lines connecting points are b-splines as visual guides</figcaption></figure>

the kVACS peaks (FWHM), by the simple relation

$$ \tau^{ph-ph} = \frac{\hbar}{\Gamma^{ph-ph}} \qquad (5) $$

For convenience with units, we use $\hbar = 5.29 \text{ ps} \cdot \text{cm}^{-1}$. In Fig. 4a) and 4b), we present the calculated temperature dependence of lifetimes of the interlayer Shear and Breathing modes respectively, for various numbers of layers. Concerning the Shear mode, it appears that within the uncertainties of the peak fitting, $\tau^{ph-ph}$ is practically insensitive to $N$. This is in agreement with some reported DFT calculations and experimental measurements, where it is suggested that the e-ph contribution in $\tau_C$ is independent of the number of layers, and the measured lifetimes changed very little with $N$ [34]. The same authors have claimed that the very low C/G peak Raman intensity ratio and their DFT calculations also indicate that the e-ph contribution to the C-band linewidth is negligible and independent of $N$ therefore the observed Raman linewidth could be directly associated with $\gamma_C^{ph-ph}$ [34]. Other authors have calculated a dominant e-ph contribution and the measured C-band linewidth appears decreasing with temperature [74], which is a sign of e-ph dominated linewidth. On the contrary, for the LBM $\tau^{ph-ph}$ increases with increasing $N$, with the exception of graphite, which appears to have values very close to 2LG. Again, the graphite LBM is an anti-phase mode and its eigenvectors resemble more those of the bilayer

LBM, which could also explain the similar $\tau^{ph-ph}(T)$ trends of the two systems. This could also explain the increase of $\tau_{LBM}$ with increasing $N$, since the eigenvectors for the true LBM are mostly zero for the inner layers and mostly the outer layers vibrate. For graphite all layers vibrate, leading to more scattering of the LBM.

Finally, the calculated lifetimes of the intralayer modes are presented in Fig. 4c) and 4d) for ZO and G, respectively. Interestingly, for both modes the lifetimes increase with increasing $N$ along the entire temperature range, with 1LG having $\tau_G^{ph-ph}=1.9 \text{ ps}$ while the calculated $\tau_G^{ph-ph}=5.6 \text{ ps}$, more than double that of the monolayer. Although to our knowledge no experimental study has so far directly measured the $\tau_G^{ph-ph}(T)$ of high-quality, monocrystalline suspended 1L graphene, one would expect intuitively that the weak interlayer interaction in FLG should not significantly affect phonon lifetimes, compared to the monolayer.

Experimentally, the total $\gamma_G$ of suspended 1L graphene at T= 300 K has been measured at 16 cm$^{-1}$ [37], giving an estimated value of the total lifetime $\tau_G= 0.33 \text{ ps}$. For graphite, $\tau_G^{ph-ph}(T)$ has been experimentally measured at around ~2.5 ps by multiple authors [40, 42], which is almost 2 times lower than our calculated values with the KC potential (~5.6 ps). For supported 1L samples, values of ~1-1.5 ps have been reported, increasing with increasing the number of layers $N$ [38, 39], while some authors have measured
$\tau_G^{ph-ph}$-2.5 ps, equal to that of graphite [41]. Values for 1L suspended on colloidal solutions as high as 4.9 ps for $\tau_G^{ph-ph}$ have also been reported [43].

<figure><img src="image_6.png" /><figcaption>FIG. 5. The calculated temperature dependence of $\tau_G^{ph-ph}$ of the monolayer Raman-active G phonon (black triangles). DFT calculations including 3ph+4ph processes and phonon renormalization [31] are also shown (crossed blue triangles). Calculations are compared with experimental data from TRIARS measurements of supported 1L graphene [39] (empty red triangles). Dashed lines connecting points are b-splines as visual guides</figcaption></figure>

We have also compared our calculations for $\tau_G^{ph-ph}$ with DFT calculations and some of the existing experimental data, shown in Fig. 5. The DFT calculations are for free 1LG and are extracted from the work of Han et al (2022) by using only the $\Gamma^{ph-ph}$ from the renormalized 3ph+4ph processes [31], and by inputing it in eq. (5). The experimental data are taken from TRIARS measurements of supported 1L graphene by Kang et al (2010) [39]. The kVACS and DFT curves agree well with each other, although the DFT absolute values are systematically ~ 40% lower. The experimental values exhibit a different trend, with $\tau_G^{ph-ph}$ remaining almost constant for a temperature $300K \le T \le 600K$. The authors have also reported that $\tau_G^{ph-ph}$ values for 1L are significantly lower than for graphite, and an increasing trend with increasing $N$, similar to the one observed in Fig. 4b). However, they have attributed this behaviour to the effect of the substrate [39].

In order to clarify the trend of $\tau^{ph-ph}(T)$, we have repeated part of the calculations replacing KC with the Lennard-Jones (LJ) potential for the interlayer interaction. Namely, the T-dependence of the $\Gamma$-point phonon energies and lifetimes for 2L and bulk graphite were repeated, for $80K \le T \le 900K$. The slopes $\chi_T$ were also calculated and are presented in Table II. To ensure that the trends for both potentials are within the linear regime, the least-square fitting was performed for $T > 200$ K for both potentials.

In Figs. 6a) and 6b) we present the comparison for the interlayer mode energies, where the trends of the phonon

<table><caption>TABLE II. The comparison of energy slopes with temperature $\chi_T$ for each phonon mode and for bilayer graphene and graphite, between the two interlayer potentials used (Lennard-Jones and Kolmogorov-Crespi)</caption><thead><tr><th rowspan="2">(×10<sup>-3</sup> cm<sup>-1</sup>K<sup>-1</sup>)</th><th colspan="2">KC</th><th colspan="2">LJ</th></tr><tr><th>Bilayer</th><th>Graphite</th><th>Bilayer</th><th>Graphite</th></tr></thead><tbody><tr><td>Shear</td><td>-6.2</td><td>-4.3</td><td>-0.4</td><td>-0.5</td></tr><tr><td>LBM</td><td>-9.7</td><td>-5.6</td><td>-11.8</td><td>-6.8</td></tr><tr><td>ZO</td><td>-35.4</td><td>-34.8</td><td>-35.8</td><td>-35.2</td></tr><tr><td>G</td><td>-47.9</td><td>-50.0</td><td>-52.7</td><td>-47.9</td></tr></tbody></table>

energies $\omega_{inter}(N)$ with increasing number of layers is as expected for both potentials (increasing for the Shear mode and decreasing for LBM). However the absolute values obtained with LJ are quite poor. The LBM energy is significantly overestimated (at $T = 300K$ the error is 30-35%, compared to only 0.5-5% for KC), while on the contrary the Shear mode is enormously underestimated (at $T = 300K$ the difference with experiment is 75%, compared to only 15% for KC), both for 2L and graphite. Also, the LJ gives for the Shear mode a slope $\chi_T$ that is an order of magnitude lower than the one by KC and, as this slope is already low for KC (~$10^{-3}$ cm$^{-1}$K$^{-1}$), LJ essentially predicts that the Shear mode energy is constant with temperature. The comparison of the slopes for all modes are presented in Table II.

As for the intralayer mode energies, presented in Figs. 6c) and 6d), the differences between the two potentials are much less pronounced, as the intralayer modes are little affected by the weak interlayer interactions. The slopes $\chi_T$ calculated with both potentials agree very well with each other, as can be clearly seen in Table II. A small but noticeable difference is that for the ZO mode, the difference in energies between 2L and graphite is much more important with LJ (~$3 \cdot 10^{-3}$) than with KC (~$2 \cdot 10^{-4}$). This and the overestimation of the LBM mode energies could possibly be related, indicating that the shear and perpendicular interlayer forces cannot be simultaneously fitted in a 2-parameter and 2-body potential like LJ, as the perpendicular components are clearly overestimated, while the shear components are underestimated.

Subsequently, we shall discuss the impact of the choice interlayer potentials on phonon lifetimes with respect to the number of layers $N$. In Figs. 7a) and 7b) we present this comparison for the interlayer mode lifetimes. It shall first of all be noted that the Shear mode SED peaks exhibit a poor Lorentzian shape and therefore the lifetimes $\tau_C^{ph-ph}(T)$ calculated with LJ have large uncertainties, but are nonetheless presented for reasons of consistency. The LJ results (Fig. 7a)), however, also demonstrate that for the Shear mode $\tau^{ph-ph}(T)$ is, within the uncertainty of the fitting, practically unaffected by $N$, as was also seen with KC. In terms of absolute values, however, the LJ potential systematically underestimates lifetimes by a factor of 2
<figure><img src="image_7.png" /><figcaption>FIG. 6. The calculated temperature dependence $\omega(T)$ of phonon energies for the interlayer (a) Shear (b) Breathing (c) ZO and (d) G mode in 2L (square symbols) and graphite (round symbols), using KC (black full symbols) and LJ (red crossed symbols) as interlayer potentials.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>FIG. 7. The calculated temperature dependence $\tau^{ph-ph}(T)$ of phonon lifetimes for the interlayer (a) Shear (b) Breathing (c) ZO and (d) G mode in 2L (square symbols) and graphite (round symbols), using KC (black full symbols) and LJ (red crossed symbols) as interlayer potentials. For the G-mode, the experimental data for graphite from Chatzakis et al [40] are also presented (hollow circles). Dashed lines connecting points are b-splines as visual guides.</figcaption></figure>

compared to KC. As for the LBM lifetimes (Fig. 7b), again LJ gives systematically underestimates $\tau_{LBM}^{ph-ph}(T)$
and predicts that for graphite and 2L the LBM lifetimes have practically the same value.

Finally, we shall focus on the comparison of the ZO and G mode lifetimes, presented in Figs. 7c) and 7d) respectively. Beginning with the ZO mode, one may first of all observe that $\tau_{ZO}^{ph-ph}(T)$ for 1L and 2L are practically equal, irrespective of the interlayer potential used. However, as for the case of the interlayer modes discussed above, LJ gives lower values than KC and predicts that the values for graphite are equal to those of 1L and 2L, in contrast to KC, which predicts that $\tau_{ZO}^{BG} > \tau_{ZO}^{1L}, \tau_{ZO}^{2L}$

Last, focusing on the G mode lifetimes, one can make the following observations: Both potentials give results that are close in terms of absolute value and with respect to the number of layers $N$, in general agreement with experiment. Namely, both potentials predict a clear monotonous increase of $\tau_G^{ph-ph}$ with $N$, although again LJ predicts lower absolute values. Also, for temperatures close to room temperature, the LJ results for 1L and graphite agree closer with TRIARS experimental data by Chatzakis et al for graphite [40], all giving a value of approximately $\tau_G^{ph-ph}(300\text{K}) \approx 2$ ps, while the KC predicted value for graphite is more than two times larger (5.6 ps). In any case, these results confirm that with increasing the number of layers, the G-mode lifetime also increases, although the amount of increase depends on the type of interlayer potential used.

<figure><img src="image_9.png" /><figcaption>FIG. 8. The effect of the coupling strength $\chi$ on the G-mode lifetimes for BG calculated with the LJ potential. The horizontal line (black) represents the calculated value for the 1LG. The vertical line (red) represents the value $\chi = 1$, which corresponds to the normal coupling strength that gives $\tau_G = 3.6$ ps. Dashed lines connecting points are b-splines as visual guides.</figcaption></figure>

To conclude this work, and in order to shed further light on how the nature of the interlayer interaction affects the lifetimes in FLG, we have performed further calculations on the G-mode lifetime of Bulk Graphite, using the LJ potential modified with a coupling strength factor $\chi$, which simply multiplies the LJ energy by $\chi$.

We have calculated the lifetime for BG with $\chi$ ranging from 0 (non-interacting single layers) to 10 (strongly interacting layers) and present it in Fig. 8. There are two main points that can be made from these results. First of all, increasing $\chi$ increases $\tau_G$, and for large enough $\chi > 2$, the lifetime reaches a plateau of ~4 ps, which could be understood as the strong-coupling limit. Secondly, the increase of the BG lifetimes with the coupling strength may indicate that the increasing trend of $\tau_G$ with $N$ is due to an increase of the interlayer interaction strength with $N$, as with increasing $N$ each layer interacts with more non-adjacent layers (mainly through the long-range part of the interlayer VdW interaction).

## IV. CONCLUSIONS:

This study presents a comprehensive systematic analysis of the phonon properties of Few-Layered Graphene. We explore the effects of the number of graphene layers, which has not been comprehensively studied computationally before, as well as the temperature dependence of phonon properties, including the first-time calculation of the temperature dependence of the Layer Breathing Mode (LBM) energy. Furthermore, we investigate the influence of the type of interlayer potential and the strength of interlayer coupling on phonon properties. Our investigation focuses particularly on phonon lifetimes, which present both experimental and computational challenges.

The behavior of inter- and intra-layer phonon modes exhibits notable differences. Inter-layer phonon energies are significantly more affected by the presence and type of interlayer interactions. As for the lifetimes, both intra-layer and LBM lifetimes increase with the number of layers, whereas for the Shear mode they are unaffected by the number of layers. One notable finding is that the Raman-active G-mode phonon lifetime $\tau_G$ increases with the number of layers $N$, regardless of the type of interlayer potential. This dependence is more pronounced with the use of the Kolmogorov-Crespi (KC) potential, while the Lennard-Jones (LJ) potential yields systematically lower values, potentially indicating increased scattering. The effect of interlayer coupling strength on G-mode lifetimes was also investigated. Increasing the LJ interlayer coupling strength leads to an increase in $\tau_G$ for graphite, up to $\chi=2$, where it reaches a plateau that indicates a strong-coupling threshold. This observation suggests that the increase in $\tau_G$ with the number of layers may be attributed to the increasing interlayer coupling strength with the addition of more layers.

The Kolmogorov-Crespi (KC) potential was proven to be highly effective in describing interlayer interactions for multi-layered graphene, particularly for the shear forces, which are challenging to model with simpler potentials. It manages to fairly reproduce the experimental Shear and Breathing mode energies, although it tends to overestimate the lifetimes of the
Raman-active G-mode in graphite by a factor of 2. Conversely, the Lennard-Jones (LJ) potential performs poorly in describing C-C interlayer interactions. It fails to simultaneously capture both shear and perpendicular interlayer forces. Particularly for the Shear mode, the LJ potential exhibits undesirable behavior, including the prediction of extremely low values of energy, which practically remained constant with temperature, and a poor Lorentzian shape of the Spectral Energy Density, leading to significant uncertainties in phonon lifetimes.

We advocate for further research combining Density Functional Theory (DFT) simulations with experimental studies to explore the effect of the number of graphene layers on phonon lifetimes due to phonon-phonon interactions. This includes DFT simulations in multilayer graphene and Time-Resolved Incoherent Anti-Stokes Raman Spectroscopy (TRIARS) experiments in high-quality suspended and monocrystalline FLG, in order to isolate the anharmonic contribution to the phonon lifetime.

## ACKNOWLEDGMENTS

The major part of the calculations presented in this work was performed in the *Jean-Zay* CPU array of the IDRIS High Performance Computing facilities (HPC) under the allocation 2022-AD010913913 made by GENCI. Part of the calculations was also performed in the *Aristotelis* HPC facilities of the Aristotle University of Thessaloniki. ENK acknowledges support for part of this research project in the framework of H.F.R.I call “Basic research Financing (Horizontal support of all Sciences)” under the National Recovery and Resilience Plan “Greece 2.0” funded by the European Union - NextGenerationEU (H.F.R.I. Project Number: 15932).

## APPENDICES:

### Appendix A: The Kolmogorov-Crespi (KC) interlayer potential:

The KC potential is given by [60]

$$V(\mathbf{r}_{ij}, \mathbf{n}_i, \mathbf{n}_j) = e^{-\lambda(r_{ij}-z_0)} [C + f(\rho_{ij}) + f(\rho_{ji})] - A \left(\frac{r_{ij}}{z_0}\right)^{-6},$$

$$\rho_{ji} = r_{ij} - (\mathbf{n}_i \mathbf{r}_{ij})^2$$

$$f(\rho) = e^{-(\rho/\delta)^2} \sum C_{2n}(\rho/\delta)^{2n}. \quad (A1)$$

where $\mathbf{n}_i$ is the local normal to the $sp^2$ plane in the vicinity of atom i, pointing to the direction of the $p_z$ orbitals of atom i and $\rho_{ji}$ the transverse distance of the $p_z$

orbitals of atoms i and j belonging to different layers. $C$, $\lambda$, $z_0$, $A$, $\delta$ and $C_{2n}$ are fitting parameters. Upon reparametrization, the KC potential has been shown to be able to reproduce very accurately the interlayer binding energy curve and sliding energy surfaces of bilayer graphene [68].

### Appendix B: The k-Space Velocity Autocorrelation (kVACS) method:

The central quantity that is calculated with this method is the Velocity Autocorrelation Sequence (k-VACS) denoted by $Z_x^p(\mathbf{k}, t_i)$. In order to arrive at its definition, we first define the k-space velocity Fourier transform $v_x^p(\mathbf{k}, t)$ given by

$$v_x^p(\mathbf{k}, t) = \frac{1}{\sqrt{N}} \sum_{l=1}^{N} v_x(\mathbf{r}_l^p, t) e^{-i\mathbf{k} \cdot \mathbf{R}_l} \quad (B1)$$

where $x$ denotes the cartesian component, $\mathbf{r}_l^p$ is a label for the atom species labeled $p$ inside the primitive unit cell labeled $l$, and $v_x(\mathbf{r}_l^p, t)$ its respective velocity at time $t$. Also $\mathbf{R}_l$ marks the position of the primitive unit cell $l$ within the computational cell. The kVACS $Z_x^p(\mathbf{k}, \omega)$ is then defined as a correlation function

$$Z_x^p(\mathbf{k}, \omega) = \int d\tau e^{i\omega\tau} \langle v_x^{*p}(\mathbf{k}, \tau) v_x^p(\mathbf{k}, 0) \rangle. \quad (B2)$$

where the brackets $\langle \cdots \rangle$ denote time averaging, the asterisk $*$ denotes complex conjugation and $\tau$ is the correlation time. It should be noted that normally the angular brackets denote ensemble averages, which are not directly computable in MD simulations, therefore we use time-averages instead. Since at the time scales of the samples we obtained (in the order of 350 ps) ergodicity is not guaranteed, we perform multiple independent simulations which we average to obtain the final $Z_x^p(\mathbf{k}, \omega)$. By adding the contributions of all atom species $p$ we obtain the polarised kVACS $Z_x(\mathbf{k}, \omega)$, and again by summing over all polarizations $x$, we get the total

$$Z_{\text{total}}(\mathbf{k}, \omega) = \sum_x \sum_p m_p Z_x^p(\mathbf{k}, \omega) \quad (B3)$$

In Appendix C we show that the kVACS is approximately a superposition of Lorentzian peaks having a central positions and FWHM's equal to the energies $\omega_{\mathbf{k},j}$ and scattering rates $2\Gamma_{\mathbf{k},j}$ respectively of the phonon mode $j$ at point $\mathbf{k}$, given by

$$Z_{\text{total}}(\mathbf{k}, \omega) = \sum_j \frac{A_j}{(\omega - \omega_{\mathbf{k},j})^2 + \Gamma_{\mathbf{k},j}^2} \quad (B4)$$

which is exactly the phonon Spectral Energy Density (SED) of the system. As it has been already mentioned before [45], taking the velocity as a Wide-Sense Stationary (WSS) random process, we can apply the Wiener-Khintchine theorem [78] to prove that
$Z_x^p(\mathbf{k}, \omega) \propto \Phi_x^p(\mathbf{k}, \omega)$, where $\Phi_x^p(\mathbf{k}, \omega)$ is the phonon SED, given by [52]

$$ \Phi_x^p(\mathbf{k}, \omega) = \left| \int dt e^{i\omega t} v_x^p(\mathbf{k}, t) \right|^2 \\ := |v_x^p(\mathbf{k}, \omega)|^2 \qquad (B5) $$

Again, by summing over all polarizations $x$ and atom species $p$, we get the total $\Phi(\mathbf{k}, \omega) = \sum_{x,p} m_p \Phi_x^p(\mathbf{k}, \omega)$. The conditions for a random process to be a WSS is to be square integrable, having a stationary mean value and an autocorrelation function that does not depend on the reference time $t_0$, which are in general fulfilled for the velocity of a stationary system in thermal equilibrium. What is more, calculating $\Phi_x(\mathbf{k}, \omega)$ instead of $Z_x(\mathbf{k}, \omega)$ is computationally more efficient, as the later via eq. (B2) is of order $O(N^2)$, while the SED in eq. (B5) is of order $O(N \log N)$, with respect to the total number $N$ of sampled time points.

Finally, $\Phi_x(\mathbf{k}, \omega)$ can be related to the phonon Density of States (pDOS) $g(\omega)$ by the relation [79]

$$ g(\omega) = \sum_x \sum_p m_p \int \Phi_x^p(\mathbf{k}, \omega) d\mathbf{k} \\ = \frac{1}{3k_B T} \sum_x \sum_p m_p \sum_l \left| \int dt e^{i\omega t} v_x(\mathbf{r}_l^p, t) \right|^2 \quad (B6) $$

## Appendix C: Proof that $Z_i(\mathbf{k}, \omega) \propto \frac{1}{(\omega - \omega_{i,0})^2 + \Gamma_i^2}$

In this appendix we shall prove that the kVACS $Z_i(\mathbf{k}, \omega)$ is approximately a superposition of Lorentzian peaks having a central positions and FWHM's equal to the energies $\omega_{\mathbf{k},j}$ and scattering rates $2\Gamma_{\mathbf{k},j}$ respectively of the phonon mode $(\mathbf{k}, j)$ where $\mathbf{k}$ is the wavevector and $j$ is the phonon branch. First, we shall recall from basic quantum mechanics that the velocity operator $v_l$ in the Heisenberg representation is given by

$$ \mathbf{u}_l(t) := e^{i\frac{H}{\hbar}t} \mathbf{u}_l e^{-i\frac{H}{\hbar}t} \qquad (C1) $$

$$ \mathbf{v}_l(t) := \frac{\mathbf{p}_l(t)}{m} = \frac{\partial \mathbf{u}_l(t)}{\partial t} = -\frac{i}{\hbar} [\mathbf{u}_l(t), H] \quad (C2) $$

where $\mathbf{p}_l(t)$ is the momentum operator, $\mathbf{u}_l(t)$ is the displacement operator and $H$ is the Hamiltonian of the system. The subscript index $l$ is an index of atoms. As a corollary, the spatial Fourier Transform (FT) of the velocity operator is given by

$$ \mathbf{v}(\mathbf{k}, t) := \sum_l \mathbf{v}_l(t) e^{-i\mathbf{k} \cdot \mathbf{R}_l} = -\frac{i}{\hbar} [\mathbf{u}_l(\mathbf{k}, t), H] \quad (C3) $$

where $\mathbf{R}_l$ is the equilibrium position of atom $l$, and is not an operator. Also $\mathbf{u}_l(\mathbf{k}, t)$ is the spatial FT of the displacement operator $\mathbf{u}_l(t)$. Now, we shall focus on the treatment of the main quantity $Z_{xy}(\mathbf{k}, \omega)$. As per

eq. (B2), it is defined as

$$ \begin{aligned} Z_{xy}(\mathbf{k}, \omega) &= \int d\tau e^{i\omega\tau} \langle v_x^\dagger(\mathbf{k}, \tau) v_y(\mathbf{k}, 0) \rangle \\ &= \int d\tau e^{i\omega\tau} \langle v_x(-\mathbf{k}, \tau) v_y^p(\mathbf{k}, 0) \rangle \quad (C4) \end{aligned} $$

where $\tau$ is the correlation time. The second line results from the relation $A^\dagger(\mathbf{k}, t) = A(-\mathbf{k}, t)$ which is true for a Hermitian operator $A_l(t)$ depending on the atom index $l$. We shall finally use the quantity $B_{xy}(\mathbf{k}, \omega)$ defined as

$$ B_{xy}(\mathbf{k}, \omega) = \int d\tau e^{i\omega\tau} \langle u_x(-\mathbf{k}, \tau) u_y(\mathbf{k}, 0) \rangle \quad (C5) $$

This quantity is closely related to the phonon propagator, which gives the phonon self-energy $\Sigma_{\mathbf{k},j}(\omega) = \Delta_{\mathbf{k},j}(\omega) + i\Gamma_{\mathbf{k},j}(\omega)$ [80]. In the seminal work by Maradudin and Fein (1962), it was proven that if the anharmonic terms in the crystal potential are small enough to be considered as perturbations, that is $\Delta_{\mathbf{k},j}, \Gamma_{\mathbf{k},j} \ll \omega_{\mathbf{k},j}$, then they can both be considered as independent of $\omega$, giving [25]

$$ \begin{aligned} B_{xy}(\mathbf{k}, \omega) = & \frac{1}{Nm} \left( \frac{\hbar}{1 - e^{-\beta\hbar\omega}} \right) \sum_j \left[ \frac{e_x^*(\mathbf{k}, j) e_y(\mathbf{k}, j)}{\omega(\mathbf{k}, j)} \right] \\ & \times \frac{\Gamma_{\mathbf{k},j}}{(\omega - \omega(\mathbf{k}, j) - \Delta_{\mathbf{k},j})^2 + \Gamma_{\mathbf{k},j}^2} \end{aligned} \quad (C6) $$

where $N$ is the number of unit cells, $m$ is the atomic mass, $e_x(\mathbf{k}, j)$ is the $x$-coordinate of the eigenvector of the phonon mode $j$ at point $\mathbf{k}$, $\beta = 1/k_B T$ and $k_B$ is the Boltzmann constant.

We shall now proceed to prove that $Z_{xy}(\mathbf{k}, \omega) = \omega^2 B_{xy}(\mathbf{k}, \omega)$. First, we recall that the thermodynamical average in the canonical ensemble of an operator $X$ is given by [80]

$$ \langle X \rangle = \frac{1}{Z} \sum_a e^{\beta E_a} \langle a | X | a \rangle \quad (C7) $$

where $Z$ is the partition function, $|a\rangle$ is the eigenstate of the Hamiltonian with energy $E_a$ and $k_B$ is the Boltzmann constant. Therefore, starting from eq. (C4) and applying the resolution of identity, we have

$$ \begin{aligned} Z_{xy}(\mathbf{k}, \omega) = & \frac{1}{Z} \sum_a e^{\beta E_a} \sum_\beta \\ & \times \int d\tau e^{i\omega\tau} \langle a | v_x(-\mathbf{k}, \tau) | \beta \rangle \langle \beta | v_y(\mathbf{k}, 0) | a \rangle \end{aligned} \quad (C8) $$

Now, focusing only on the integral of the second line in eq. (C8), if we expand the commutator of eq. (C3) and apply the Hamiltonian on the eigenstates, we get

$$ \left(-\frac{i}{\hbar}\right)^2 \int d\tau e^{i\omega\tau} \langle a | u_x(-\mathbf{k}, \tau) H - H u_x(-\mathbf{k}, \tau) | \beta \rangle \\ \times \langle \beta | u_y(\mathbf{k}, 0) H - H u_y(\mathbf{k}, 0) | a \rangle $$
$$
\begin{align*}
&= \left(-\frac{i}{\hbar}\right)^2 \int d\tau e^{i\omega\tau} \langle a|u_x(-\mathbf{k}, \tau)|\beta\rangle \langle\beta|u_y(\mathbf{k}, 0)|a\rangle \\
&\quad \times \left[-(E_a - E_\beta)^2\right] \\
&= -\left[\frac{i(E_a - E_\beta)}{\hbar}\right]^2 \int d\tau e^{i\omega\tau} \langle a|e^{i\frac{H}{\hbar}\tau}\mathbf{u}_x(-\mathbf{k})e^{-i\frac{H}{\hbar}\tau}|\beta\rangle \\
&\quad \times \langle\beta|u_y(\mathbf{k}, 0)|a\rangle \\
&= \left(\frac{E_\beta - E_a}{\hbar}\right)^2 \langle a|\mathbf{u}_x(-\mathbf{k})|\beta\rangle \langle\beta|u_y(\mathbf{k})|a\rangle \\
&\quad \times \int d\tau e^{i\left(\omega - \frac{E_\beta - E_a}{\hbar}\right)} \tag{C9} \\
&= \omega^2 \langle a|\mathbf{u}_x(-\mathbf{k})|\beta\rangle \langle\beta|u_y(\mathbf{k})|a\rangle \\
&\quad \times \delta\left(\omega - \frac{E_\beta - E_a}{\hbar}\right) \tag{C10}
\end{align*}
$$

As the integral in eq. (C9) is a Dirac delta function, the only non-zero terms in the integral of eq. (C8) that survive are those where $\omega = \frac{E_\beta - E_a}{\hbar}$, therefore we can directly replace $(\frac{E_\beta - E_a}{\hbar})^2$ by $\omega^2$, obtaining eq. (C10). From this point it is straightforward to rewrite eq. (C10) backwards as

$$ \omega^2 \int d\tau e^{i\omega\tau} \langle a|u_x(-\mathbf{k}, \tau)|\beta\rangle \langle\beta|u_y(\mathbf{k}, 0)|a\rangle \quad (C11) $$

which is exactly the integral of eq. (C8). Finally, as the $\omega^2$ factor stays outside the double summation $\sum_a \sum_\beta$, we have that

$$
\begin{align}
Z_{xy}(\mathbf{k}, \omega) &= \omega^2 \int d\tau e^{i\omega\tau} \langle u_x(-\mathbf{k}, \tau)u_y(\mathbf{k}, 0) \rangle \nonumber \\
&= \omega^2 B_{xy}(\mathbf{k}, \omega) \tag{C12}
\end{align}
$$

Therefore, we have proven that relation for the quantum-mechanical correlation functions (CF). The usual 'recipe' to obtain the classical CF from the quantum mechanical CF is to apply the limit $\hbar \to 0$ [81]. By applying this limit to the expression for $B_{xy}(\mathbf{k}, \omega)$ in eq. (C6), we obtain for the prefactor

$$ \lim_{\hbar \to 0} \frac{\hbar}{1 - e^{-\beta \hbar \omega}} = \lim_{\hbar \to 0} \frac{\hbar}{1 - (1 - \beta \hbar \omega)} = \frac{1}{\beta \omega} \quad (C13) $$

Another important point to be mentioned is the behaviour of the self energy components $\Delta_{kj}(\omega)$ and $\Gamma_{kj}(\omega)$ in the limit $\hbar \to 0$. As was pointed out in eqs (5.11 a,b) in Maradudin and Fein (1962) [25], this limit is equivalent to taking the limit $T \to \infty$ and is found to be finite and non-zero. Therefore, combining eqs (C6), (C12) and (C13), we can finally write

$$
\begin{align}
Z_{xy}(\mathbf{k}, \omega) = \frac{\omega k_B T}{Nm} & \sum_j \left[ \frac{e_x^*(\mathbf{k}, j) e_y(\mathbf{k}, j)}{\omega(\mathbf{k}, j)} \right] \nonumber \\
& \times \frac{\Gamma_{kj}}{(\omega - \omega(\mathbf{k}, j) - \Delta_{kj})^2 + \Gamma_{kj}^2} \tag{C14}
\end{align}
$$

which is not exactly Lorentzian, as it has a prefactor that depends on $\omega$. However, as per the assumption that $\Delta_{kj}, \Gamma_{kj} \ll \omega_{kj}$, the prefactor is varying much more slowly than the Lorentzians, which almost behave as delta functions, therefore we can approximate it as a constant $\omega \approx \omega_{kj}$, finally getting

$$
Z_{xy}(\mathbf{k}, \omega) = \frac{k_B T}{NM} \sum_j \frac{\Gamma_{kj} e_x^*(\mathbf{k}, j) e_y(\mathbf{k}, j)}{(\omega - \omega(\mathbf{k}, j) - \Delta_{kj})^2 + \Gamma_{kj}^2} \quad (C15)
$$

which is now a superposition of Lorentzians, with positions $(\omega(\mathbf{k}, j) + \Delta_{kj})$ and FWHM $2\Gamma_{kj}$ as was to be proven.

As a sidenote, the above derivation was done for a crystal having a single atom species in its unit cell. It is straightforward to extend this to multiple species, by simply replacing $m \to m_p$ and $e_x(\mathbf{k}, j) \to e_x^p(\mathbf{k}, j)$, where now $p$ is an index to atom species. Then, by remembering that for each phonon mode $(\mathbf{k}, j)$ the eigenvectors are orthonormal $\sum_{x,p} e_x^p(\mathbf{k}, j) e_x^p(\mathbf{k}, j) = 1$, we can write $Z_{total}(\mathbf{k}, j)$ as

$$
\begin{align}
Z_{\text{total}}(\mathbf{k}, \omega) &= \sum_{x,p} m_p Z_{xx}(\mathbf{k}, \omega) \nonumber \\
&= \frac{k_B T}{N} \sum_j \frac{\Gamma_{kj}}{(\omega - \omega(\mathbf{k}, j) - \Delta_{kj})^2 + \Gamma_{kj}^2} \tag{C16}
\end{align}
$$

where we have used the correspondence between eqs (B3) and (C4) as $Z_x(\mathbf{k}, \omega) = Z_{xx}(\mathbf{k}, \omega)$

[1] A. R. Urade, I. Lahiri, and K. S. Suresh, Graphene properties, synthesis and applications: A review, JOM 75, 614 (2023).

[2] K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, Y. Zhang, S. V. Dubonos, I. V. Grigorieva, and A. A. Firsov, Electric field effect in atomically thin carbon films, Science 306, 666 (2004).

[3] J. Wang, X. Mu, M. Sun, and T. Mu, Optoelectronic properties and applications of graphene-based hybrid nanomaterials and van der Waals heterostructures, Applied Materials Today 16, 1 (2019).

[4] C. Lee, X. Wei, J. W. Kysar, and J. Hone, Measurement of the elastic properties and intrinsic strength of monolayer graphene, Science 321, 385 (2008).
[5] M. Sang, J. Shin, K. Kim, and K. J. Yu, Electronic and thermal properties of graphene and recent advances in graphene based electronics applications, Nanomaterials **9**, 10.3390/nano09030374 (2019).

[6] K. S. Novoselov, D. Jiang, F. Schedin, T. J. Booth, V. V. Khotkevich, S. V. Morozov, and A. K. Geim, Two-dimensional atomic crystals, Proceedings of the National Academy of Sciences **102**, 10451 (2005).

[7] H. Murata, Y. Nakajima, N. Saitoh, N. Yoshizawa, T. Suemasu, and K. Toko, High-electrical-conductivity multilayer graphene formed by layer exchange with controlled thickness and interlayer, Scientific Reports **9**, 4068 (2019).

[8] D. Singh, J. Y. Murthy, and T. S. Fisher, Mechanism of thermal conductivity reduction in few-layer graphene, Journal of Applied Physics **110**, 044317 (2011).

[9] J. Liu, T. Wang, S. Xu, P. Yuan, X. Xu, and X. Wang, Thermal conductivity of giant mono- to few-layered cvd graphene supported on an organic substrate, Nanoscale **8**, 10298 (2016).

[10] S. Hamze, N. Berrada, D. Cabaleiro, A. Desforges, J. Ghanbaja, J. Gleize, D. Bégin, F. Michaux, T. Maré, B. Vigolo, and P. Estellé, Few-layer graphene-based nanofluids with enhanced thermal conductivity, Nanomaterials **10**, 10.3390/nano10071258 (2020).

[11] H. Fu, G. Yan, M. Li, H. Wang, Y. Chen, C. Yan, C.-T. Lin, N. Jiang, and J. Yu, Graphene as a nanofiller for enhancing the tribological properties and thermal conductivity of base grease, RSC Adv. **9**, 42481 (2019).

[12] C. Androulidakis, E. N. Koukaras, G. Paterakis, G. Trakakis, and C. Galiotis, Tunable macroscale structural superlubricity in two-layer graphene via strain engineering, Nature Communications **11**, 1595 (2020).

[13] K. M. Shahil and A. A. Balandin, Thermal properties of graphene and multilayer graphene: Applications in thermal interface materials, Solid State Communications **152**, 1331 (2012), exploring Graphene, Recent Research Advances.

[14] Y. Cao, V. Fatemi, S. Fang, K. Watanabe, T. Taniguchi, E. Kaxiras, and P. Jarillo-Herrero, Unconventional superconductivity in magic-angle graphene superlattices, Nature **556**, 43 (2018).

[15] S. M. Davis, Y.-Z. Chou, F. Wu, and S. Das Sarma, Phonon-limited resistivity of multilayer graphene systems, Phys. Rev. B **107**, 045426 (2023).

[16] D. L. Nika, E. P. Pokatilov, and A. A. Balandin, Theoretical description of thermal transport in graphene: The issues of phonon cut-off frequencies and polarization branches, physica status solidi (b) **248**, 2609 (2011).

[17] S. Ghosh, I. Calizo, D. Teweldebrhan, E. P. Pokatilov, D. L. Nika, A. A. Balandin, W. Bao, F. Miao, and C. N. Lau, Extremely high thermal conductivity of graphene: Prospects for thermal management applications in nanoelectronic circuits, Applied Physics Letters **92**, 151911 (2008).

[18] L. Lindsay, D. A. Broido, and N. Mingo, Flexural phonons and thermal transport in multilayer graphene and graphite, Phys. Rev. B **83**, 235428 (2011).

[19] R. Nemanich, G. Lucovsky, and S. Solin, Infrared active optical vibrations of graphite, Solid State Communications **23**, 117 (1977).

[20] A. S. Davydov, *Theory of Molecular Excitons* (Springer New York, NY, 1971).

[21] P. N. Ghosh, Davydov splitting and multipole interactions, Solid State Communications **19**, 639 (1976).

[22] A. Molina-Sánchez, K. Hummer, and L. Wirtz, Vibrational and optical properties of mos2: From monolayer to bulk, Surface Science Reports **70**, 554 (2015).

[23] C. Gibaja, D. Rodriguez-San-Miguel, P. Ares, J. Gómez-Herrero, M. Varela, R. Gillen, J. Maultzsch, J. Angewandte Chemie International Edition **55**, 14345 (2016).

[24] S. K. Saha, U. V. Waghmare, H. R. Krishnamurthy, and A. K. Sood, Phonons in few-layer graphene and interplanar interaction: A first-principles study, Phys. Rev. B **78**, 165421 (2008).

[25] A. A. Maradudin and A. E. Fein, Scattering of neutrons by an anharmonic crystal, Phys. Rev. **128**, 2589 (1962).

[26] O. Hellman, P. Steneteg, I. A. Abrikosov, and S. I. Simak, Temperature dependent effective potential method for accurate free energy calculations of solids, Phys. Rev. B **87**, 104111 (2013).

[27] X. Gu, Z. Fan, H. Bao, and C. Y. Zhao, Revisiting phonon-phonon scattering in single-layer graphene, Phys. Rev. B **100**, 064306 (2019).

[28] J. Maultzsch, S. Reich, C. Thomsen, H. Requardt, and P. Ordejón, Phonon dispersion in graphite, Phys. Rev. Lett. **92**, 075501 (2004).

[29] S. Sullivan, A. Vallabhaneni, I. Kholmanov, X. Ruan, J. Murthy, and L. Shi, Optical generation and detection of local nonequilibrium phonons in suspended graphene, Nano Letters **17**, 2049 (2017).

[30] I. Calizo, A. A. Balandin, W. Bao, F. Miao, and C. N. Lau, Temperature dependence of the raman spectra of graphene and graphene multilayers, Nano Letters **7**, 2645 (2007).

[31] Z. Han, X. Yang, S. E. Sullivan, T. Feng, L. Shi, W. Li, and X. Ruan, Raman linewidth contributions from four-phonon and electron-phonon interactions in graphene, Phys. Rev. Lett. **128**, 045901 (2022).

[32] A. C. Ferrari, J. C. Meyer, V. Scardaci, C. Casiraghi, M. Lazzeri, F. Mauri, S. Piscanec, D. Jiang, K. S. Novoselov, S. Roth, and A. K. Geim, Raman spectrum of graphene and graphene layers, Phys. Rev. Lett. **97**, 187401 (2006).

[33] D. J. Late, U. Maitra, L. S. Panchakarla, U. V. Waghmare, and C. N. R. Rao, Temperature effects on the raman spectra of graphenes: dependence on the number of layers and doping, Journal of Physics: Condensed Matter **23**, 055303 (2011).

[34] P. H. Tan, W. P. Han, W. J. Zhao, Z. H. Wu, K. Chang, H. Wang, Y. F. Wang, N. Bonini, N. Marzari, N. Pugno, G. Savini, A. Lombardo, and A. C. Ferrari, The shear mode of multilayer graphene, Nature Materials **11**, 10.1038/nmat3245 (2012).

[35] C. H. Lui, L. M. Malard, S. Kim, G. Lantz, F. E. Laverge, R. Saito, and T. F. Heinz, Observation of layer-breathing mode vibrations in few-layer graphene through combination raman scattering, Nano Letters **12**, 5539 (2012), pMID: 22963681.

[36] C. H. Lui, Z. Ye, C. Keiser, X. Xiao, and R. He, Temperature-activated layer-breathing vibrations in few-layer graphene, Nano Letters **14**, 4615 (2014), pMID: 25019216.

[37] J. Lin, L. Guo, Q. Huang, Y. Jia, K. Li, X. Lai,
and X. Chen, Anharmonic phonon effects in ra-
man spectra of unsupported vertical graphene sheets,
Phys. Rev. B **83**, 125430 (2011).

[38] S. Katsiaounis, A. V. Sharkov, E. V. Khoroshilov, G. Paterakis, J. Parthenios, and K. Papagelis, Time-resolved raman scattering in exfoliated and cvd graphene crystals, The Journal of Physical Chemistry C **125**, 21003 (2021).

[39] K. Kang, D. Abdula, D. G. Cahill, and M. Shim, Lifetimes of optical phonons in graphene and graphite by time-resolved incoherent anti-stokes raman scattering, Phys. Rev. B **81**, 165405 (2010).

[40] I. Chatzakis, H. Yan, D. Song, S. Berciaud, and T. F. Heinz, Temperature dependence of the anharmonic decay of optical phonons in carbon nanotubes and graphite, Phys. Rev. B **83**, 205411 (2011).

[41] H. Wang, J. H. Strait, P. A. George, S. Shivaraman, V. B. Shields, M. Chandrashekhar, J. Hwang, F. Rana, M. G. Spencer, C. S. Ruiz-Vargas, and J. Park, Ultrafast relaxation dynamics of hot optical phonons in graphene, Applied Physics Letters **96**, 081917 (2010).

[42] H. Yan, D. Song, K. F. Mak, I. Chatzakis, J. Maultzsch, and T. F. Heinz, Time-resolved raman spectroscopy of optical phonons in graphite: Phonon anharmonic coupling and anomalous stiffening, Phys. Rev. B **80**, 121403 (2009).

[43] S. Kumar, M. Anija, N. Kamaraju, K. S. Vasu, K. S. Subrahmanyam, A. K. Sood, and C. N. R. Rao, Femtosecond carrier dynamics and saturable absorption in graphene suspensions, Applied Physics Letters **95**, 191911 (2009).

[44] N. Bonini, M. Lazzeri, N. Marzari, and F. Mauri, Phonon anharmonicities in graphite and graphene, Phys. Rev. Lett. **99**, 176802 (2007).

[45] E. N. Koukaras, G. Kalosakas, C. Galiotis, and K. Papagelis, Phonon properties of graphene derived from molecular dynamics simulations, Scientific Reports **5**, 12923 (2015).

[46] D. W. Noid, M. L. Koszykowski, and R. A. Marcus, A spectral analysis method of obtaining molecular spectra from classical trajectories, The Journal of Chemical Physics **67**, 404 (2008).

[47] A. Rahman, Correlations in the motion of atoms in liquid argon, Phys. Rev. **136**, A405 (1964).

[48] J. Kohanoff, Phonon spectra from short non-thermally equilibrated molecular dynamics simulations, Computational Materials Science **2**, 221 (1994).

[49] C. Z. Wang, C. T. Chan, and K. M. Ho, Molecular-dynamics study of anharmonic effects in silicon, Phys. Rev. B **40**, 3390 (1989).

[50] N. Papanicolaou, I. Lagaris, and G. Evangelakis, Modification of phonon spectral densities of the (001) copper surface due to copper adatoms by molecular dynamics simulation, Surface Science **337**, L819 (1995).

[51] P. D. Ditlevsen, P. Stoltze, and J. K. No/rskov, Anharmonicity and disorder on the cu(110) surface, Phys. Rev. B **44**, 13002 (1991).

[52] J. A. Thomas, J. E. Turney, R. M. Iutzi, C. H. Amon, and A. J. H. McGaughey, Predicting phonon dispersion relations and lifetimes from the spectral energy density, Phys. Rev. B **81**, 081411 (2010).

[53] M. Mafakheri and A. A. S. Dodaran, Graphene phonons lifetime and mean free path using ab initio molecular dynamics and spectral energy density analysis, Physica Scripta **96**, 065702 (2021).

[54] J. Chen and X. Zhang, Nanoscale size effect and phonon properties of silicon material through simple spectral energy density analysis based on molecular dynamics, Journal of Physics: Condensed Matter **31**, 425701 (2019).

[55] P. Desmarchelier, A. Tanguy, and K. Termentzidis, Thermal rectification in asymmetric two-phase nanowires, Phys. Rev. B **103**, 014202 (2021).

[56] Y. Zhang, Q. Lv, H. Wang, S. Zhao, Q. Xiong, R. Lv, and X. Zhang, Simultaneous electrical and thermal rectification in a monolayer lateral heterojunction, Science **378**, 169 (2022).

[57] X.-K. Chen, M. Pang, T. Chen, D. Du, and K.-Q. Chen, Thermal rectification in asymmetric graphene/hexagonal boron nitride van der waals heterostructures, ACS Applied Materials & Interfaces **12**, 15517 (2020), pMID: 32153173.

[58] M. Born and K. Huang, *Dynamical theory of crystal lattices*, International series of monographs on physics (Clarendon Press Oxford, Oxford, 1956).

[59] J. Adams, Bonding energy models, in *Encyclopedia of Materials: Science and Technology*, edited by K. J. Buschow, R. W. Cahn, M. C. Flemings, B. Ilschner, E. J. Kramer, S. Mahajan, and P. Veyssière (Elsevier, Oxford, 2001) pp. 763–767.

[60] A. N. Kolmogorov and V. H. Crespi, Registry-dependent interlayer potential for graphitic systems, Phys. Rev. B **71**, 235415 (2005).

[61] I. Leven, I. Azuri, L. Kronik, and O. Hod, Inter-layer potential for hexagonal boron nitride, The Journal of Chemical Physics **140**, 104106 (2014).

[62] M. Wen, S. Carr, S. Fang, E. Kaxiras, and E. B. Tadmor, Dihedral-angle-corrected registry-dependent interlayer potential for multilayer graphene structures, Phys. Rev. B **98**, 235404 (2018).

[63] I. Lebedeva, A. Knizhnik, A. Popov, Y. Lozovik, and B. Potapkin, Modeling of graphene-based nems, Physica E: Low-dimensional Systems and Nanostructures **44**, 949 (2010) (the proceedings of the European Materials Research Symposium on Science and Technology of Nanotubes, Nanowires and Graphene).

[64] M. Mohr, J. Maultzsch, E. Dobardžić, S. Reich, I. Milošević, M. Damnjanović, A. Bosak, M. Krisch, and C. Thomsen, Phonon dispersion of graphite by inelastic x-ray scattering, Phys. Rev. B **76**, 035439 (2007).

[65] I. I. Al-Qasir, A. A. Campbell, G. Sala, J. Y. Lin, Y. Cheng, F. F. Islam, D. L. Abernathy, and M. B. Stone, Vacancy-driven variations in the phonon density of states of fast neutron irradiated nuclear graphite, Carbon **168**, 42 (2020).

[66] L. Lindsay and D. A. Broido, Optimized tersoff and brenner empirical potential parameters for lattice dynamics and phonon thermal transport in carbon nanotubes and graphene, Phys. Rev. B **81**, 205441 (2010).

[67] P. Rowe, G. Csányi, D. Alfe, and A. Michaelides, Development of a machine learning potential for graphene, Phys. Rev. B **97**, 054303 (2018).

[68] W. Ouyang, D. Mandelli, M. Urbakh, and O. Hod, Nanoserpers: Graphene nanoribbon motion on two-dimensional hexagonal materials, Nano Letters **18**, 6009 (2018).

[69] T. Gould, Z. Liu, J. Z. Liu, J. F. Dobson, Q. Zheng, and S. Lebègue, Binding and in-
terlayer force in the near-contact region of
two graphite slabs: Experiment and theory,
The Journal of Chemical Physics **139**, 224704 (2013).

[70] W. Wang, S. Dai, X. Li, J. Yang, D. J. Srolovitz, and Q. Zheng, Measurement of the cleavage energy of graphite, Nature Communications **6**, 7853 (2015).

[71] See Supplemental Material at URL-will-be-inserted-by-publisher for the data of the experiments.

[72] A. P. Thompson, H. M. Aktulga, R. Berger, D. S. Bolintineanu, W. M. Brown, P. S. Crozier, P. J. in 't Veld, A. Kohlmeyer, S. G. Moore, T. D. Nguyen, R. Shan, M. J. Stevens, J. Tranchida, C. Trott, and S. J. Plimpton, Lammps - a flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales, Computer Physics Communications **271**, 108171 (2022).

[73] G. Pizzi, S. Milana, A. C. Ferrari, N. Marzari, and M. Gibertini, Shear and breathing modes of layered materials, ACS Nano **15**, 12509 (2021).

[74] C. Cong and T. Yu, Enhanced ultra-low-frequency interlayer shear modes in folded graphene layers, Nature Communications **5**, 4709 (2014).

[75] M. Yang, L. Wang, X. Qiao, Y. Liu, Y. Liu, Y. Shi, H. Wu, B. Liang, X. Li, and X. Zhao,

Temperature dependence of g and d' phonons in
monolayer to few-layer graphene with vacancies,
Nanoscale Research Letters **15**, 189 (2020).

[76] H.-N. Liu, X. Cong, M.-L. Lin, and P.-H. Tan, The intrinsic temperature-dependent raman spectra of graphite in the temperature range from 4k to 1000k, Carbon **152**, 451 (2019).

[77] M. Balkanski, R. F. Wallis, and E. Haro, Anharmonic effects in light scattering due to optical phonons in silicon, Phys. Rev. B **28**, 1928 (1983).

[78] R. Blackman and J. Tukey, *The Measurement of Power Spectra, from the Point of View of Communications Engineering*, Dover books on electronics, electricity, computers, electrical engineering (Dover Publications, 1959).

[79] H. Chakraborty, S. Mogurampelly, V. K. Yadav, U. V. Waghmare, and M. L. Klein, Phonons and thermal conducting properties of borocarbonitride (bcn) nanosheets, Nanoscale **10**, 22148 (2018).

[80] G. Mahan, *Many-Particle Physics*, Physics of Solids and Liquids (Springer US, 1990).

[81] R. Zwanzig, Time-correlation functions and transport coefficients in statistical mechanics, Ann. Rev. Phys. Chem. **16**, 67 (1965).