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

Markos Poulos<EMAIL_ADDRESS>INSA Lyon, CNRS, CETHIL, UMR5008,
69621 Villeurbanne, France  Konstantinos Papagelis  Aristotle University of
Thessaloniki, School of Physics, 54124 Thessaloniki, Greece  Institute of
Chemical Engineering Sciences, FORTH/ICE-HT, 26504 Patras, Greece  Emmanuel N.
Koukaras  Institute of Chemical Engineering Sciences, FORTH/ICE-HT, 26504
Patras, Greece  Aristotle University of Thessaloniki, School of Chemistry,
54124 Thessaloniki, Greece  George Kalosakas  University of Patras, Department
of Materials Science, GR-26504 Rio, Greece  Giorgia Fugallo  CNRS, Université
de Nantes, LTEN, UMR6607, 44306 Nantes, France  Konstantinos Termentzidis
CNRS, INSA Lyon, CETHIL, UMR5008, 69621 Villeurbanne, France

###### Abstract

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 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.

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

††preprint: APS/123-QED

##  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, 9, 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 E2g 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
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 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 readers to the work of Saha et al. [24].

The dependence on temperature $T$ of the phonon energies $\omega(T)$ and
lifetimes $\tau(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 $\Sigma(\omega)$ of an interacting phonon
system [25], one can obtain from its real $\Delta(T)$ and imaginary part
$\Gamma(T)$, respectively, the anharmonic phonon frequency shift and finite
linewidth in the energies $\omega_{0}$ of the harmonic non-interacting
phonons. The T-dependence then mainly arises due to the dependence of $\Gamma$
and $\Delta$ on the phonon populations $n(\omega,T)$. Lattice thermal
expansion also contributes to $\Sigma(\omega)$ by the implicit inclusion of
higher order anharmonic terms in the renormalization of the harmonic energy
[26] leading to a temperature-dependent $\omega_{0}(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 ($\omega_{G}\approx 1580$ cm-1). For graphite, the
phonon dispersion curves by means of Inelastic X-ray Scattering (IXS) have
been reported [28], and for 1L graphene studies on the T-dependence of
$\omega_{G}(T)$ by means of conventional Raman spectroscopy are also present
in the literature, with excellent agreement between them [29, 30, 31]. Studies
of the effect of the number of layers $N$ on the phonon energies of FLG have
also been conducted [32, 33, 34, 35, 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 [37, 31]. The lifetime of a
Raman-active mode can be estimated by the linewidth $\gamma$ in cm-1 as
$\tau=\left(2\gamma\pi c\right)^{-1}$ [38, 39]. The major contributions to the
linewidth are due to el-ph and ph-ph interactions $\gamma_{G}=\gamma_{G}^{ph-
ph}+\gamma_{G}^{el-ph}$. Most notably for the G-mode, the major contribution
to $\gamma_{G}$ is due to the electron-phonon part [40, 41, 31]. 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 $\tau^{ph-ph}_{G}$ ranges between 2-2.5 ps
at room temperature [40, 42]. For 1L, however, reported values are dispersed,
from $\sim$1-1.5 ps [39, 38] 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 calculations in the
literature included only the electron-phonon (e$\mathchar 45\relax$ph)
andthree-phonon (3ph) contributions to the temperature-dependence of
$\tau_{G}(T)$ in single-layer graphene [44]. These calculations showed that
the e$\mathchar 45\relax$ph contribution overwhelmingly dominates the
temperature dependence,causing $\gamma_{G}$ to decrease with temperature. This
finding starkly contrasts with all the expermental 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$\mathchar 45\relax$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$\mathchar 45\relax$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 $\Gamma$-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 [46, 45]. 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 refered 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 MoSe2-WSe2 [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 interaction 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 $\pi$-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 potenials 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 $\Gamma$-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.

Figure 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.

##  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

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

For the KC potential, a global cutoff of rc = 12 $\AA$ was chosen. For the LJ
potential, the parameters $\epsilon$= 4.6 meV and $\sigma$ = 3.276 $\AA$ 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 intdividual 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 $\sim$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 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}}\textbf{b}_{1}$ and $\frac{1}{N_{2}}\textbf{b}_{2}$, where
$\textbf{b}_{1}$ and $\textbf{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 vibrational 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 accross 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-1, 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-1). 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 $\sim$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 $\sim$5-20 ps and increasing
close to the $\Gamma$-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 overstimated by $\sim$10%, while the ZO
branch is significantly overestimated at about 45%. Also, the description of
the LO and TO branches is satisfactory on the $\Gamma$-M path (armchair), but
significantly overestimates about 27% the splitting of the LO branch at the
$\Gamma$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 ($\Gamma$-A) of the Brillouin Zone of graphite
was taken, those modes are therefore missing from the PDOS.

Figure 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).

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 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
experimental 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]

