Free energy barrier and thermal-quantum behavior of sliding bilayer graphene

Abstract

In multilayer graphene, the stacking order of the layers plays a crucial role in the electronic properties and the manifestation of superconductivity. By applying shear stress, it is possible to induce sliding between different layers, altering the stacking order. Here, focusing on bilayer graphene, we analyze how ionic fluctuations alter the free energy barrier between different stacking equilibria. We calculate the free energy barrier through the state-of-the-art self-consistent harmonic approximation, which can be evaluated at unstable configurations. We find that above 100 K there is a large reduction of the barrier of more than 30% due to thermal vibrations, which significantly improves the agreement between previous first-principles calculations and experiments on a single graphite crystal. As the temperature increases, the barrier remains nearly constant up to around 500 K, with a more pronounced decrease only at higher temperatures. Our approach is general and paves the way for systematically accounting for thermal effects in free energy barriers of other macroscopic systems.

---

I. Introduction

There has been intense research on graphene since its discovery in 2004, and also on multilayer graphene systems. Among many other interesting properties, flat bands are present in multilayer rhombohedral (ABC-stacked) graphene (RG),¹ and superconductivity has been observed in RG in 2021.² Several works recently studied, both experimentally and theoretically, under which conditions RG can form,^3–8 given that AB-stacked graphene is the more commonly occurring phase. In particular, some of us unveiled how the stacking of graphene layers changes under shear stress,⁹ and produces many consecutive layers of RG. In fact, all samples with long-range RG that have been observed so far have been subject to some form of shear stress.^5,10–14

To move graphene layers with respect to each other, they have to overcome a free energy barrier that separates different minima. When doing so, they exhibit “stick–slip” motion: they move gradually with increasing shear stress, followed by sudden rearrangements and a decrease of the stress.¹⁵ The maximum shear–stress required to “unlock” layers is correlated with the height of the barrier. Such a barrier corresponds to the activation energy of a metastable state, where a graphene flake sitting on top of a large graphene layer is displaced relative to the equilibrium position.

Superlubricity is another phenomenon where the energy landscape corresponding to the relative displacement of layers is also central. In this mechanism, which has attracted a lot of attention,^16,17 particularly in carbon-based systems like multilayer graphene and graphite, incommensurate surfaces (i.e. atomic lattices without matching periodicity) lead to extremely low friction. Yet another phenomenon involving sliding surfaces is thermolubricity, which pertains to the assistance from thermal excitations to overcome barriers. This term has been coined by some works that have studied the transition from stick–slip motion into stochastic fluctuations in friction at sufficiently high temperatures.¹⁸ Despite the interest in sliding surfaces, there are no studies that account for both temperature and quantum effects of free-energy barriers in macroscopic systems; instead, existing works have focused on systems involving hydrogen molecules. For example, zero-point renormalization of the free energy barrier or tunneling have been examined in the dissociative absorption¹⁹ or diffusion of hydrogen.^20–22

The variation of the potential along a periodic crystal surface, combined with the attraction between layers, prevents free sliding between them. This phenomenon has been extensively studied in various systems, like graphite,²³ MoS₂,²⁴ silicon on gold,²⁵ among others. When using a tip with the goal of making layers slide, normal forces need not be large. For pressures much lower than 1 GPa, the alteration of the free-energy landscape is minimal (see the SI of ref. 9), so the relaxed barrier remains the dominant contribution that impedes sliding. In addition, in some works on superlubricity it has been observed that friction has an extremely weak dependence with the load, attributed for example²³ to a contact area that remains nearly constant in the force range applied. Here, we focus on calculating the height of the barrier, which should be approximately proportional to the force – given by the maximum slope of the free energy profile – that is needed for a layer to move from one minimum to another in commensurate systems.

A quantity closely related to the free energy barrier is the shear frequency, corresponding to the mode in which layers move rigidly in-plane and out of phase around the equilibrium position. The free energy barrier separates two equivalent minima, so the barrier and shear frequency are expected to have a similar temperature dependence (a lower curvature at the minima corresponds to a lower shear frequency and a lower free energy, at least close to the minima). The value of the shear frequency in graphite is a relatively well-established quantity, with experimental values²⁶ oscillating between 42 and 45 cm⁻¹, using for example inelastic X-ray scattering,²⁷ neutron-coherent inelastic scattering,²⁸ Raman spectra,^29,30 and coherent phonon spectroscopy.³¹ In bilayer graphene, there are some experiments with values around 30 cm⁻¹, for example, 28 ± 3 cm⁻¹ in ref. 31 and 32, and 32 cm⁻¹ in ref. 29. The temperature dependence of the shear modes, however, has been barely studied. Few experiments measure the temperature dependence of the shear mode, using Raman scattering in folded 2+2 (4 layer) graphene,³³ or femtosecond pump–probe spectroscopy in graphite.³⁴

Another aspect of multilayer graphene and graphite that still needs further investigation and that serves as a consistency check on the adequacy of our approach is thermal expansion (TE). In ref. 35 it was measured in bilayer and trilayer graphene for the first time. In ref. 36 the quasi-harmonic approximation (QHA) is used, which underestimates the out-of-plane TE at low temperatures³⁷ and exhibits a very rapid increase at high temperatures; in-plane, it also underestimates the TE at low temperatures³⁷ and overestimates it at high temperatures.³⁸ Much better out-of-plane results were obtained by using path integral molecular dynamics (PIMD),³⁹ but while the TE flattens at 500 K in graphite, experimental values continue to increase.^38,40 However, in-plane PIMD results of ref. 39 vastly overestimate the TE.

Motivated by the aforementioned phenomena, especially by the possibility of facilitating sliding through temperature, we focus on thermal effects on the free-energy barrier. This is a very challenging problem due to the difficulty of determining the free energy in a non-equilibrium position. Here we tackle it using the state-of-the-art stochastic self-consistent harmonic approximation (SSCHA).^41,42 First, in Section II, we briefly introduce SSCHA and the QHA, and explain how the barrier is determined. Then, in Section III, we present our main results and compare them to the height of the barrier with fixed atoms. Subsequently, also using the QHA, we calculate the temperature dependence of the shear frequency and compare it to the limited experimental data available. Finally, we determine the in-plane and out-of-plane TE coefficient in bilayer graphene and graphite, and compare our results to previous theoretical and experimental work. In Section IV, we summarize the conclusions of our work.

II. Theoretical and computational framework

Here, we briefly describe the QHA, SSCHA, and a new SSCHA interpolation method used to determine the barrier. A more detailed description of SSCHA can be found in ref. 41 or 42.

A. Quasi-harmonic approximation (QHA)

In the QHA, at each value of the lattice parameters {a_i}, the free energy F(T,{a_i}) = U({a_i}) − TS(T,{a_i}), where U is the internal energy, S the entropy and T the temperature, is given by the standard harmonic expression. That is,where U_lat{a_i} is the internal energy of the lattice for a given value of the lattice parameters, N_q the number of phonon wavevectors q, ν the mode, and ω_qν the corresponding phonon frequency.

The values of the lattice parameters at a given temperature, a_i(T), are determined by minimizing the free energy. This also leads to temperature dependent phonon frequencies ω_qν({a_i(T)}). Using an interatomic potential, the phonon frequencies at each value of the lattice parameters can be determined with little computational effort, and thus the QHA offers a quick method to determine the lattice expansion. Using first-principles calculations, phonons can be determined using density functional perturbation theory, which uses a primitive cell and is much faster than frozen phonon (finite difference) calculations in supercells. However, it does not account for the “true anharmonicity”,⁴³ present at fixed values of the parameters when varying temperature. For example, it gives an increasing as opposed to decreasing temperature dependence of the G mode in graphene, due to the negative TE of the lattice parameter. Also, if the TE is negative and the lattice parameter reduces below the classical lattice value, as in the case of monolayer graphene at about 400 K (see the SI, Fig. S12), then an acoustic mode becomes unstable and the QHA is ill-defined. These issues are not present in SSCHA.

B. Stochastic self-consistent harmonic approximation (SSCHA)

In general, the free energy of an ionic Hamiltonian H = T + V is given bywhere ρ_H is the density matrix. In SSCHA, the density matrix is restricted to a trial harmonic Hamiltonian , and the free energy is given by

Due to a variational principle, . The parameters of the harmonic Hamiltonian include the centroids of the atoms (their average positions), which correspond to the lattice parameters a_i of the QHA, but it also includes the interatomic force constant matrix, which contains the information of all the phonon frequencies and accounts for the true anharmonicities. Thus, SSCHA has additional parameters, and is expected in general to give better results than QHA. Also, QHA does not satisfy a variational principle. It can be shown that

where R indicates a general ionic configuration and is the probability density of finding the crystal in a generic configuration is the partition function). The second term includes the difference between the true potential V and the trial harmonic potential . By stochastically sampling configurations R for a given density matrix , and determining the gradient of the free energy with respect the harmonic parameters (centroids and force constant matrix elements), the free energy can be minimized self-consistently. SSCHA is thus a nonperturbative method, as it does not rely on a perturbative expansion: it explores configurations with amplitudes that depend on temperature, without being limited to small displacements.