$\displaystyle\omega^{(i,N)}_{\nu}$ $\displaystyle=$
$\displaystyle\sqrt{\frac{2k_{\nu}}{\mu}\left\\{1-\cos\left(\frac{\left(i-1\right)\pi}{N}\right)\right\\}}$
(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\dotsc 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$\cdot\AA^{-2}$ . The model
has been succesfull 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

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

$\displaystyle\omega_{LBM}(N)$ $\displaystyle=$ $\displaystyle\frac{1}{\pi
c}\sqrt{\frac{b}{\mu}}\sin\left(\frac{\pi}{2N}\right)$ (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$\cdot$s-1 is the speed of light. It should also be
noted that in eqs. (3) and (4) the pre-factors of the sinusoidals 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 Supplemental
Material [71]), which for $N\to\infty$ tends to zero.

Figure 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$.

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$\cdot$m-2 for the interlayer interaction,
which is very close to the value of $\alpha_{exp}=1.27\times 10^{19}$
N$\cdot$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$\cdot$m-2, while
$b_{exp}=8.84\times 10^{19}$ N$\cdot$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$\cdot$m-2 and
$b_{LD}=10.18\times 10^{19}$ N$\cdot$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 of 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 1.

Table 1:  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-1K-1 # of Layers | Shear | LBM | ZO | G
---|---|---|---|---
($\times$10-3 cm-1K-1)
1 | – | – | -36 | -51
2 | -6.2 | -9.7 | -35 | -47
3 | -5.9 | -6.3 | -34 | -48
5 | -5.7 | -3.6 | -35 | -49
BG | -4.4 | -5.6 | -35 | -50
Exper. | – | – | – | -21

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}$ cm-1 T-1 for bilayer, down to $-4.4\cdot
10^{-3}$ cm-1 T-1 for graphite. The same trend is observed for the LBM, where
the slopes are ranging from $-9.7\cdot 10^{-3}$ cm-1 T-1 for bilayer, down to
$-3.6\cdot 10^{-3}$ cm-1 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}$ cm-1 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.

Figure 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

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 $\sim$1-5 cm -1 going from 1L to graphite. Experimental
measurements on the G mode also confirm these findings [32, 33, 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 then for the
interlayer modes (see Table 1). 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 1).

By comparing with the experimental measurements, we notice again that the
Tersoff-2010 potential significantly overestimates the G-mode energies by a
factor of $\sim$5.5 %, and the slope $\chi_{T}$ by a factor of $\sim$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 the kVACS peaks (FWHM), by
the simple relation

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