Thanks to the variational ansatz, the SSCHA quantum free energy is rigorously defined, even in out-of-equilibrium and unstable configurations like SP. In contrast, the QHA cannot determine the barrier, since the harmonic frequencies are not well defined at SP.

The logarithm term in the harmonic term eqn (1) can be hard to converge for low frequencies. Separating the total free energy as in eqn (4), the harmonic term can be interpolated into a very fine q-mesh, which is fundamental to avoid oversampling the energy of the long-ranged shear modes, while the anharmonic contribution can be properly converged already with small simulation cells. So rather than calculating directly with the grid corresponding to the SC used in the minimization procedure, the force constant matrix can be Fourier interpolated in the usual way (building the force constant in real space from the existing q-grid, and then Fourier transforming back to any desired value of q), and can be calculated using a dense mesh.

This procedure is followed both at SP and the equilibrium position, and the barrier is obtained by taking the difference of the corresponding free energies. More technical details are provided in the SI, Section S1.

The SSCHA frequencies are in principle auxiliary frequencies that minimize the free energy. In a purely harmonic system, such frequencies correspond to the physical frequencies. But in general, they have to be corrected with the so-called “bubble correction” (and also an additional higher order term in very anharmonic cases),⁴² which we do in Section III B to obtain the temperature dependence of the shear frequency. Convergence with the supercell size of both the physical and auxiliary frequencies can be seen in the SI, Fig. S7.

Since the centroids of the atoms are also parameters in SSCHA, from the minimization of the free energy we additionally obtain the lattice parameters as a function of temperature, allowing us to study their thermal expansion. This is discussed in Section III C.

III. Results

III.A. Barrier

We first focus on studying bilayer graphene, where the in-plane and (most importantly) out-of-plane interactions are given by the interatomic potential of ref. 44 and 45, respectively. The out-of-plane interactions determine the barrier. The energy profile of a layer moving with respect to the other along the bond direction, shown in Fig. 1, is obtained by fixing the in-plane positions of the layers and relaxing the interlayer distance until the out-of-plane component of the atomic forces vanish, and gives good agreement with theory and experiment (see Table 1).

Fig. 1 Energy profile of bilayer graphene when moving one layer relative to the other along the bond direction. The lower layer is fixed in position A. Configuration AA (one layer on top of each other) is the least favorable configuration. There are two minima of the same energy, AB and AC, which correspond to the two stable configurations. To go from one minima to the other one applying shear stress, the upper layer has to overcome the barrier V_SP, which has a value of 1.24 meV per atom relative to the minima. The diagrams below illustrate the different configurations, with the upper layer (darker blue) moving relative to the lower one in position A (lighter blue).

Table 1 Comparison of the barrier V_SP and the shear frequency at fixed ionic positions (classical at T = 0 K) between the interatomic potential we used in this work,^44,45 LDA, and theoretical calculations and experiments from previous works. The values are a little bit lower but close to that of the adiabatic-connection fluctuation–dissipation theorem within the random phase approximation (ACFDT-RPA)⁴⁶ and experimental data²⁶

	2 layers this work	2 layers LDA	Previous works
Small barrier VSP (meV per atom)	1.24	1.58	1.53 (RPA)^46
Shear frequency (cm^−1) (bulk)	37	42	42–45^26

---

The lower layer is fixed in position A, and there are two minima of the same energy, AB and AC, at a distance of a bond length d = 1.42 Å. For layers to slide with respect to each other, they have to overcome a barrier V_SP. It is worth noting that although the profile is for bilayer graphene, some of us previously showed in first-principles calculations that using 6 layers – 3 layers displaced relative to the other 3 – the profile is almost identical (see the SI of ref. 9). Since interactions with subsequent layers are even weaker, the barrier should stay essentially unchanged with additional layers. In addition, in ref. 47 we obtained that energy differences between different functionals for different stacking sequences, including van-der-Waals corrections, differ by around 0.01 meV per atom. Thus, the barrier in the bulk limit should correspond closely to the bilayer one.

SP is a saddle point in the full energy landscape, so the shear mode in the perpendicular direction (along the direction of two nearby AA maxima) is stable. In order to determine temperature effects on the barrier, we use SSCHA to determine the free energy at the equilibrium position and SP, using the interpolation method described in Section II B. The barrier as a function of temperature can be seen in blue in Fig. 2.

Fig. 2 Barrier as a function of temperature. It is reduced from 1.24 meV per atom, the value at fixed nuclei (V_SP in Fig. 1), to about 0.8 meV per atom, corresponding to a large reduction of about 35% due to phonons. The barrier does not change much up to about 500 K, but the value at 1000 K is about 20% smaller relative to the value at 100 K. The lines are a guide to the eye, as in the other figures of this work. The standard quantum SSCHA result (blue) and its classical limit (red) give almost identical results, showing quantum effects are negligible above 100 K.

The barrier at 100 K is reduced by a substantial 35% with respect to the value at fixed nuclei. We expect intralayer modes to contribute little to the reduction of the barrier, since they should not change much when displacing the layers. On the other hand, modes involving the relative motion of the layers should be much more affected. Indeed, intralayer modes are virtually identical, so they do not contribute to a change in the barrier. Regarding the interlayer modes, the breathing mode changes little, while the shear modes do differ significantly (see Fig. S4). Thus, the shear modes are responsible for the large reduction of the barrier. This makes sense intuitively: modes that involve a relative in-plane displacement of the layers should facilitate sliding. In the fixed ion calculation, the system is exactly at SP, but in reality, the system is also exploring surrounding configurations, which have a lower energy. In an analogous fashion, there is a barrier in real space, since the interlayer distance at AB is lower relative to SP. At 100 K, this barrier gets similarly reduced by 38%, from 0.029 Å to 0.018 Å, and also shows a temperature dependence similar to that of the free energy barrier (see Fig. S5 and S6).

Remarkably, our results show improved agreement with experiments. In a previous work,⁹ using density functional theory with an LDA functional (see Table 1), some of us obtained that the maximum shear stress in simulated stick–slip motion is 0.2 GPa, and we mentioned a good agreement with experimental values 0.14 GPa (ref. 15) – a benchmark for defect free single crystal graphite⁴⁸ – and 0.1 GPa (ref. 23). Assuming that the maximum shear stress decreases by 35% in line with the barrier, we now obtain 0.13 GPa, in excellent agreement with experiments.^15,23

Since 100 K – about 70 cm⁻¹ in wavenumber units – is a relatively low temperature, quantum fluctuations may play a role. To assess their impact on the free energy barrier, we also performed calculations with the classical distribution, corresponding to the classical limit ħω/kT → 0. In this way, 2n_s + 1 → 2T/ω, and the quantum (Bose–Einstein) probability distribution for ionic configurations reduces to the classical Boltzmann distribution (an explicit expression for the distribution is included in the SI). We can see the SSCHA results in the classical limit in red in Fig. 2, which are virtually identical to the SSCHA quantum calculations, showing that quantum effects are negligible above 100 K. The reason is that shear modes responsible for the barrier reduction have very low energies, and even a modest temperature can populate such modes. This is not a general feature, as other systems present significant differences between classical and quantum free energy barriers, as in materials formed by light hydrogen atoms. Such quantum effects have been observed at 100 K (ref. 19), 250 K (ref. 21), and even higher temperatures.²⁰

The size of the graphene flakes considered in our work in principle corresponds from nanometers up to possibly micrometer scales. In clean samples with micrometer sized flakes, different regions can have different stackings (AB and AC),^12,49 and a domain wall in the metastable state separates both regions. The displacement of layers in this case is more complex and is not considered in our work. We expect shear stress to be lower in this case, since the metastable state is already present. Defects and stacking faults also lead to lower friction due to incommensurabilities, so the shear stress of 0.14 GPa corresponds to the maximum expected value.⁴⁸

When it comes to displacing many layers in multilayer graphene or graphite, the change of the barrier could in principle differ, since it depends on the phonon frequencies (and interlayer modes depend noticeably on the number of layers). For example, graphite exhibits a larger shear frequency at equilibrium, since a layer experiences primarily the restoring force of two layers instead of one (see the next subsection for more details). However, layers are typically not simultaneously in unstable positions relative to both nearest layers, so, as noted earlier, the barrier remains approximately unchanged in bulk. Since the frequency at AB increases while remaining similar at SP, the barrier reduction in bulk might be even greater than in the bilayer case. Other variations may arise from the choice of the interlayer interatomic potential: as shown in Table 1, the RPA barrier is higher, indicating that the curvature around SP could be larger. Thus, vibrations around SP with more accurate interlayer potentials might yield a larger absolute barrier reduction, although the relative change may still be similar.

In other materials, the anharmonicities associated with the unstable modes at a saddle point should also contribute negatively to the free energy. However, the overall change of the barrier will depend on differences in the entire SSCHA phonon spectra at the minima and saddle point.