For convenience with units, we use $\hbar=5.29$ ps$\cdot$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 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^{ph-ph}_{C}$ [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 ps while the calculated $\tau_{G}^{ph-
ph}$=5.6 ps, more than double that of the monolayer. Although to our knowledge
no experimental study has so far directly measured the $\tau^{ph-ph}_{G}(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 ps. For graphite, $\tau^{ph-ph}_{G}(T)$ has been
experimentally measured at around $\sim$2.5 ps by multiple authors [40, 42],
which is almost 2 times lower than our calculated values with the KC potential
($\sim$5.6 ps). For supported 1L samples, values of $\sim$1$\mathchar
45\relax$1.5 ps have been reported, increasing with increasing the number of
layers $N$ [39, 38], while some authors have measured $\tau^{ph-
ph}_{G}\mathchar 45\relax$2.5 ps, equal to that of graphite [41]. Values for
1L suspended on colloidal solutions as high as $4.9\ ps$ for $\tau^{ph-
ph}_{G}$ have also been reported [43].

Figure 5:  The calculated temperature dependence of $\tau^{ph-ph}$ of the
monolayer Raman-active G phonon (black triangles). DFT calculations including
3ph+4ph processes and phonon renormalizaton [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 Figure 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.

We have also compared our calculations for $\tau^{ph-ph}_{G}$ 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 $\sim 40\%$ lower. The experimental values exhibit a
different trend, with $\tau^{ph-ph}_{G}$ remaining almost constant for a
temperature $300K\leq T\leq 600K$. The authors have also reported that
$\tau^{ph-ph}_{G}$ 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\leq
T\leq 900K$. The slopes $\chi_{T}$ were also calculated and are presented in
Table 2. 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.

Table 2:  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) ($\times$10-3 cm-1K-1) | KC | LJ
---|---|---
Bilayer | Graphite | Bilayer | Graphite
Shear | -6.2 | -4.3 | -0.4 | -0.5
LBM | -9.7 | -5.6 | -11.8 | -6.8
ZO | -35.4 | -34.8 | -35.8 | -35.2
G | -47.9 | -50.0 | -52.7 | -47.9
Figure 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.

In Figs. 6a) and 6b) we present the comparison for the interlayer mode
energies, where the trends of the phonon 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 ($\sim 10^{-3}$ cm-1K-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 2.

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 2. 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 ($\sim 3\cdot 10^{-3}$) than with KC ($\sim 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^{ph-ph}_{C}(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 compared to KC. As for the LBM lifetimes (Fig. 7b), again LJ gives
systematically underestimates $\tau^{ph-ph}_{LBM}(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^{ph-ph}_{ZO}(T)$ for 1L and 2L are
practically equal, irrespectful 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^{BG}_{ZO}>\tau^{1L}_{ZO},\tau^{2L}_{ZO}$

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^{ph-ph}_{G}$ 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^{ph-
ph}_{G}(300K)\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 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.

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
$\sim$4 ps, which could be understood as the strong-coupling limit. Secondly,
the increae 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
ressearch 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]

$\displaystyle V(\textbf{r}_{ij},\textbf{n}_{i},\textbf{n}_{j})$
$\displaystyle=$ $\displaystyle
e^{-\lambda(r_{ij}-z_{0})}\left[C+f(\rho_{ij})+f(\rho_{ji})\right]$
$\displaystyle-$ $\displaystyle A\left(\frac{r_{ij}}{z_{0}}\right)^{-6},$
$\displaystyle\rho_{ji}$ $\displaystyle=$ $\displaystyle
r_{ij}-\left(\textbf{n}_{i}\textbf{r}_{ij}\right)^{2}$ $\displaystyle f(\rho)$
$\displaystyle=$ $\displaystyle e^{-(\rho/\delta)^{2}}\sum
C_{2n}(\rho/\delta)^{2n}.$ (6)

where $\textbf{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}(\textbf{k},t_{i})$. In
order to arrive at its definition, we first define the k-space velocity
Fourier transform $v_{x}^{p}(\textbf{k},t)$ given by

$\displaystyle v_{x}^{p}(\textbf{k},t)$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{N}}\sum_{l=1}^{N}v_{x}(\textbf{r}_{l}^{p},t)e^{-i\textbf{k}\cdot\textbf{R}_{l}}$
(7)

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

$\displaystyle Z_{x}^{p}(\textbf{k},\omega)$ $\displaystyle=$
$\displaystyle\int d\tau e^{i\omega\tau}\langle
v_{x}^{*p}(\textbf{k},\tau)v_{x}^{p}(\textbf{k},0)\rangle.$ (8)

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}(\textbf{k},\omega)$. By adding the contributions of all atom
species $p$ we obtain the polarised kVACS $Z_{x}(\textbf{k},\omega)$, and
again by summing over all polarizations $x$, we get the total

$\displaystyle Z_{total}(\textbf{k},\omega)$ $\displaystyle=$
$\displaystyle\sum_{x}\sum_{p}m_{p}Z_{x}^{p}(\textbf{k},\omega)$ (9)

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_{\textbf{k},j}$ and scattering rates $2\Gamma_{\textbf{k},j}$
respectively of the phonon mode $j$ at point k, given by

$\displaystyle Z_{total}(\textbf{k},\omega)$ $\displaystyle=$
$\displaystyle\sum_{j}\frac{A_{j}}{\left(\omega-\omega_{\textbf{k},j}\right)^{2}+\Gamma_{\textbf{k},j}^{2}}$
(10)

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}(\textbf{k},\omega)\propto\Phi_{x}^{p}(\textbf{k},\omega)$, where
$\Phi_{x}^{p}(\textbf{k},\omega)$ is the phonon SED, given by [52]