Regarding the evolution of the barrier as the temperature increases, there is a weak dependence up to about 500 K, since the barrier remains close to around 0.8 meV per atom. However, when temperature increases to 1000 K, the barrier is further reduced by 20% relative to the value at 100 K. Thus, high temperatures should facilitate sliding between layers. In particular, it could aid the production of ABC-stacked multilayer graphene from its more common AB-stacked version through shear stress, as mentioned in the Introduction. More in general, high temperatures could also help engineer specific stacking arrangements in other layered materials.

III.B. Shear frequency

At the equilibrium position, one of the vibrational modes at Γ is the shear mode in both directions of the plane (doubly degenerate, see Fig. 3). Its temperature dependence can be determined using both SSCHA and QHA, as described in Section II.

Fig. 3 Bilayer graphene phonon dispersion. The inset shows a zoomed in image close to Γ to better visualize the doubly degenerate shear modes S at Γ, and the breathing mode B. Due to the high speed velocity in graphene, the region of the shear mode in the BZ is small (a fraction of ΓM), and is overrepresented in small supercell calculations. The breathing mode, on the other hand, is flat, so nearby sampling points are not necessary to converge the free energy.

The third out of the phase modes of the layers is the breathing mode, in the out-of-plane direction. A discussion of its temperature dependence can be found in the SI.

The temperature dependence of the shear frequency can be seen in Fig. 4. The quantum effects evaluated within the SSCHA harden the shear frequency energy at low temperatures compared with the classical value at 0 K (dashed line). On the other hand, the zero-point renormalization in QHA softens the shear frequency. The temperature dependence of the shear mode is similar to that of the barrier, varying by about 20% from 0 to 1000 K, and having a larger dependence at temperatures above 500 K.

Fig. 4 Shear frequency as a function of temperature using SSCHA (blue circles) and QHA (red crosses). There is no experimental data for the temperature dependence of the shear frequency as a function of temperature for bilayer graphene. Experimental values correspond to bulk (blue crosses) and folded 2+2 graphene (orange circles), which are corrected by the factor of a nearest neighbor model⁵⁰ to compare to the bilayer case (and rigidly shifted by −5 cm⁻¹ to facilitate the visual comparison). The agreement with the temperature dependence of the bulk values is good.

Regarding the experimental data, in ref. 33 the measurements are done in 2+2 folded graphene (4 layers), and in ref. 34 they are done in the bulk. To compare with our bilayer calculations, we rely on a nearest layer model which has simple conversion factors between the different number of layers.²⁹ Since interlayer interactions in graphite are known to be weak, this should provide a good approximation. Indeed, there are multiple references that have measured the shear frequency for different number of layers and obtained that the model describes the trend very well, like in multilayer graphene²⁹ and NbSe₂.⁵⁰ The bulk case can be readily understood: while in the bilayer case the restitutive lateral force is exerted only by the other layer, in the bulk case each layer is surrounded by two of them, so the force is doubled. Since the interatomic force constant matrix is proportional to the square of the phonon frequencies, the shear frequency in bilayer is smaller than the bulk frequency in this model. Using these factors,²⁹ we obtain the adjusted experimental values of Fig. 4.

The temperature dependence of both SSCHA and QHA curves is similar to the experimental adjusted bulk values. Folded graphene values have a steeper dependence but are likely unsuitable for comparison, because layers fold with different rotational angles, while layers are aligned in standard graphite. The shear frequency has also been measured in other materials like bilayer NbSe₂ in ref. 50, or bulk h-BN in ref. 51, which show a similar temperature dependence (adjusting the bulk value by a factor) of roughly −0.5 cm⁻¹/100 K.

III.C. Thermal expansion

To further check the adequacy of our method (the use of SSCHA with the in-plane⁴⁴ and out-of-plane⁴⁵ interatomic potentials), we look at the TE in the previous calculations and also in bulk graphite, and compare to experiments. The TE coefficients are given by and , with a and c the in-plane and out-of-plane lattice parameters, respectively.

In Fig. 5, we show the TE coefficient calculated with SSCHA and QHA in the bulk, and compare with experimental data^37,40 and PIMD calculations.³⁹ SSCHA results underestimate the TE, which we attribute to the interlayer potential.⁴⁵ This potential is fitted to energies and forces that use the many-body dispersion method in conjunction with a PBE functional, as it reproduces binding energies and interlayer distances reasonably well. However, it yields an elastic modulus along the c-axis about 15% lower than experimental values (see more details within ref. 45). Unsurprisingly, ref. 39 obtains better results with the LCBOPII potential⁵² – although the TE remains somewhat below experiments^37,40 – since it uses the experimental out of plane compressibility as a fitting parameter.⁵³ The bilayer TE is also similarly underestimated (see SI). Thus, the variation of the potential around the interlayer equilibrium distance seems to be the main factor influencing the TE, as opposed to the weaker long-range interactions involving subsequent layers. It is also possible that anharmonicities, captured more accurately by PIMD relative to SSCHA, play a role at high temperatures. In any case, for the purpose of studying the barrier – the focus of our work – the problem with the LCBOPII potential is that it includes no registry effects, making the barrier essentially non-existent.

Fig. 5 (a) Graphite out-of-plane coefficient of TE as a function of temperature using SSCHA and QHA. PIMD calculations³⁹ and experimental data³⁷ are also included. PIMD overall agrees well, while QHA and SSCHA (with the potential of ref. 45) underestimates the TE. (b) Graphite in-plane coefficient of TE using SSCHA and QHA. QHA again underestimates the dependence at low temperatures, and also at high temperatures. On the other hand, PIMD gives a dependence that is surprisingly very far from experimental values. SSCHA gives excellent results both at low and high temperatures, considering also the disagreement within experimental data.

While the out-of-plane TE we obtained is not optimal, in-plane SSCHA calculations give excellent results. QHA works well at low temperatures but departs from experiments at higher temperatures. The negative TE in graphite and other layered systems is well-known and arises from flexural (out-of-plane) modes.⁵⁴ We briefly summarize the mechanism here. In flexural modes, atoms are displaced out-of-plane, reducing the average projected in-plane distance, even though the 3D bond lengths remain essentially unchanged. This can also be understood intuitively via the “membrane effect”: when a layer is stretched, the out-of-plane amplitude decreases and the frequency increases, which corresponds to a negative Grüneisen parameter. See, for example, ref. 54 for more details. We conducted a careful convergence of these results, reported in the SI, Section S3. PIMD calculations of ref. 39, although they work very well out-of-plane, greatly overestimate in-plane TE, presumably due to the inadequacy of their interatomic potential. It would be interesting in future work to compare SSCHA and PIMD using the same potentials.

IV. Conclusions

To conclude, we demonstrated that ionic vibrations play a fundamental role in determining properties of layered materials that depend on free energy barriers. Using the self-consistent harmonic approximation, which can calculate the free energy at unstable configurations, and a novel interpolation method, we unveiled the suppression of the barrier separating stable configurations in bilayer graphene. Phonon vibrations reduce the free energy barrier by a significant 35%, and get further suppressed at high temperatures, so barriers at fixed ions cannot be trusted. Now, correcting the shear stress of 0.2 GPa reported in our previous work, we arrive at a value of 0.13 GPa, in excellent agreement with experiments. The shear modes, which occupy a small region of the phase space, have to be adequately sampled to correctly determine the free energy barrier as well as the thermal expansion. Although our approach underestimates the out-of-plane thermal expansion compared with experimental values, in-plane agreement with experiments is excellent. We also observed that the temperature dependence of the barrier is similar to that of the shear frequency, which agrees well with existing measurements.

Conflicts of interest

There are no conflicts to declare.

Data availability

Additional data supporting this article has been included as part of the supplementary information (SI). Supplementary information is available. See DOI: https://doi.org/10.1039/d5cp02188j.

The SSCHA⁴² code used in this article is publicly available at https://sscha.eu/.

Acknowledgements

This work received funding from the European Union's Horizon 2020 research and innovation programme under Grant Agreement No. 881603 and from the MORE-TEM ERC-SYN project, Grant Agreement No. 951215. J. P. N. is currently supported by the European Union under a Horizon Europe Marie Skłodowska-Curie Individual Fellowship, Project GreenNP No. 101151380.

References