$\displaystyle\Phi_{x}^{p}(\textbf{k},\omega)$ $\displaystyle=$
$\displaystyle\left|\int dte^{i\omega t}v_{x}^{p}(\textbf{k},t)\right|^{2}$
(11) $\displaystyle\coloneqq$
$\displaystyle\left|v_{x}^{p}(\textbf{k},\omega)\right|^{2}$

Again, by summing over all polarizations $x$ and atom species $p$, we get the
total
$\Phi(\textbf{k},\omega)=\sum_{x,p}m_{p}\Phi_{x}^{p}(\textbf{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}(\textbf{k},\omega)$ instead of $Z_{x}(\textbf{k},\omega)$ is
computationally more efficient, as the later via eq. (8) is of order
$\mathcal{O}(N^{2})$, while the the SED in eq. (11) is of order
$\mathcal{O}(NlogN)$, with respect to the total number $N$ of sampled time
points.

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

$\displaystyle
g(\omega)=\sum_{x}\sum_{p}m_{p}\int\Phi_{x}^{p}(\textbf{k},\omega)d\textbf{k}$
$\displaystyle=\frac{1}{3k_{B}T}\sum_{x}\sum_{p}m_{p}\sum_{l}\left|\int
dte^{i\omega t}v_{x}(\textbf{r}_{l}^{p},t)\right|^{2}$ (12)

##  Appendix C Proof that
$Z_{i}(\textbf{k},\omega)\propto\frac{1}{\left(\omega-\omega_{i,0}\right)^{2}+\Gamma_{i}^{2}}$

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

$\displaystyle\textbf{u}_{l}(t)$ $\displaystyle\coloneq$ $\displaystyle
e^{i\frac{H}{\hbar}t}\textbf{u}_{l}e^{-i\frac{H}{\hbar}t}$ (13)
$\displaystyle\textbf{v}_{l}(t)$ $\displaystyle\coloneq$
$\displaystyle\frac{\textbf{p}_{l}(t)}{m}=\frac{\partial\textbf{u}_{l}(t)}{\partial
t}=-\frac{i}{\hbar}\left[\textbf{u}_{l}(t),H\right]$ (14)

where $\textbf{p}_{l}(t)$ is the momentum operator, $\textbf{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

$\textbf{v}(\textbf{k},t)\coloneq\sum_{l}\textbf{v}_{l}(t)e^{-i\textbf{k}\cdot\textbf{R}_{l}}=-\frac{i}{\hbar}\left[\textbf{u}_{l}(\textbf{k},t),H\right]$
(15)

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

$\displaystyle Z_{xy}(\textbf{k},\omega)$ $\displaystyle=$ $\displaystyle\int
d\tau e^{i\omega\tau}\langle
v_{x}^{\dagger}(\textbf{k},\tau)v_{y}(\textbf{k},0)\rangle$ (16)
$\displaystyle=$ $\displaystyle\int d\tau e^{i\omega\tau}\langle
v_{x}(-\textbf{k},\tau)v_{y}^{p}(\textbf{k},0)\rangle$

where $\tau$ is the correlation time. The second line results from the
relation $A^{\dagger}(\textbf{k},t)=A(-\textbf{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}(\textbf{k},\omega)$ defined as

$B_{xy}(\textbf{k},\omega)=\int d\tau e^{i\omega\tau}\langle
u_{x}(-\textbf{k},\tau)u_{y}(\textbf{k},0)\rangle$ (17)

This quantity is closely related to the phonon propagator, which gives the
phonon self-energy
$\Sigma_{\textbf{k}j}(\omega)=\Delta_{\textbf{k}j}(\omega)+i\Gamma_{\textbf{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_{\textbf{k}j},\Gamma_{\textbf{k}j}\ll\omega_{\textbf{k}j}$, then they
can both be considered as independent of $\omega$, giving [25]

$\displaystyle B_{xy}(\textbf{k},\omega)$ $\displaystyle=$
$\displaystyle\frac{1}{Nm}\left(\frac{\hbar}{1-e^{-\beta\hbar\omega}}\right)\sum_{j}\left[\frac{e^{*}_{x}(\textbf{k},j)e_{y}(\textbf{k},j)}{\omega(\textbf{k},j)}\right]$
(18) $\displaystyle\times$
$\displaystyle\frac{\Gamma_{\textbf{k}j}}{\left(\omega-\omega(\textbf{k},j)-\Delta_{\textbf{k}j}\right)^{2}+\Gamma^{2}_{\textbf{k}j}}$

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

We shall now proceed to prove that
$Z_{xy}(\textbf{k},\omega)=\omega^{2}B_{xy}(\textbf{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$ (19)

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. (16) and applying the resolution of identity, we
have

$\displaystyle Z_{xy}$
$\displaystyle(\textbf{k},\omega)=\frac{1}{Z}\sum_{a}e^{\beta
E_{a}}\sum_{\beta}$ (20) $\displaystyle\times$ $\displaystyle\int d\tau
e^{i\omega\tau}\langle
a|v_{x}(-\textbf{k},\tau)|\beta\rangle\langle\beta|v_{y}(\textbf{k},0)|a\rangle$

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

$\displaystyle\left(-\frac{i}{\hbar}\right)^{2}\int d\tau e^{i\omega\tau}$
$\displaystyle\langle
a|u_{x}(-\textbf{k},\tau)H-Hu_{x}(-\textbf{k},\tau)|\beta\rangle$
$\displaystyle\times\langle\beta|u_{y}(\textbf{k},0)H-Hu_{y}(\textbf{k},0)|a\rangle$
$\displaystyle=\left(-\frac{i}{\hbar}\right)^{2}\int d\tau e^{i\omega\tau}$
$\displaystyle\langle
a|u_{x}(-\textbf{k},\tau)|\beta\rangle\langle\beta|u_{y}(\textbf{k},0)|a\rangle$
$\displaystyle\times$
$\displaystyle\left[-\left(E_{a}-E_{\beta}\right)^{2}\right]$
$\displaystyle=-\left[\frac{i(E_{a}-E_{\beta})}{\hbar}\right]^{2}\int d\tau
e^{i\omega\tau}$ $\displaystyle\langle
a|e^{i\frac{H}{\hbar}\tau}\textbf{u}_{x}(-\textbf{k})e^{-i\frac{H}{\hbar}\tau}|\beta\rangle$
$\displaystyle\times$ $\displaystyle\langle\beta|u_{y}(\textbf{k},0)|a\rangle$
$\displaystyle=\left(\frac{E_{\beta}-E_{a}}{\hbar}\right)^{2}$
$\displaystyle\langle
a|\textbf{u}_{x}(-\textbf{k})|\beta\rangle\langle\beta|u_{y}(\textbf{k})|a\rangle$
$\displaystyle\times$ $\displaystyle\int d\tau
e^{i\left(\omega-\frac{E_{\beta}-E_{a}}{\hbar}\right)}$ (21)
$\displaystyle=\omega^{2}\langle
a|\textbf{u}_{x}(-\textbf{k})|\beta\rangle\langle\beta|u_{y}(\textbf{k})$
$\displaystyle|a\rangle$ $\displaystyle\times$
$\displaystyle\delta\left(\omega-\frac{E_{\beta}-E_{a}}{\hbar}\right)$ (22)

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

$\omega^{2}\int d\tau e^{i\omega\tau}\langle
a|u_{x}(-\textbf{k},\tau)|\beta\rangle\langle\beta|u_{y}(\textbf{k},0)|a\rangle$
(23)

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

$\displaystyle Z_{xy}(\textbf{k},\omega)$ $\displaystyle=\omega^{2}\int d\tau
e^{i\omega\tau}\langle u_{x}(-\textbf{k},\tau)u_{y}(\textbf{k},0)\rangle$
$\displaystyle=\omega^{2}B_{xy}(\textbf{k},\omega)$ (24)

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\rightarrow 0$ [81]. By applying
this limit to the expression for $B_{xy}(\textbf{k},\omega)$ in eq. (18), we
obtain for the prefactor

$\lim_{\hbar\rightarrow
0}\frac{\hbar}{1-e^{-\beta\hbar\omega}}=\lim_{\hbar\rightarrow
0}\frac{\hbar}{1-\left(1-\beta\hbar\omega\right)}=\frac{1}{\beta\omega}$ (25)

Another important point to be mentioned is the behaviour of the self energy
components $\Delta_{\textbf{k}j}(\omega)$ and $\Gamma_{\textbf{k}j}(\omega)$
in the limit $\hbar\rightarrow 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\rightarrow\infty$ and is found to be finite and non-zero. Therefore,
combining eqs (18), (24) and (25), we can finally write

$\displaystyle Z_{xy}(\textbf{k},\omega)$ $\displaystyle=\frac{\omega
k_{B}T}{Nm}\sum_{j}\left[\frac{e^{*}_{x}(\textbf{k},j)e_{y}(\textbf{k},j)}{\omega(\textbf{k},j)}\right]$
$\displaystyle\times\frac{\Gamma_{\textbf{k}j}}{\left(\omega-\omega(\textbf{k},j)-\Delta_{\textbf{k}j}\right)^{2}+\Gamma^{2}_{\textbf{k}j}}$
(26)

which is not exactly Lorentzian, as it has a prefactor that depends on
$\omega$. However, as per the assumption that
$\Delta_{\textbf{k}j},\Gamma_{\textbf{k}j}\ll\omega_{\textbf{k}j}$, 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_{\textbf{k},j}$, finally getting

$Z_{xy}(\textbf{k},\omega)=\frac{k_{B}T}{NM}\sum_{j}\frac{\Gamma_{\textbf{k}j}e^{*}_{x}(\textbf{k},j)e_{y}(\textbf{k},j)}{\left(\omega-\omega(\textbf{k},j)-\Delta_{\textbf{k}j}\right)^{2}+\Gamma^{2}_{\textbf{k}j}}$
(27)

which is now a superposition of Lorentzians, with positions
$\left(\omega(\textbf{k},j)+\Delta_{\textbf{k}j}\right)$ and FWHM
$2\Gamma_{\textbf{k}j}$ 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\rightarrow m_{p}$ and
$e_{x}(\textbf{k},j)\rightarrow e_{x}^{p}(\textbf{k},j)$, where now $p$ is an
index to atom species. Then, by remembering that for each phonon mode
$(\textbf{k},j)$ the eigenvectors are orthonormal
$\sum_{x,p}e_{x}^{p}(\textbf{k},j)e_{x}^{*p}(\textbf{k},j)=1$, we can write
$Z_{total}(\textbf{k},j)$ as

$\displaystyle
Z_{total}(\textbf{k},\omega)=\sum_{x,p}m_{p}Z_{xx}(\textbf{k},\omega)$
$\displaystyle=\frac{k_{B}T}{N}\sum_{j}\frac{\Gamma_{\textbf{k}j}}{\left(\omega-\omega(\textbf{k},j)-\Delta_{\textbf{k}j}\right)^{2}+\Gamma^{2}_{\textbf{k}j}}$
(28)

where we have used the correspondence between eqs (9) and (16) as
$Z_{x}(\textbf{k},\omega)=Z_{xx}(\textbf{k},\omega)$

## References

  * Urade  _et al._ [2023] A. R. Urade, I. Lahiri, and K. S. Suresh, Graphene properties, synthesis and applications: A review, JOM 75, 614 (2023).
  * Novoselov  _et al._ [2004] 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).
  * Wang  _et al._ [2019] 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).
  * Lee  _et al._ [2008] 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).
  * Sang  _et al._ [2019] 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/nano9030374 (2019).
  * Novoselov  _et al._ [2005] 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).
  * Murata  _et al._ [2019] 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).
  * Singh  _et al._ [2011] 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).
  * Liu  _et al._ [2016] 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).
  * Hamze  _et al._ [2020] 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).
  * Fu  _et al._ [2019] 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).
  * Androulidakis  _et al._ [2020] 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).
  * Shahil and Balandin [2012] 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.
  * Cao  _et al._ [2018] 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).
  * Davis  _et al._ [2023] 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).
  * Nika  _et al._ [2011] 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).
  * Ghosh  _et al._ [2008] 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).
  * Lindsay  _et al._ [2011] L. Lindsay, D. A. Broido, and N. Mingo, Flexural phonons and thermal transport in multilayer graphene and graphite, Phys. Rev. B 83, 235428 (2011).
  * Nemanich  _et al._ [1977] R. Nemanich, G. Lucovsky, and S. Solin, Infrared active optical vibrations of graphite, Solid State Communications 23, 117 (1977).
  * Davydov [1971] A. S. Davydov,  _Theory of Molecular Excitons_ (Springer New York, NY, 1971).
  * Ghosh [1976] P. N. Ghosh, Davydov splitting and multipole interactions, Solid State Communications 19, 639 (1976).
  * Molina-Sánchez  _et al._ [2015] 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).
  * Gibaja  _et al._ [2016] C. Gibaja, D. Rodriguez-San-Miguel, P. Ares, J. Gómez-Herrero, M. Varela, R. Gillen, J. Maultzsch, F. Hauke, A. Hirsch, G. Abellán, and F. Zamora, Few-layer antimonene by liquid-phase exfoliation, Angewandte Chemie International Edition 55, 14345 (2016).
  * Saha  _et al._ [2008] 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).
  * Maradudin and Fein [1962] A. A. Maradudin and A. E. Fein, Scattering of neutrons by an anharmonic crystal, Phys. Rev. 128, 2589 (1962).
  * Hellman  _et al._ [2013] 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).
  * Gu  _et al._ [2019] X. Gu, Z. Fan, H. Bao, and C. Y. Zhao, Revisiting phonon-phonon scattering in single-layer graphene, Phys. Rev. B 100, 064306 (2019).
  * Maultzsch  _et al._ [2004] J. Maultzsch, S. Reich, C. Thomsen, H. Requardt, and P. Ordejón, Phonon dispersion in graphite, Phys. Rev. Lett. 92, 075501 (2004).
  * Sullivan  _et al._ [2017] 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).
  * Calizo  _et al._ [2007] 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).
  * Han  _et al._ [2022] 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).
  * Ferrari  _et al._ [2006] 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).
  * Late  _et al._ [2011] 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).
  * Tan  _et al._ [2012] 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).
  * Lui  _et al._ [2012] 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.
  * Lui  _et al._ [2014] 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.
  * Lin  _et al._ [2011] J. Lin, L. Guo, Q. Huang, Y. Jia, K. Li, X. Lai, and X. Chen, Anharmonic phonon effects in raman spectra of unsupported vertical graphene sheets, Phys. Rev. B 83, 125430 (2011).
  * Katsiaounis  _et al._ [2021] 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).
  * Kang  _et al._ [2010] 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).
  * Chatzakis  _et al._ [2011] 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).
  * Wang  _et al._ [2010] 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).
  * Yan  _et al._ [2009] 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).
  * Kumar  _et al._ [2009] 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).
  * Bonini  _et al._ [2007] N. Bonini, M. Lazzeri, N. Marzari, and F. Mauri, Phonon anharmonicities in graphite and graphene, Phys. Rev. Lett. 99, 176802 (2007).
  * Koukaras  _et al._ [2015] E. N. Koukaras, G. Kalosakas, C. Galiotis, and K. Papagelis, Phonon properties of graphene derived from molecular dynamics simulations, Scientific Reports 5, 12923 (2015).
  * Noid  _et al._ [2008] 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).
  * Rahman [1964] A. Rahman, Correlations in the motion of atoms in liquid argon, Phys. Rev. 136, A405 (1964).
  * Kohanoff [1994] J. Kohanoff, Phonon spectra from short non-thermally equilibrated molecular dynamics simulations, Computational Materials Science 2, 221 (1994).
  * Wang  _et al._ [1989] C. Z. Wang, C. T. Chan, and K. M. Ho, Molecular-dynamics study of anharmonic effects in silicon, Phys. Rev. B 40, 3390 (1989).
  * Papanicolaou  _et al._ [1995] 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).
  * Ditlevsen  _et al._ [1991] P. D. Ditlevsen, P. Stoltze, and J. K. No/rskov, Anharmonicity and disorder on the cu(110) surface, Phys. Rev. B 44, 13002 (1991).
  * Thomas  _et al._ [2010] 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).
  * Mafakheri and Dodaran [2021] 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).
  * Chen and Zhang [2019] 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).
  * Desmarchelier  _et al._ [2021] P. Desmarchelier, A. Tanguy, and K. Termentzidis, Thermal rectification in asymmetric two-phase nanowires, Phys. Rev. B 103, 014202 (2021).
  * Zhang  _et al._ [2022] 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).
  * Chen  _et al._ [2020] 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.
  * Born and Huang [1956] M. Born and K. Huang,  _Dynamical theory of crystal lattices_ , International series of monographs on physics (Clarendon Press Oxford, Oxford, 1956).
  * Adams [2001] 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.
  * Kolmogorov and Crespi [2005] A. N. Kolmogorov and V. H. Crespi, Registry-dependent interlayer potential for graphitic systems, Phys. Rev. B 71, 235415 (2005).
  * Leven  _et al._ [2014] I. Leven, I. Azuri, L. Kronik, and O. Hod, Inter-layer potential for hexagonal boron nitride, The Journal of Chemical Physics 140, 104106 (2014).
  * Wen  _et al._ [2018] 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).
  * Lebedeva  _et al._ [2012] 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 (2012), the proceedings of the European Materials Research Symposium on Science and Technology of Nanotubes, Nanowires and Graphene.
  * Mohr  _et al._ [2007] 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).
  * Al-Qasir  _et al._ [2020] 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).
  * Lindsay and Broido [2010] 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).
  * Rowe  _et al._ [2018] P. Rowe, G. Csányi, D. Alfè, and A. Michaelides, Development of a machine learning potential for graphene, Phys. Rev. B 97, 054303 (2018).
  * Ouyang  _et al._ [2018] W. Ouyang, D. Mandelli, M. Urbakh, and O. Hod, Nanoserpents: Graphene nanoribbon motion on two-dimensional hexagonal materials, Nano Letters 18, 6009 (2018).
  * Gould  _et al._ [2013] T. Gould, Z. Liu, J. Z. Liu, J. F. Dobson, Q. Zheng, and S. Lebègue, Binding and interlayer force in the near-contact region of two graphite slabs: Experiment and theory, The Journal of Chemical Physics 139, 224704 (2013).
  * Wang  _et al._ [2015] 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.
  * Thompson  _et al._ [2022] 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).
  * Pizzi  _et al._ [2021] G. Pizzi, S. Milana, A. C. Ferrari, N. Marzari, and M. Gibertini, Shear and breathing modes of layered materials, ACS Nano 15, 12509 (2021).
  * Cong and Yu [2014] C. Cong and T. Yu, Enhanced ultra-low-frequency interlayer shear modes in folded graphene layers, Nature Communications 5, 4709 (2014).
  * Yang  _et al._ [2020] 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).
  * Liu  _et al._ [2019] 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).
  * Balkanski  _et al._ [1983] 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).
  * Blackman and Tukey [1959] 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).
  * Chakraborty  _et al._ [2018] 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).
  * Mahan [1990] G. Mahan,  _Many-Particle Physics_ , Physics of Solids and Liquids (Springer US, 1990).
  * Zwanzig [1965] R. Zwanzig, Time-correlation functions and transport coefficients in statistical mechanics, Ann. Rev. Phys. Chem. 16, 67 (1965).