1. M. Koshino, Phys. Rev. B: Condens. Matter Mater. Phys., 2010, 81, 125304.
2. H. Zhou, T. Xie, T. Taniguchi, K. Watanabe and A. F. Young, Nature, 2021, 598, 434.
3. M. Yankowitz, J. I.-J. Wang, A. G. Birdwell, Y.-A. Chen, K. Watanabe, T. Taniguchi, P. Jacquod, P. San-Jose, P. Jarillo-Herrero and B. J. LeRoy, Nat. Mater., 2014, 13, 786.
4. F. R. Geisenhof, F. Winterer, S. Wakolbinger, T. D. Gokus, Y. C. Durmaz, D. Priesack, J. Lenz, F. Keilmann, K. Watanabe, T. Taniguchi, R. Guerrero-Aviles, M. Pelc, A. Ayuela and R. T. Weitz, ACS Appl. Nano Mater., 2019, 2, 6067–6075.
5. C. Bouhafs, S. Pezzini, F. R. Geisenhof, N. Mishra, V. Miseikis, Y. Niu, C. Struzzi, R. T. Weitz, A. A. Zakharov, S. Forti and C. Coletti, Carbon, 2021, 177, 282–290.
6. Z. Gao, S. Wang, J. Berry, Q. Zhang, J. Gebhardt, W. M. Parkin, J. Avila, H. Yi, C. Chen, S. Hurtado-Parra, M. Drndic, A. M. Rappe, D. J. Srolovitz, J. M. Kikkawa, Z. Luo, M. C. Asensio, F. Wang and A. T. C. Johnson, Nat. Commun., 2020, 11, 546.
7. H. Li, M. I. B. Utama, S. Wang, W. Zhao, S. Zhao, X. Xiao, Y. Jiang, L. Jiang, T. Taniguchi, K. Watanabe, A. Weber-Bargioni, A. Zettl and F. Wang, Nano Lett., 2020, 20, 3106–3112.
8. A. Kerelsky, C. Rubio-Verdu, L. Xian, D. M. Kennes, D. Halbertal, N. Finney, L. Song, S. Turkel, L. Wang, K. Watanabe, T. Taniguchi, J. Hone, C. Dean, D. N. Basov, A. Rubio and A. N. Pasupathy, Proc. Natl. Acad. Sci. U. S. A., 2021, 118, e2017366118.
9. J. P. Nery, M. Calandra and F. Mauri, Nano Lett., 2020, 20, 5017–5023.
10. Q. Lin, T. Li, Z. Liu, Y. Song, L. He, Z. Hu, Q. Guo and H. Ye, Carbon, 2012, 50, 2369–2371.
11. Y. Henni, H. P. O. Collado, K. Nogajewski, M. R. Molas, G. Usaj, C. A. Balseiro, M. Orlita, M. Potemski and C. Faugeras, Nano Lett., 2016, 16, 3710–3716.
12. H. Henck, J. Avila, Z. B. Aziza, D. Pierucci, J. Baima, B. Pamuk, J. Chaste, D. Utt, M. Bartos, K. Nogajewski, B. A. Piot, M. Orlita, M. Potemski, M. Calandra, M. C. Asensio, F. Mauri, C. Faugeras and A. Ouerghi, Phys. Rev. B, 2018, 97, 245421.
13. Y. Yang, Y.-C. Zou, C. R. Woods, Y. Shi, J. Yin, S. Xu, S. Ozdemir, T. Taniguchi, K. Watanabe, A. K. Geim, K. S. Novoselov, S. J. Haigh and A. Mishchenko, Nano Lett., 2019, 19, 8526–8532.
14. Y. Shi, S. Xu, Y. Yang, S. Slizovskiy, S.-K. Son, S. Ozdemir, C. Mullan, J. Barrier, J. Yin, A. I. Berdyugin, B. A. Piot, T. Taniguchi, K. Watanabe, V. I. Fal’ko, K. S. Novoselov, A. K. Geim and A. Mishchenko, Nature, 2020, 584, 210.
15. Z. Liu, S.-M. Zhang, J.-R. Yang, J. Z. Liu, Y.-L. Yang and Q.-S. Zheng, Acta Mech. Sin., 2012, 28, 978–982.
16. Y. Liu, F. Grey and Q. Zheng, Sci. Rep., 2014, 4, 4875.
17. R. C. Sinclair, J. L. Suter and P. V. Coveney, Adv. Mater., 2018, 30, 1705791.
18. M. Z. Baykara, M. R. Vazirisereshk and A. Martini, Appl. Phys. Rev., 2018, 5, 041102.
19. G. Mills and H. Jonsson, Phys. Rev. Lett., 1994, 72, 1124–1127.
20. A. Poma, M. Monteferrante, S. Bonella and G. Ciccotti, Phys. Chem. Chem. Phys., 2012, 14, 15458.
21. E. M. McIntosh, K. T. Wikfeldt, J. Ellis, A. Michaelides and W. Allison, J. Phys. Chem. Lett., 2013, 4, 1565–1569.
22. J. R. Cendagorta, A. Powers, T. J. H. Hele, O. Marsalek, Z. Bacic and M. E. Tuckerman, Phys. Chem. Chem. Phys., 2016, 18, 32169–32177.
23. M. Dienwiebel, G. S. Verhoeven, N. Pradeep, J. W. M. Frenken, J. A. Heimberg and H. W. Zandbergen, Phys. Rev. Lett., 2004, 92, 126101.
24. M. Li, J. Shi, L. Liu, P. Yu, N. Xi and Y. Wang, Sci. Technol. Adv. Mater., 2016, 17, 189–199.
25. X.-Z. Liu, Z. Ye, Y. Dong, P. Egberts, R. Carpick and A. Martini, Phys. Rev. Lett., 2015, 114, 146102.
26. I. V. Lebedeva, A. V. Lebedev, A. M. Popov and A. A. Knizhnik, Comput. Mater. Sci., 2017, 128, 45–58.
27. M. Mohr, J. Maultzsch, E. Dobardzic, S. Reich, I. Milosevic, M. Damnjanovic, A. Bosak, M. Krisch and C. Thomsen, Phys. Rev. B: Condens. Matter Mater. Phys., 2007, 76, 035439.
28. R. Nicklow, N. Wakabayashi and H. G. Smith, Phys. Rev. B, 1972, 5, 4951–4962.
29. 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, Nat. Mater., 2012, 11, 294–300.
30. M. Hanfland, H. Beister and K. Syassen, Phys. Rev. B: Condens. Matter Mater. Phys., 1989, 39, 12598.
31. D. Boschetto, L. Malard, C. H. Lui, K. F. Mak, Z. Li, H. Yan and T. F. Heinz, Nano Lett., 2013, 13, 4620.
32. M. H. S. Amelinckx, P. Delavignette and M. Heerschap, Chem. Phys. Carbon, 1966, 1, 1.
33. C. Cong and T. Yu, Nat. Commun., 2014, 5, 1–7.
34. T. Mishina, K. Nitta and Y. Masumoto, Phys. Rev. B: Condens. Matter Mater. Phys., 2000, 62, 2908.
35. G. A. McQuade, A. S. Plaut, A. Usher and J. Martin, Appl. Phys. Lett., 2021, 118, 203101.
36. R. Mittal, M. K. Gupta, B. Singh, S. Mishra and S. L. Chaplot, Solid State Commun., 2021, 332, 114324.
37. A. C. Bailey and B. Yates, J. Appl. Phys., 1970, 41, 5088.
38. W. Morgan, Carbon, 1972, 10, 73–79.
39. C. P. Herrero and R. Ramirez, Phys. Rev. B, 2020, 101, 035405.
40. B. Marsden, A. Mummery and P. Mummery, Proc. R. Soc. A, 2018, 474, 20180075.
41. I. Errea, M. Calandra and F. Mauri, Phys. Rev. B: Condens. Matter Mater. Phys., 2014, 89, 064302.
42. L. Monacelli, R. Bianco, M. Cherubini, M. Calandra, I. Errea and F. Mauri, J. Phys.: Condens. Matter, 2021, 33, 363001.
43. P. B. Allen, Mod. Phys. Lett. B, 2019, 34, 2050025.
44. F. Libbi, N. Bonini and N. Marzari, 2D Mater., 2020, 8, 015026.
45. M. Wen, S. Carr, S. Fang, E. Kaxiras and E. B. Tadmor, Phys. Rev. B, 2018, 98, 235404.
46. S. Zhou, J. Han, S. Dai, J. Sun and D. J. Srolovitz, Phys. Rev. B: Condens. Matter Mater. Phys., 2015, 92, 155438.
47. J. Paul Nery, M. Calandra and F. Mauri, 2D Mater., 2021, 8, 035006.
48. N. Blomquist, M. Alimadadi, M. Hummelgard, C. Dahlstrom, M. Olsen and H. Olin, Sci. Rep., 2019, 9, 8966.
49. L. Jiang, S. Wang, Z. Shi, C. Jin, M. I. B. Utama, S. Zhao, Y.-R. Shen, H.-J. Gao, G. Zhang and F. Wang, Nat. Nanotechnol., 2018, 13, 204–208.
50. R. He, J. van Baren, J.-A. Yan, X. Xi, Z. Ye, G. Ye, I.-H. Lu, S. M. Leong and C. H. Lui, 2D Mater., 2016, 3, 031008.
51. R. Cusco, J. H. Edgar, S. Liu, G. Cassabois, B. Gil and L. Artus, Phys. Rev. B, 2019, 99, 085428.
52. R. Ramirez, E. Chacon and C. P. Herrero, Phys. Rev. B, 2016, 93, 235419.
53. L. Karssemeijer and A. Fasolino, Surf. Sci., 2011, 605, 1611–1615.
54. N. Mounet and N. Marzari, Phys. Rev. B: Condens. Matter Mater. Phys., 2005, 71, 205214.

---