# Atmospheric Retrievals with petitRADTRANS

Evert Nasedkin  Paul Mollière  Doriann Blain

††margin:  DOI:
https://doi.org/10.5281/zenodo.8337871https://doi.org10.5281/zenodo.8337871
Software • https://doi.orgReview •
https://gitlab.com/mauricemolli/petitRADTRANShttps://doi.orgRepository •
http://dx.doi.org/10.5281/zenodo.8337871https://doi.orgArchive Submitted: 12
September 2023
Published: Unpublished License
Authors of papers retain copyright and release the work under a Creative
Commons Attribution 4.0 International License
(https://creativecommons.org/licenses/by/4.0/https://doi.orgCC BY 4.0).

## Summary

petitRADTRANS (pRT) is a fast radiative transfer code used for computing
emission and transmission spectra of exoplanet atmospheres (Mollière et al.
2019), combining a FORTRAN back end with a Python based user interface. It is
widely used in the exoplanet community with 161 references in the literature
to date, and has been benchmarked against numerous similar tools, including
many listed in (MacDonald and Batalha 2023). The spectra calculated with pRT
can be used as a forward model for fitting spectroscopic data using Monte
Carlo techniques, commonly referred to as an atmospheric retrieval
(Madhusudhan and Seager 2009). The new retrieval module combines fast forward
modelling with nested sampling codes, allowing for atmospheric retrievals on a
large range of different types of exoplanet data. Thus it is now possible to
use pRT to easily and quickly infer the atmospheric properties of exoplanets
in both transmission and thermal emission.

## Statement of need

Atmospheric retrievals are a cornerstone of exoplanet atmospheric
characterisation. pRT provides a powerful and user-friendly tool for
researchers to fit exoplanet spectra with a range of built in or custom
atmospheric models. Various thermal structures, chemistry and cloud
parameterisations and opacity calculation methods can be combined and used to
perform parameter estimation and model comparison for a given atmospheric
spectrum. With increasing volumes of both ground- and space-based spectra
available, it is necessary for exoplanet researchers to have access to a range
of characterisation tools.

## petitRADTANS Retrieval Module

The Retrieval module combines the Radtrans forward modelling class with a
nested sampler via a likelihood function to perform an atmospheric retrieval.
Both MultiNest (F. Feroz and Hobson 2008; F. Feroz, Hobson, and Bridges 2009;
Farhan Feroz et al. 2019; J. Buchner et al. 2014) and Ultranest (J. Buchner et
al. 2014; Johannes Buchner 2019) samplers are available, with both offering
MPI implementations that allow for easy parallelisation.

Datasets, priors and other retrieval hyper parameters are set through the
RetrievalConfig class, while the models module includes a range of complete
atmospheric models that can be fit to the data. Users can also define their
own model function, either by making use of temperature profiles from the
physics module and chemistry parameterisations from the chemistry module or by
implementing their own forward model.

Figure 1: Typical example of default pRT outputs. This highlights the fit of a
transmission spectrum model to JWST/NIRISS/SOSS data of WASP 39 b as part of
the Transiting Early Release Science program.

Multiple datasets can be included into a single retrieval, with each dataset
receiving its own Radtrans object used for the radiative transfer calculation
where some or all forward model parameters may be shared between the different
data sets. This allows for highly flexible retrievals where multiple spectral
resolutions, wavelength ranges and even atmospheric models can be combined in
a single retrieval. Each dataset can also receive scaling factors (for the
flux, uncertainties or both), error inflation factors and offsets. The model
functions are used to compute a spectrum $\vec{S}$, which is convolved to the
instrumental resolution and binned to the wavelength bins of the data using a
custom binning function to account for non-uniform bin sizes. The resulting
spectrum compared to the data with flux $\vec{F}$ and covariance $\mathbf{C}$
in the likelihood function:

$-2\log\mathcal{L}=\left(\vec{S}-\vec{F}\right)^{T}\mathbf{C}^{-1}\left(\vec{S}-\vec{F}\right)+\log\left(2\pi\det\left(\mathbf{C}\right)\right).$
(1)

The second term is included in the likelihood to allow for uncertainties to
vary as a free parameter during the retrieval, and penalizes overly large
uncertainties. An likelihood function for high resolution data based on that
of (Gibson et al. 2020) is also available.

pRT can compute spectra either using line-by-line calculations, or using
correlated-k (c-k) tables for defining the opacities of molecular species. We
include up-to-date correlated-k line lists from Exomol (Tennyson and Yurchenko
2012; McKemmish, Yurchenko, and Tennyson 2016; Polyansky et al. 2018; Chubb et
al. 2020) and HITEMP Hargreaves et al. (2020), with the full set of available
opacities listed in the online documentation. The `exo-k` package is used to
resample the the correlated-k opacity tables to a lower spectral resolution in
order to reduce the computation time (Leconte 2021). Combining the c-k
opacities of multiple species requires mixing the distributions in $g$ space.
This operation is necessary when calculating emission spectra and accounting
for multiple scattering in the clouds. Previously, this was accomplished by
taking 1000 samples of each distribution. This sampling process resulted in
non-deterministic spectral calculations with a small (up to 1%) scatter about
the expected result, which could lead to unexpected behaviour from the nested
sampler as the same set of parameters could result in different log-
likelihood. pRT has been updated to fully mix the c-k distributions,
iteratively mixing in any species with a significant opacity contribution: a
species is only mixed in if the highest opacity value in a given spectral bin
is larger than a threshold value. This theshold value is obtained by listing
the smallest opacity value for every species in a given spectral bin, and then
setting the threshold to 1% of the largest value from the list for each
spectral bin. The resulting grid is linearly interpolated back to the 16 $g$
points at each pressure and frequency bin as required by pRT. This fully
deterministic process produces stable log-likelihood calculations, and
resulted in a 5$\times$ improvement in the speed of the c-k mixing function.

Various thermal, chemical and cloud parameterisations are available in pRT.
Built in temperature profiles range from interpolated splines to physically
motivated profiles as in Guillot (2010) and Mollière et al. (2020).
Equilibrium and disequilibrium chemistry can be interpolated from a pre-
computed grid on-the-fly. Chemical abundances can also be freely retrieved,
with the additional possibility of using a combination of free and chemically
consistent abundances. Cloud parametersiations range from a ‘grey’ continuum
opacity applied at all wavelengths, to clouds parameterised as in Ackerman and
Marley (2001), using log-normal or Hansen (1971) particle size distributions
with real optical opacities for different compositions and particle shapes,
and including self-scattering. Users can also pass functions to the Radtrans
object that encode any absorption or scattering opacity as a function of
wavelength and atmospheric position, allowing for a generic cloud
implementation. Included in pRT is an option to use an adaptive pressure grid
with a higher resolution around the location of the cloud base, allowing for
more precise positioning of the cloud layers within the atmosphere.

Photometric data are fully incorporated into the retrieval process. The
spectral model is multiplied by a filter transmission profile from the SVO
database using the species package (Stolker et al. 2020). This results in
accurate synthetic photometry, which can be compared to the values specied by
the user with the add_photometry function.

Publication-ready summary plots of best fits, temperature and abundance
profiles and corner plots can be automatically generated. Multiple retrieval
results can be combined in the plots for model intercomparisons. Such results
have been benchmarked against other widely used retrieval codes, in particular
as part of the JWST Early Release Science (ERS) program (Welbanks et al, in
prep). Figure 1 shows the fit of a transmission model to the JWST/NIRISS/SOSS
observations of WASP 39 b (Feinstein et al. 2023) from the ERS program.

## Acknowledgements

We gratefully acknowledge contributions to petitRADTRANS from Eleonara Alei,
Karan Molaverdikhani, Francois Rozet, Aaron David Schneider, Tomas Stolker,
Nick Wogan and Mantas Zilinskas.

## References

Ackerman, Andrew S., and Mark S. Marley. 2001. “Precipitating Condensation
Clouds in Substellar Atmospheres” 556 (2): 872–84.
https://doi.org/10.1086/321540.

Buchner, J., A. Georgakakis, K. Nandra, L. Hsu, C. Rangel, M. Brightman, A.
Merloni, M. Salvato, J. Donley, and D. Kocevski. 2014. “X-Ray Spectral
Modelling of the AGN Obscuring Region in the CDFS: Bayesian Model Selection
and Catalogue” 564 (April): A125. https://doi.org/10.1051/0004-6361/201322971.

Buchner, Johannes. 2019. “Collaborative Nested Sampling: Big Data versus
Complex Physical Models” 131 (1004): 108005.
https://doi.org/10.1088/1538-3873/aae7fc.

Chubb, K. L., M. Rocchetto, A. F. Al-Refaie, I. Waldmann, M. Min, J. Barstow,
P. Molliére, M. W. Phillips, J. Tennyson, and S. N. Yurchenko. 2020\. “The
ExoMolOP Database: Cross-Sections and k-Tables for Molecules of Interest in
High-Temperature Exoplanet Atmospheres.” _A &A_.

Feinstein, Adina D., Michael Radica, Luis Welbanks, Catriona Anne Murray,
Kazumasa Ohno, Louis-Philippe Coulombe, Néstor Espinoza, et al. 2023\. “Early
Release Science of the exoplanet WASP-39b with JWST NIRISS” 614 (7949):
670–75. https://doi.org/10.1038/s41586-022-05674-1.

Feroz, Farhan, Michael P. Hobson, Ewan Cameron, and Anthony N. Pettitt. 2019\.
“Importance Nested Sampling and the MultiNest Algorithm.” _The Open Journal of
Astrophysics_ 2 (1): 10. https://doi.org/10.21105/astro.1306.2144.

Feroz, F., and M. P. Hobson. 2008. “Multimodal Nested Sampling: An Efficient
and Robust Alternative to Markov Chain Monte Carlo Methods for Astronomical
Data Analyses” 384 (2): 449–63.
https://doi.org/10.1111/j.1365-2966.2007.12353.x.

Feroz, F., M. P. Hobson, and M. Bridges. 2009. “MULTINEST: An Efficient and
Robust Bayesian Inference Tool for Cosmology and Particle Physics” 398 (4):
1601–14. https://doi.org/10.1111/j.1365-2966.2009.14548.x.

Gibson, Neale P., Stephanie Merritt, Stevanus K. Nugroho, Patricio E.
Cubillos, Ernst J. W. de Mooij, Thomas Mikal-Evans, Luca Fossati, et al.
2020\. “Detection of Fe I in the atmosphere of the ultra-hot Jupiter
WASP-121b, and a new likelihood-based approach for Doppler-resolved
spectroscopy” 493 (2): 2215–28. https://doi.org/10.1093/mnras/staa228.

Guillot, T. 2010. “On the radiative equilibrium of irradiated planetary
atmospheres” 520 (September): A27.
https://doi.org/10.1051/0004-6361/200913396.

Hansen, James E. 1971. “Multiple Scattering of Polarized Light in Planetary
Atmospheres Part II. Sunlight Reflected by Terrestrial Water Clouds.” _Journal
of Atmospheric Sciences_ 28 (8): 1400–1426.
https://doi.org/10.1175/1520-0469(1971)028%3C1400:MSOPLI%3E2.0.CO;2.

Hargreaves, Robert J., Iouli E. Gordon, Michael Rey, Andrei V. Nikitin,
Vladimir G. Tyuterev, Roman V. Kochanov, and Laurence S. Rothman. 2020. “An
Accurate, Extensive, and Practical Line List of Methane for the HITEMP
Database” 247 (2): 55. https://doi.org/10.3847/1538-4365/ab7a1a.

Leconte, Jérémy. 2021. “Spectral binning of precomputed correlated-k
coefficients” 645 (January): A20. https://doi.org/10.1051/0004-6361/202039040.

MacDonald, Ryan J., and Natasha E. Batalha. 2023. “A Catalog of Exoplanet
Atmospheric Retrieval Codes.” _Research Notes of the American Astronomical
Society_ 7 (3): 54. https://doi.org/10.3847/2515-5172/acc46a.

Madhusudhan, N., and S. Seager. 2009. “A Temperature and Abundance Retrieval
Method for Exoplanet Atmospheres” 707 (1): 24–39.
https://doi.org/10.1088/0004-637X/707/1/24.

McKemmish, L. K., S. N. Yurchenko, and J. Tennyson. 2016. “ExoMol Molecular
Linelists – XVIII. The Spectrum of Vanadium Oxide.” _Mon. Not. R. Astron.
Soc._ 463: 771–93. https://doi.org/10.1093/mnras/stw1969.

Mollière, P., T. Stolker, S. Lacour, G. P. P. L. Otten, J. Shangguan, B.
Charnay, T. Molyarova, et al. 2020. “Retrieving scattering clouds and
disequilibrium chemistry in the atmosphere of HR 8799e” 640 (August): A131.
https://doi.org/10.1051/0004-6361/202038325.

Mollière, P., J. P. Wardenier, R. van Boekel, Th. Henning, K. Molaverdikhani,
and I. A. G. Snellen. 2019. “petitRADTRANS: A Python Radiative Transfer
Package for Exoplanet Characterization and Retrieval.” _A &A_ 627 (July): A67.
https://doi.org/10.1051/0004-6361/201935470.

Polyansky, O. L., A. A. Kyuberis, N. F. Zobov, J. Tennyson, S. N. Yurchenko,
and L. Lodi. 2018. “ExoMol Molecular Line Lists XXX: A Complete High-Accuracy
Line List for Water.” _Mon. Not. R. Astron. Soc._ 480: 2597–2608.
https://doi.org/10.1093/mnras/sty1877.

Rothman, L. S., I. E. Gordon, R. J. Barber, H. Dothe, R. R. Gamache, A.
Goldman, V. I. Perevalov, S. A. Tashkun, and J. Tennyson. 2010. “HITEMP, the
High-Temperature Molecular Spectroscopic Database” 111 (October): 2139–50.
https://doi.org/10.1016/j.jqsrt.2010.05.001.

Stolker, T., S. P. Quanz, K. O. Todorov, J. Kühn, P. Mollière, M. R. Meyer, T.
Currie, S. Daemgen, and B. Lavie. 2020. “MIRACLES: atmospheric
characterization of directly imaged planets and substellar companions at 4-5
$\mu$m. I. Photometric analysis of $\beta$ Pic b, HIP 65426 b, PZ Tel B, and
HD 206893 B” 635 (March): A182. https://doi.org/10.1051/0004-6361/201937159.

Tennyson, Jonathan, and Sergei N. Yurchenko. 2012. “ExoMol: Molecular Line
Lists for Exoplanet and Other Atmospheres” 425 (1): 21–33.
https://doi.org/10.1111/j.1365-2966.2012.21440.x.


# Interlayer interaction, shear vibrational mode, and tribological properties
of two-dimensional bilayers with a commensurate moiré pattern

Alexander S. Minkin<EMAIL_ADDRESS>Keldysh Institute of Applied Mathematics of
Russian Academy of Sciences, 4 Miusskaya sq., Moscow, 125047, Russia  Irina V.
Lebedeva<EMAIL_ADDRESS>Simune Atomistics, Avenida de Tolosa 76, San
Sebastian 20018, Spain  Andrey M. Popov<EMAIL_ADDRESS>Institute for
Spectroscopy of Russian Academy of Sciences, Fizicheskaya str. 5, Troitsk,
Moscow 108840, Russia  Sergey A. Vyrko<EMAIL_ADDRESS>Physics Department,
Belarusian State University, Nezavisimosti Ave. 4, Minsk 220030, Belarus
Nikolai A. Poklonski  Corresponding author<EMAIL_ADDRESS>Physics
Department, Belarusian State University, Nezavisimosti Ave. 4, Minsk 220030,
Belarus  Yurii E. Lozovik<EMAIL_ADDRESS>Institute for Spectroscopy
of Russian Academy of Sciences, Fizicheskaya str. 5, Troitsk, Moscow 108840,
Russia  Moscow Institute of Electronics and Mathematics, National Research
University Higher School of Economics, Bol. Trekhsvjatitel’skij per., 1-3/12,
build. 8, Moscow, 101000, Russia

###### Abstract

The potential energy surface (PES) of interlayer interaction of infinite
twisted bilayer graphene is calculated for a set of commensurate moiré
patterns using the registry-dependent Kolmogorov–Crespi empirical potential.
The calculated PESs have the same shape for all considered moiré patterns with
the unit cell size of the PES which is inversely related to the unit cell size
of the moiré pattern. The amplitude of PES corrugations is found to decrease
exponentially upon increasing the size of the moiré pattern unit cell. An
analytical expression for such a PES including the first Fourier harmonics
compatible with the symmetries of both layers is derived. It is shown that the
calculated PESs can be approximated by the derived expression with the
accuracy within 1%. This means that different physical properties associated
with relative in-plane motion of graphene layers are interrelated and can be
expressed analytically as functions of the amplitude of PES corrugations. In
this way, we obtain the shear mode frequency, shear modulus, shear strength
and barrier for relative rotation of the commensurate twisted layers to a
fully incommensurate state for the considered moiré patterns. This barrier may
possibly lead to the macroscopic robust superlubricity for twisted graphene
bilayer with a commensurate moiré pattern. The conclusions made should be
valid for diverse 2D systems of twisted commensurate layers.

##  I Introduction

Structural superlubricity, i.e. the mode of relative motion of the layers with
vanishing or nearly vanishing friction Hirano and Shinjo (1990); Hirano  _et
al._ (1991), has attracted a considerable attention in the context of
discovery of graphene and other 2D materials, see Ref. Hod  _et al._ , 2018
for a review. First this phenomenon was observed for nanoscale contacts
between graphene flakes at the tip of a microscope probe and graphite surface
Verhoeven  _et al._ (2004); Dienwiebel  _et al._ (2005); Filippov  _et al._
(2008). A wide set of atomistic simulations has been devoted to superlubricity
for 2D systems with a finite size of the contact area where the edge or rim
contribution to the static friction is dominant Verhoeven  _et al._ (2004);
Filippov  _et al._ (2008); Xu  _et al._ (2013); Koren and Duerig (2016);
Bonelli  _et al._ (2009); van Wijk  _et al._ (2013); Guo  _et al._ (2007);
Shibuta and Elliott (2011); Wang  _et al._ (2019); Zhang  _et al._ (2015);
Zhang and Chang (2018); Zhang  _et al._ (2022). Recently not only nanoscale
but also micro- and macroscale superlubricity has been found in systems of 2D
layers Androulidakis  _et al._ (2020); Liu  _et al._ (2012); Vu  _et al._
(2016); Song  _et al._ (2018). These studies raise interest to possible
factors which cause the static friction and can restrict superlubricity for a
macroscale incommensurate contact area Liu  _et al._ (2012); Mandelli  _et
al._ (2017); Hod  _et al._ (2018); Koren and Duerig (2016); Minkin  _et al._
(2021, 2022). The following possible reasons of very low but nevertheless
nonzero static friction have been considered: 1) contribution of incomplete
unit cells located at the rim area of one of the layers forming a moiré
pattern (rim contribution), Koren and Duerig (2016) 2) incomplete static
friction force cancellation within complete unit cells of a commensurate moiré
pattern (area contribution) Koren and Duerig (2016), 3) motion of domain walls
of large commensurate domains formed upon relaxation of moiré patterns Liu
_et al._ (2012); Mandelli  _et al._ (2017); Hod  _et al._ (2018), and 4)
contribution of atomic-scale defects Liu  _et al._ (2012); Minkin  _et al._
(2021, 2022). The present paper is devoted to the detailed study of the area
contribution to the static friction by the example of twisted graphene
bilayer.

Whereas for very small twist angles, the size of the moiré pattern unit cell
is large Campanera  _et al._ (2007) and, therefore, formation of commensurate
domains separated by incommensurate domain boundaries occurs during the
structural relaxation Lebedeva and Popov (2020), for twist angles far from the
coaligned orientation (0∘, 60∘ and so on), the size of the moiré pattern unit
cell is smaller or comparable with the width of commensurate domain walls
(about 10 nm for bilayer graphene Popov  _et al._ (2011)). For such angles,
relaxation to the commensurate domains is not possible and a set of
commensurate moiré patterns can be observed Mele (2012); Campanera  _et al._
(2007). Twisted bilayer graphene with a commensurate moiré pattern has the
interlayer interaction energy slightly lower than in a fully incommensurate
state Xu  _et al._ (2013). Thus, one can expect that than such patterns can be
formed preferably for the corresponding range of twist angles. Here we propose
that the energetic preference of commensurate moiré patterns can also lead to
the robust superlubricity. This is why the study of tribological properties of
2D systems with commensurate moiré patterns is of high interest.

The tribological properties are determined by the potential energy surface
(PES) of interlayer interaction that is the dependence of this energy on
coordinates describing the relative in-plane displacement of 2D layers.
Previous atomistic calculations allowed to distinguish the area contribution
into the static friction of twisted graphene bilayer for a few commensurate
moiré patterns Xu  _et al._ (2013); Koren and Duerig (2016); Kabengele and
Johnson (2021). Cancellation of the static friction force within complete unit
cells of commensurate moiré patterns of graphene bilayer Koren and Duerig
(2016); Minkin  _et al._ (2021, 2022) and double-walled carbon nanotubes
Kolmogorov and Crespi (2000); Belikov  _et al._ (2004); Bichoutskaia  _et al._
(2006) has been also demonstrated. However, the symmetry and shape of the PES
of interlayer interaction for commensurate moiré patterns of infinite twisted
graphene bilayers has not been studied yet. Recently we proposed a hypothesis
that such PESs in diverse 2D materials with layers aligned in the same or
opposite directions can be universally described by the first spatial Fourier
harmonics Lebedev  _et al._ (2020). This hypothesis has been confirmed by
calculations of PESs for different 2D materials Ershova  _et al._ (2010);
Lebedeva  _et al._ (2011a); Popov  _et al._ (2012a); Lebedeva  _et al._
(2012); Zhou  _et al._ (2015); Reguzzoni  _et al._ (2012); Lebedev  _et al._
(2016); Zhou  _et al._ (2015); Lebedev  _et al._ (2020) and 2D
heterostructures Popov  _et al._ (2013); Jung  _et al._ (2015); Kumar  _et
al._ (2015); Lebedev  _et al._ (2017). Moreover, this hypothesis is valid also
for double-walled carbon nanotubes Vuković  _et al._ (2003); Belikov  _et al._
(2004); Bichoutskaia  _et al._ (2005, 2009a); Popov  _et al._ (2009, 2012b),
where only Fourier harmonics compatible with the symmetry of both wall
contribute into the PES of interwall interaction Vuković  _et al._ (2003);
Popov  _et al._ (2009). By analogy with double-walled nanotubes, one might
expect that PES for interlayer interaction of twisted layers of infinite
commensurate moiré pattern is determined by Fourier harmonics compatible with
symmetries of the both layers, i.e. with the symmetry of the whole moiré
pattern. In such a case, the approximated PES depends on a single parameter
and a set of physical quantities determined by the PES are interrelated Popov
_et al._ (2012a); Lebedev  _et al._ (2016).

Here we calculate PES for a wide set of infinite commensurate moiré patterns
of twisted graphene bilayers using the registry-dependent Kolmogorov–Crespi
potential Kolmogorov and Crespi (2005) and show that these PESs can be
excellently approximated by the first Fourier harmonics which are compatible
with the symmetry of the whole moiré pattern. The PES approximation derived is
used to obtain analytical expressions for a set of physical quantities of
moiré patterns determined by the PES such as the shear frequency, shear
modulus, shear strength and barrier for relative rotation of the layers to a
fully incommensurate state.

The paper is organized in the following way. In Sec. II, the model of the
superlubric system and calculation methods are described. Sec. III is devoted
to our results on the PES calculation and approximation by Fourier harmonics
as well as estimates of physical quantities determined by the PES. The
conclusions and discussion are presented in Sec. IV.

##  II Methodology

Figure 1: (Color online) (a) and (b) schemes of twin commensurate moiré
patterns (2,1) and (4,1) of twisted graphene bilayer, respectively, with the
same size of the unit cell. Lattice vectors $\bm{\mathrm{a}}_{1}$ and
$\bm{\mathrm{a}}_{2}$ of the bottom graphene layer and
$\bm{\mathrm{a}}_{1}^{\prime}$ and $\bm{\mathrm{a}}_{2}^{\prime}$ of the top
layer, lattice vectors $\bm{\mathrm{L}}_{1}$ and $\bm{\mathrm{L}}_{2}$ of the
commensurate moiré pattern, angle $\theta$ of relative rotation of the
graphene layers and angle $\varphi$ between the lattice vector
$\bm{\mathrm{a}}_{1}$ of the bottom layer and lattice vector of commensurate
moiré pattern $\bm{\mathrm{L}}_{1}$ are indicated. The translational
displacement $\bm{\mathrm{t}}$ of the upper layer that converts one moiré
pattern into the other are shown by black arrows.

###  II.1 Structure of commensurate twisted graphene bilayer

Let us consider first the structure of commensurate twisted graphene bilayer.
The commensurate moiré pattern $(n_{1},n_{2})$ is defined by the indices
$n_{1}$ and $n_{2}$ which are coprime numbers Mele (2012). It has been shown
that for each moiré pattern $(n_{1},n_{2})$, there is a twin pattern
$(n_{1}^{\prime},n_{2}^{\prime})$ with greater indices $n_{1}^{\prime}$ and
$n_{2}^{\prime}$ with the same size of the unit cell Mele (2012); Campanera
_et al._ (2007). In the pairs of indices of twin moiré patterns,
$(n_{1}-n_{2})/3$ is not integer for the smaller pair of indices, whereas
$(n_{1}^{\prime}-n_{2}^{\prime})/3$ is integer for the greater indices. The
examples of twin commensurate moiré patterns are shown in Fig. 1.

If $(n_{1}-n_{2})/3$ is not integer, the unit cell of the commensurate moiré
pattern $(n_{1},n_{2})$ is defined by lattice vectors $\bm{\mathrm{L}}_{1}$
and $\bm{\mathrm{L}}_{2}$ [see Fig. 1(a)]

$\bm{\mathrm{L}}_{1}=n_{1}\bm{\mathrm{a}}_{1}+n_{2}\bm{\mathrm{a}}_{2},\quad\bm{\mathrm{L}}_{2}=-n_{2}\bm{\mathrm{a}}_{1}+(n_{1}+n_{2})\bm{\mathrm{a}}_{2},$

where $\bm{\mathrm{a}}_{1}$ and $\bm{\mathrm{a}}_{2}$ are lattice vectors of
the bottom graphene layer.

If $(n_{1}^{\prime}-n_{2}^{\prime})/3$ is integer, the same equations give the
vectors one of which corresponds to the diagonal of the moiré pattern unit
cell [see Fig. 1(b)]

$\bm{\mathrm{L}}_{1}^{\prime}=n_{1}^{\prime}\bm{\mathrm{a}}_{1}+n_{2}^{\prime}\bm{\mathrm{a}}_{2},\quad\bm{\mathrm{L}}_{2}^{\prime}=-n_{2}^{\prime}\bm{\mathrm{a}}_{1}+(n_{1}^{\prime}+n_{2}^{\prime})\bm{\mathrm{a}}_{2}.$

In this case the indices $n_{1}$ and $n_{2}$ which determine the moiré pattern
unit cell

$\bm{\mathrm{L}}_{1}^{\text{(u)}}=n_{1}\bm{\mathrm{a}}_{2}+n_{2}\bm{\mathrm{a}}_{1},\quad\bm{\mathrm{L}}_{2}^{\text{(u)}}=-n_{2}\bm{\mathrm{a}}_{2}+(n_{1}+n_{2})\bm{\mathrm{a}}_{1}$

for given $n_{1}^{\prime}$ and $n_{2}^{\prime}$ can be found from any of the
equations

$\bm{\mathrm{L}}_{1}^{\prime}=\bm{\mathrm{L}}_{1}^{\text{(u)}}+\bm{\mathrm{L}}_{2}^{\text{(u)}},\quad\bm{\mathrm{L}}_{2}^{\prime}=2\bm{\mathrm{L}}_{1}^{\text{(u)}}-\bm{\mathrm{L}}_{2}^{\text{(u)}}$

in the form

$n_{1}=\frac{{n_{1}}^{\prime}+2{n_{2}}^{\prime}}{3},\quad
n_{2}=\frac{n_{1}^{\prime}-n_{2}^{\prime}}{3}.$

The angle $\theta$ of relative rotation of graphene layers of the commensurate
moiré pattern (that is the angle between the vectors $\bm{\mathrm{a}}_{1}$ and
$\bm{\mathrm{a}}_{1}^{\prime}$) is defined as

$\cos\theta=\frac{n_{1}^{2}+4n_{1}n_{2}+n_{2}^{2}}{2(n_{1}^{2}+n_{1}n_{2}+n_{2}^{2})}.$

In the case where $(n_{1}-n_{2})/3$ is not integer, the angle $\varphi$
between the lattice vector $\bm{\mathrm{a}}_{1}$ and the lattice vector of
commensurate moiré pattern
$\bm{\mathrm{L}}_{1}=n_{1}\bm{\mathrm{a}}_{1}+n_{2}\bm{\mathrm{a}}_{2}$ is

$\varphi=30^{\circ}-\frac{\theta}{2}.$

The area of the unit cell of moiré pattern $(n_{1},n_{2})$ is

$S=S_{g}N_{c}=\frac{\sqrt{3}a^{2}(n_{1}^{2}+n_{1}n_{2}+n_{2}^{2})}{2R},$

where $S_{g}=\sqrt{3}a^{2}/2$ is the area of the unit cell of graphene,
$a=|\bm{\mathrm{a}}_{1}|=|\bm{\mathrm{a}}_{2}|$ is the graphene lattice
constant, $N_{c}=(n_{1}^{2}+n_{1}n_{2}+n_{2}^{2})/R$ is the number of unit
cells of graphene per unit cell of the commensurate moiré pattern, the
parameter $R=3$ if $(n_{1}-n_{2})/3$ is integer and $R=1$ otherwise.

Pairs of twin commensurate moiré patterns with the same size of the unit cell
and with different symmetry of the stacking just after the relative rotation
of layers were considered as different in the original work Mele (2012).
However, these commensurate moiré pattern can be obtained one from another by
the translational displacement $\bm{\mathrm{t}}$ of one of the layers in the
layer plane (see Fig. 1). Here we study the PES of the interlayer interaction
energy as a function of the coordinates describing the in-plane relative
displacement of the layers. Evidently, the moiré patterns related by the
translational displacement $\bm{\mathrm{t}}$ correspond to the same PES. Thus,
only one of pair moiré patterns with the indices $n_{1}$ and $n_{2}$, where
$(n_{1}-n_{2})/3$ is not integer, are considered here for the PES
calculations.

###  II.2 Computational details

Table 1: Calculation details for the considered commensurate moiré patterns with coprime indices $(n_{1},n_{2})$: the angle $\theta$ of relative rotation of graphene layers, the number of atoms $N_{c}$ in the moiré pattern unit cell, simulation cell size in the units of moiré pattern unit cells, total number of atoms $N_{a}$ in the simulation cell and cutoff radius $R_{c}$ of the Kolmogorov–Crespi potential. $(n_{1},n_{2})$ |  $\theta$ (∘) | $N_{c}$ | cell size | $N_{a}$ |  $R_{c}$ (Å)
---|---|---|---|---|---
(2,1) | 21.787 | 7 | $18\times 18$ | 9072 | 16
(3,1) | 32.204 | 13 | $18\times 18$ | 16848 | 16
(3,2) | 13.174 | 19 | $10\times 10$ | 7600 | 16
(5,1) | 42.103 | 31 | $11\times 11$ | 15004 | 70
(5,3) | 16.426 | 49 | $9\times 9$ | 15876 | 70
(7,2) | 35.567 | 67 | $6\times 6$ | 9648 | 50
(7,3) | 26.008 | 79 | $6\times 6$ | 11376 | 50
(7,5) | 10.993 | 109 | $5\times 5$ | 10900 | 50

The ratio of the PES corrugations to the average interlayer interaction energy
is extremely small for twisted graphene bilayers Xu  _et al._ (2013); Koren
and Duerig (2016); Minkin  _et al._ (2022). Simultaneously the size of the
simulation cell and the number of computational runs for each considered moiré
pattern are too high to study the PES of interlayer interaction by ab initio
methods. Thus, we use classical potentials in the present study. At this
moment, there are no experimental data on physical properties of systems of
twisted layers that can be used to fit parameters of classical potentials for
description of the interlayer interaction (or check the adequacy of existing
potentials for twisted layers). The parameters of the popular
Kolmogorov–Crespi and Lebedeva potentials (Tables S1 and S2 in Supplemental
Material) for the interaction between graphene layers were fitted to the PES
of interlayer interaction of coaligned layers (with zero twist angle) obtained
by density functional theory (DFT) calculations Kolmogorov and Crespi (2005);
Lebedeva  _et al._ (2011a). In the case of the Lebedeva potential, the
experimental data on the frequency of in-plane interlayer vibrations of
coaligned layers were also taken into account Popov  _et al._ (2012a).
Nevertheless, corrugations of the PES computed for infinite commesurate moiré
patterns using the Lebedeva potential do not exceed the calculation accuracy
Minkin  _et al._ (2021, 2022). At the same time, they are finite and well-
defined for commensurate moiré patterns with the smallest unit cells when the
Kolmogorov–Crespi potential is usedXu  _et al._ (2013); Koren and Duerig
(2016); Minkin  _et al._ (2022). The explanation of this discrepancy between
the Kolmogorov–Crespi and Lebedeva potentials is discussed in Section IIIB. In
the present study all calculations are performed using the Kolmogorov–Crespi
potential which allows to analyze the shape of the PES determined by the
symmetry of commensurate twisted graphene bilayer. However, we emphasize that
the results obtained here are only of qualitative nature.

Figure 2: (Color online) (a-e) Potential energy $U$ (per atom of the upper
layer) of interlayer interaction of twisted graphene bilayer as a function of
the relative displacement of the layers in the zigzag ($x$, in Å) and armchair
($y$, in Å) directions of the lower layer computed at the optimal interlayer
distance of 3.46 Å for commensurate moiré patterns with coprime indices (a)
(2,1), (b) (3,1), (c) (3,2), (d) (5,1) and (e) (5,3). The energy is given
relative to the minimum. (f) Amplitude of PES corrugations, $\Delta U$ (per
atom of the upper layer), as a function of the number $N_{c}$ of unit cells of
the PES per unit cell of graphene. The indices of the considered moiré
patterns are indicated.

The PES calculations have been carried out under the periodic boundary
conditions. The simulation cells of height 100 Å have been used for all the
considered moiré patterns. The bond length between atoms in the graphene
layers is taken equal 1.42 Å. The upper graphene layer is placed at the
interlayer distance 3.46 Å (which is determined here on the example of (2,1)
moiré pattern to be optimal for the Kolmogorov–Crespi potential) and is
rigidly shifted with respect to the bottom layer with steps of 0.0168 Å and
0.0193 Å in the zigzag and armchair directions of the bottom layer,
respectively. Further calculation details which are different for 8 considered
moiré patterns with the smallest sizes of unit cell are listed in Table 1.

##  III Results

###  III.1 PES of twisted graphene bilayer

The amplitude $\Delta U$ of PES corrugation (i.e. the difference between
maximum and minimum values of the interlayer interaction energy, $\Delta
U=U_{\mathrm{max}}-U_{\mathrm{min}}$) exceeds the calculation accuracy only
for 5 out of 8 considered moiré patterns with smaller sizes of the unit cells.
Calculated PESs for these moiré patterns are shown in Fig. 2. For 3 out of 8
considered moiré patterns with larger sizes of the unit cells (equivalent to
the smaller sizes of the unit cells of the PES as shown below), the amplitude
of PES corrugations is lower than the artifacts related with the finite value
of the cutoff radius of the potential. Calculated PESs for these moiré
patterns are shown in Fig. S1 in Supplemental Material.

Two types of the PESs have been found. The PESs of the first type have a
triangular lattice of minima and honeycomb lattice of maxima whereas the PESs
of the second type, on the contrary, have a triangular lattice of maxima and
honeycomb lattice of minima. As discussed in Section IIIB, the PES shape is
described for both PES types by the same expression which contains only the
first spatial Fourier harmonics and, therefore, only a single energetic
parameter. The difference between two types of the PES is determined by the
sign of this parameter.

It should be noted that the number of the unit cells of the PES per an unit
cell of graphene is the same as the number $N_{c}$ of unit cells of graphene
per the unit cell of the moiré pattern (see also Section IIIB). The dependence
of the amplitude $\Delta U$ of PES corrugations on the number $N_{c}$ is shown
in Fig. 2(f). The amplitude $\Delta U$ decreases nearly exponentially with the
decrease of the size of the unit cell of the PES or, equivalently, with the
increase of the unit cell of the moiré pattern. The analogous exponential
decrease of the amplitude of PES corrugations with the increase of the unit
cell of the moiré pattern was observed previously for rigid finite graphene
layers Xu  _et al._ (2013) where the rim contribution to static friction is
dominant. Note that extremely low values of the amplitude $\Delta U$ for the
moiré patterns with the smaller sizes of the PES unit cell make evident the
necessity of using classical potentials for the PES calculations.

###  III.2 Approximation of PES by the first Fourier harmonics

PESs of interlayer interaction in diverse hexagonal 2D materials can be
closely approximated by the expressions containing only the first spatial
Fourier harmonics determined by the system symmetry. The adequacy of such an
approximation was demonstrated for coaligned graphene layers Ershova  _et al._
(2010); Lebedeva  _et al._ (2011a); Popov  _et al._ (2012a); Zhou  _et al._
(2015); Reguzzoni  _et al._ (2012), hexagonal boron nitride (h-BN) Lebedev
_et al._ (2016); Zhou  _et al._ (2015), hydrofluorinated graphene Lebedev  _et
al._ (2020), graphene/h-BN heterostructure Jung  _et al._ (2015); Kumar  _et
al._ (2015); Lebedev  _et al._ (2017) and double-layer graphene with krypton
spacer Popov  _et al._ (2013).

These approximations are based on the following considerations. The
translational symmetry of the PES for an atom adsorbed on a triangular lattice
is the same as of the triangular lattice, that is
$U_{a}(\bm{\mathrm{r}})=U_{a}(\bm{\mathrm{r}}+n_{1}\bm{\mathrm{a}}_{1}+n_{2}\bm{\mathrm{a}}_{1})$,
where $\bm{\mathrm{a}}_{1}$ and $\bm{\mathrm{a}}_{2}$ are the lattice vectors
($|\bm{\mathrm{a}}_{1}|=|\bm{\mathrm{a}}_{2}|=a$ and the angle between the
vectors is 60∘, Fig. 3), for any integer $n_{1}$ and $n_{2}$. This means that
the Fourier transform of $U_{a}(\bm{\mathrm{r}})$ consists of harmonics
corresponding to vertices of the reciprocal lattice with the lattice vectors
$\bm{\mathrm{b}}_{1}$ and $\bm{\mathrm{b}}_{2}$ such that
$\bm{\mathrm{a}}_{i}\cdot\bm{\mathrm{b}}_{j}=2\pi\delta_{ij}$
($|\bm{\mathrm{b}}_{1}|=|\bm{\mathrm{b}}_{2}|=4\pi/\sqrt{3}a$ and the angle
between these vectors is 120∘, Fig. 3). Taking into account only the first
Fourier harmonics with wavevectors $\bm{\mathrm{b}}_{1}$,
$\bm{\mathrm{b}}_{2}$ and $\bm{\mathrm{b}}_{1}+\bm{\mathrm{b}}_{2}$, the PES
for an atom on a triangular lattice can be approximated as

$\begin{split}\delta
U_{a}(\bm{\mathrm{r}})=U_{a,1}\mathrm{Re}\bigg{[}e^{i\bm{\mathrm{b}}_{1}\bm{\mathrm{r}}}+e^{i\bm{\mathrm{b}}_{2}\bm{\mathrm{r}}}+e^{i(\bm{\mathrm{b}}_{1}+\bm{\mathrm{b}}_{2})\bm{\mathrm{r}}}\bigg{]},\end{split}$
(1)

where $\delta U_{a}$ is the deviation from the average interaction energy
between the atom and lattice and point $\bm{\mathrm{r}}=0$ corresponds to the
case when the atom is located on top of one of the lattice atoms.

Figure 3: Triangular lattice (orange dots) and its reciprocal lattice (blue
dots). The lattice vectors $\bm{\mathrm{a}}_{1}$ and $\bm{\mathrm{a}}_{2}$
($|\bm{\mathrm{a}}_{1}|=|\bm{\mathrm{a}}_{2}|=a$) as well as reciprocal
lattice vectors $\bm{\mathrm{b}}_{1}$ and $\bm{\mathrm{b}}_{2}$ are shown.
Graphene honeycomb lattice is shown by grey lines. Coordinate $x$ corresponds
to zigzag direction of graphene layer.

For $x$ and $y$ axes chosen along one of the lattice vectors
($\bm{\mathrm{a}}_{1}$) and in the perpendicular direction Verhoeven  _et al._
(2004),

$\begin{split}\delta
U_{a}(x,y)=U_{a,1}\Big{[}2\cos{(k_{y}y)}\cos{(k_{x}x)}+\cos{(2k_{y}y)}\Big{]},\end{split}$
(2)

where $k_{x}=2\pi/a$ and $k_{y}=2\pi/\sqrt{3}a$.

To get the potential energy of an atom on a honeycomb lattice, it is needed to
sum up the expressions for two sublattices separated by $a/\sqrt{3}$ along the
$y$ axis (in the armchair direction). This leads to a sign change for the
second term of Eq. (2). For two coaligned honeycomb layers, it is needed to
sum up once more the contributions of two sublattices of the adsorbed layer.
This leads to another sign change in the equation. Thus, the PES of coaligned
honeycomb layers is described the equation similar to Eq. (2) (Refs. Ershova
_et al._ , 2010; Lebedeva  _et al._ , 2011a; Popov  _et al._ , 2012a).

Now let us consider twisted honeycomb lattices. Let $\bm{\mathrm{a}}_{1}$ and
$\bm{\mathrm{a}}_{2}$ be the lattice vectors of the bottom layer and
$\bm{\mathrm{a}}^{\prime}_{1}$ and $\bm{\mathrm{a}}^{\prime}_{2}$ of the upper
one. The PES of twisted honeycomb lattices is periodic with respect to
translation along any of these lattice vectors:
$U(\bm{\mathrm{r}})=U(\bm{\mathrm{r}}+n_{1}\bm{\mathrm{a}}_{1}+n_{2}\bm{\mathrm{a}}_{2})=U(\bm{\mathrm{r}}+n^{\prime}_{1}\bm{\mathrm{a}}^{\prime}_{1}+n^{\prime}_{2}\bm{\mathrm{a}}^{\prime}_{2})$,
where $n_{1}$, $n_{2}$, $n^{\prime}_{1}$ and $n^{\prime}_{2}$ are any integer
numbers. Thus, the harmonics that contribute to the Fourier transform if the
PES of twisted honeycomb lattices should comply with the condition
$\bm{\mathrm{G}}=m_{1}\bm{\mathrm{b}}_{1}+m_{2}\bm{\mathrm{b}}_{2}=m^{\prime}_{1}\bm{\mathrm{b}}^{\prime}_{1}+m^{\prime}_{2}\bm{\mathrm{b}}^{\prime}_{2}$,
where $\bm{\mathrm{b}}_{1}$ and $\bm{\mathrm{b}}_{2}$ are the lattice vectors
of the reciprocal lattice of the bottom layer, $\bm{\mathrm{b}}^{\prime}_{1}$
and $\bm{\mathrm{b}}^{\prime}_{2}$ of the upper one and $m_{1}$, $m_{2}$,
$m^{\prime}_{1}$ and $m^{\prime}_{2}$ are some integer numbers. This means
these Fourier harmonics correspond to overlapping vertices of the reciprocal
lattices of the twisted layers. The reciprocal lattices of the twisted
honeycomb layers are also two twisted honeycomb lattices forming a
commensurate moiré pattern similar to one in real space. Therefore, the first
Fourier terms contributing to the PES in this case correspond to the lattice
vectors of this moiré pattern of the reciprocal lattices. They have the length
$B=bL/a$, where $L$ is the period of the moiré pattern, and are rotated with
respect to the reciprocal lattice vectors by the same angle $\varphi$ as the
moiré pattern vectors are rotated with respect to the lattice vectors in real
space. As a result, the PES of twisted honeycomb layers can be approximated in
the same form as Eq. (2) if we consider $x^{\prime}$ and $y^{\prime}$ axes
aligned along one of the moiré pattern vectors and in the perpendicular
direction as well as increase the wavevectors by the factor of
$L/a=\sqrt{N_{c}}$:

$\delta
U(x^{\prime},y^{\prime})=U_{1}\bigg{(}2\cos{(k^{\prime}_{y}y^{\prime})}\cos{(k^{\prime}_{x}x^{\prime})}+\cos{(2k^{\prime}_{y}y^{\prime})}\bigg{)},$
(3)

where $\delta U=U-U_{\text{av}}$ is the deviation from the average PES energy
$U_{\text{av}}={}$const, $x^{\prime}=x\cos\varphi-y\sin\varphi$,
$y^{\prime}=y\cos\varphi+x\sin\varphi$, $k^{\prime}_{x}=\sqrt{N_{c}}k_{x}$ and
$k^{\prime}_{y}=\sqrt{N_{c}}k_{y}$. Thus the PES of interlayer interaction of
an infinite graphene bilayer with a commensurate moiré pattern has the same
shape as the PES of for graphene bilayer with coaligned layers (presented in
Ref. Ershova  _et al._ , 2010; Lebedeva  _et al._ , 2011a; Popov  _et al._ ,
2012a; Reguzzoni  _et al._ , 2012) and differs only by the PES amplitude and
the period which is lower by the factor of $\sqrt{N_{c}}$.

The derived Eq. (3) is used here for the approximation of the calculated PESs.
To fit the amplitude $\Delta U$ of PES corrugations, the _single_ parameter of
the approximation is chosen as $U_{1}=2\Delta U/9$. The relative deviation
$\epsilon$ of the approximated and computed PESs is found as the root-mean-
square deviation divided by $\Delta U$. As can be seen in Table 2, this
relative deviation varies from 0.02% to the maximum of 1% for all the
considered moiré patterns. The relative deviations in previous studies of PESs
of coaligned layers are about 1% for graphene Popov  _et al._ (2012a);
Lebedeva  _et al._ (2011a), 0.1–0.3% for h-BN Lebedev  _et al._ (2016), 0.3%
for graphene/h-BN heterostructure Lebedev  _et al._ (2017) and 3% for
hydrofluorinated graphene Lebedev  _et al._ (2020). We believe that the simple
shape of the PES obtained here is a universal property of commensurate twisted
bilayers consisting of various 2D materials analogous to that of coaligned
commensurate bilayers.

It should be noted that for coaligned graphene layers, $U_{a,1}$ in Eq. (2) is
positive and corresponds to the repulsion between the atoms of the upper and
lower layers.Ershova  _et al._ (2010); Lebedeva  _et al._ (2011a); Popov  _et
al._ (2012a); Zhou  _et al._ (2015); Reguzzoni  _et al._ (2012) In the case of
a moiré pattern, the PES is determined by the sum of contributions from many
atoms within the moiré pattern unit cell and $U_{1}$ in Eq. (3) can be
positive or negative (Fig. 2 and Fig. S1 in Supplemental Material) depending
in which symmetry points there is more repulsion between the layers.

Let us discuss the discrepancy of the PESs obtained here using the
Kolmogorov–Crespi potential and the results of the calculations using the
Lebedeva potential, where no corrugations which exceed the calculation
accuracy are observed for the (2,1) moiré pattern Minkin  _et al._ (2021,
2022). The PES of interlayer interaction for coaligned graphene layers
obtained by the DFT-D calculations is excellently approximated by the first
Fourier harmonics Ershova  _et al._ (2010); Lebedeva  _et al._ (2011a); Popov
_et al._ (2012a); Reguzzoni  _et al._ (2012). The Lebedeva potential was
specifically designed to reproduce this property of the PES so that the
relative root-mean-square deviation of the approximated and computed PESs is
within several percents for the Lebedeva potential. A similar deviation for
the Kolmogorov–Crespi potential is 20 times greater Lebedeva  _et al._
(2011a). The greater deviation of the Kolmogorov–Crespi potential is related
with considerably larger amplitudes of other Fourier harmonics except the
first one (which exactly reproduces the approximated PES) including those
compatible with the symmetry of commensurate twisted graphene bilayers and,
therefore, responsible for the shape of corresponding PESs. This explains the
discrepancy of the results obtained using the Kolmogorov–Crespi and Lebedeva
potentials.

Eqs. (2) and (3) are derived based on the system symmetry and make sense not
only for the interlayer interaction energy but also for other properties of 2D
materials. For example, Eq. (2) was used to approximate the interlayer
tunneling contribution to the Hamiltonian of graphene layers in Ref. Jung  _et
al._ , 2014. Considering the Brillouin zone for twisted layers, the authors
derived for them the explicit Hamiltonian and analyzed their band structures.
A similar study was also performed for graphene/h-BN heterostructure Jung  _et
al._ (2014, 2015). The use of Eq. (3) taking into account the symmetry of
commensurate twisted bilayers might be used to further simplify such models.

###  III.3 Properties related with PES

A number of physical properties associated with relative in-plane motion of
the layers are determined by the PES at a constant interlayer distance Popov
_et al._ (2012a); Lebedev  _et al._ (2016, 2017, 2020). Since the PES is
described by a simple expression involving just one energetic parameter [see
Eq. (3)], all these properties can be described analytically as functions of
this parameter. Below we use the PES for the considered moiré patterns to
estimate the shear mode frequency, shear modulus, shear strength and barrier
for relative rotation of the commensurate twisted layers to an incommensurate
state.

Table 2: Approximation parameters $U_{1}$, relative root-mean-square deviations $\epsilon$, shear mode frequencies $f$, shear moduli $C_{44}$, shear strengths $\tau$ and barriers $\Delta U_{\mathrm{rot}}$ for relative rotation of the commensurate twisted layers to an incommensurate state estimated for different moiré patterns based on calculations with the Kolmogorov–Crespi potential. moiré pattern |  $U_{1}$ (eV/atom111per atom of the upper layer) | $\epsilon$ |  $f$ (cm-1) |  $C_{44}$ (Pa) |  $\tau$ (Pa) |  $\Delta U_{\mathrm{rot}}$ (eV/atoma)
---|---|---|---|---|---|---
(2,1) |  $-$2.02$\cdot 10^{-5}$ | 1.04$\cdot 10^{-2}$ | 9.153 | 3.91$\cdot 10^{8}$ | 1.67$\cdot 10^{7}$ | 6.07$\cdot 10^{-5}$
(3,1) | 4.66$\cdot 10^{-6}$ | 5.34$\cdot 10^{-3}$ | 4.233 | 8.37$\cdot 10^{7}$ | 1.12$\cdot 10^{6}$ | 6.99$\cdot 10^{-6}$
(3,2) | 4.34$\cdot 10^{-7}$ | 3.33$\cdot 10^{-3}$ | 1.561 | 1.14$\cdot 10^{7}$ | 1.26$\cdot 10^{5}$ | 6.50$\cdot 10^{-7}$
(5,1) | 1.15$\cdot 10^{-9}$ | 2.04$\cdot 10^{-3}$ | 0.103 | 4.91$\cdot 10^{4}$ | 4.25$\cdot 10^{2}$ | 1.72$\cdot 10^{-9}$

The frequency of the shear mode $E_{2g}$, in which adjacent layers slide
rigidly in the opposite in-plane directions, can be found from the PES
curvature in a given energy minimum Popov  _et al._ (2012a); Lebedev  _et al._
(2016, 2017, 2020) as

$\begin{split}f=\frac{1}{2\pi}\sqrt{\frac{1}{\mu}\frac{\partial^{2}U}{\partial
x^{\prime 2}}}=\frac{1}{a}\sqrt{\frac{1}{\mu}U_{\mathrm{eff}}},\end{split}$
(4)

where $U_{\mathrm{eff}}=(a/2\pi)^{2}\partial^{2}U/\partial x^{\prime 2}$ is
the second-order derivative of the energy per carbon atom of the upper layer
in energy units and $\mu$ is the reduced mass. The latter for bilayer graphene
is as $\mu=m_{\mathrm{C}}/2$, where $m_{\mathrm{C}}$ is the mass of a carbon
atom.

From Eq. (3), it follows that the PES curvature corresponds to
$U_{\mathrm{eff}}=N_{c}U_{1}=(2/9)N_{c}\Delta U$ for moiré patterns with
$U_{1}>0$ and $-2N_{c}U_{1}$ for $U_{1}<0$. The shear mode frequencies
estimated for different moiré patterns using the values of the parameter
$U_{1}$ derived from the calculations with the Kolmogorov–Crespi potential are
listed in Table 2. They all are within 10 cm-1 and are considerably smaller
than the shear mode frequency for the coaligned graphene layers: 35 cm-1 (Ref.
Lebedeva  _et al._ , 2011a, 2012) and 21–34 cm-1 (Ref. Lebedeva  _et al._ ,
2017) according to DFT calculations and $28\pm 3$ cm-1 (Ref. Boschetto  _et
al._ , 2013) and 32 cm-1 (Ref. Tan  _et al._ , 2012) according to the
experiments. Indeed, as follows from Eq. (4), the frequency depends on the
square root $\sqrt{\Delta UN_{c}}$ of the product of the amplitude $\Delta U$
of PES corrugations and the number $N_{c}$ of unit cells of graphene per the
unit cell of the commensurate moiré pattern. Since $\Delta U$ decreases
exponentially with the growth of $N_{c}$ [Fig. 2(f)], the shear mode frequency
also gets reduced upon increasing $N_{c}$.

The PES curvature also determines the shear modulus Lebedev  _et al._ (2016,
2020)

$\begin{split}C_{44}=\frac{d}{\sigma}\frac{\partial^{2}U}{\partial x^{\prime
2}}=\frac{16\pi^{2}d}{\sqrt{3}a^{4}}U_{\mathrm{eff}},\end{split}$ (5)

where $\sigma=\sqrt{3}a^{2}/4$ is the area per carbon atom and $d=3.46$ Å is
the interlayer distance. The shear moduli estimated for different moiré
patterns do not exceed 0.4 GPa (Table 2). This is an order of magnitude
smaller than the DFT result for the coaligned graphene bilayer of 3.8–4.1 GPa
(Ref. Lebedeva  _et al._ , 2017) because the shear modulus is proportional to
$\Delta UN_{c}$.

The PES also determines the static friction force $f_{s}$ for moving the
layers as a whole, i.e. the maximal first derivative of the potential energy
along the minimum energy path (MEP) between adjacent energy minima and,
correspondingly, the shear strength $\tau$ related to it as
$\tau=f_{s}/\sigma$ (here the force is taken per atom of the upper layer).
Analogous estimates of shear strength have been performed first for
commensurate double-walled nanotubes Bichoutskaia  _et al._ (2009b). For
$U_{1}>0$, the MEP between adjacent energy minima corresponds to the line
$x^{\prime}=0$ and $y^{\prime}$ from $a/\sqrt{3N_{c}}$ to $2a/\sqrt{3N_{c}}$
in Eq. (3). The force along the MEP is given by

$-\frac{\partial U}{\partial
y^{\prime}}\Big{|}_{x^{\prime}=0}=2k^{\prime}_{y}U_{1}\bigg{(}\sin{(k^{\prime}_{y}y^{\prime})}+\sin{(2k^{\prime}_{y}y^{\prime})}\bigg{)}$
(6)

and the force extrema are determined by the equation
$\cos{(k^{\prime}_{y}y^{\prime})}+2\cos{(2k^{\prime}_{y}y^{\prime})}=0$. The
latter equation gives that the maximal absolute force is achieved for
$\cos{(k^{\prime}_{y}y^{\prime})}=-(1+\sqrt{33})/8$ for the considered MEP.
From this, we find that the shear strength is
$\tau=6.183\sqrt{N_{c}}U_{1}/a^{3}$ for $U_{1}>0$.

For $U_{1}<0$, the MEP corresponds to $y^{\prime}=0$ and $x^{\prime}$ from $0$
to $a/\sqrt{N_{c}}$ in Eq. (3). In this case,

$-\frac{\partial U}{\partial
x^{\prime}}\Big{|}_{y^{\prime}=0}=2k^{\prime}_{x}U_{1}\sin{(k^{\prime}_{x}x^{\prime})}.$
(7)

The maximal absolute force is achieved for $x^{\prime}=a/4$. Correspondingly,
the shear strength is $\tau=16\pi\sqrt{N_{c}}U_{1}/\sqrt{3}a^{3}$.

The shear strength $\tau$ values estimated for different moiré patterns are
within 0.02 GPa (Table 2). From typical DFT values for the amplitude of PES
corrugations for coaligned graphene layers of about 15 meV per atom of the
upper layer Popov  _et al._ (2012a); Lebedeva  _et al._ (2017), we deduce that
the shear strength in that case should be about 0.22 GPa, i.e. an order of
magnitude greater. As seen from the above equations, the shear strength for
moiré patterns is proportional to $\sqrt{N_{c}}\Delta U$.

When the graphene layers are rotated with respect to each other by an
arbitrary angle that does not correspond to any commensurate moiré pattern,
the area contribution to the PES vanishes and the PES of an infinite
incommensurate twisted bilayer becomes flat. Therefore, the interaction energy
in such a fully incommensurate state can be found as an average over the PES:
$U_{\mathrm{rot}}=\langle U\rangle_{x,y}$ (Refs. Popov  _et al._ , 2012a;
Lebedev  _et al._ , 2016, 2017, 2020). The barrier $\Delta U_{\mathrm{rot}}$
for relative rotation of the layers to a fully incommensurate state can, thus,
be obtained by subtracting the energy in the minimum from $U_{\mathrm{rot}}$.
From Eq. (3), one gets $\Delta U_{\mathrm{rot}}=1.5U_{1}$ for $U_{1}>0$ and
$-3U_{1}$ for $U_{1}<0$. The values of the barrier estimated for the moiré
patterns considered are within 0.06 meV per atom of the upper layer (Table 2).
Obviously they are much smaller than the previous predictions for the
coaligned graphene bilayer of 4 meV/atom (Refs. Lebedeva  _et al._ , 2010,
2011b) and 5 meV/atom (Ref. Popov  _et al._ , 2012a).

Structural superlubricity can be lost via rotation of the layers with the same
lattice constant to the commensurate ground state with coaligned layers Hirano
and Shinjo (1990); Verhoeven  _et al._ (2004); Dienwiebel  _et al._ (2005);
Filippov  _et al._ (2008); Bonelli  _et al._ (2009); Guo  _et al._ (2007);
Shibuta and Elliott (2011); Xu  _et al._ (2013); Wang  _et al._ (2019); Feng
_et al._ (2013). The robust superlubricity has been recently achieved for
systems with a lattice mismatch such as heterostructures composed of layers of
different 2D materials Song  _et al._ (2018) or layers of the same 2D material
under different tension applied Wang  _et al._ (2019); Androulidakis  _et al._
(2020). For such robust superlubric systems, the relative rotation of the
layers to a commensurate interface with the loss of superlubricity is not
possible. Here we propose that the robust superlubricity can be also achieved
for systems in which rotation of the layers to the commensurate ground state
with coaligned layers is possible but hindered by a barrier. For twisted
commensurate layers, such a rotation should occur through a fully
incommensurate state. Although the barriers between a local minimum of twisted
commensurate state and a fully incommensurate state calculated here are rather
small, they might result sufficient to prevent the relative rotation of layers
for a sufficiently large contact area and thus can ensure the macroscopic
robust superlubricity. Further macroscale investigations would be needed to
confirm this hypothesis.

It should be kept in mind that the values of the physical quantities obtained
in the present Section and listed in Table 2 are based on the calculations
with the Kolmogorov-Crespi potential Kolmogorov and Crespi (2005), while its
adequacy for twisted graphene layers has not been proven. Nevertheless, the
potential gives a reasonable dependence of the amplitude $\Delta U$ of PES
corrugations on the size of the unit cell of the moiré pattern [Fig. 2(f)] and
thus should properly describe the trend in the evaluated physical quanities
for different moiré patterns. Once a more reliable potential for twisted
graphene layers is available, it can be used to obtain more accurate estimates
based on the equations given above. On the other hand, as soon as any of these
physical quantities is accessed experimentally (by analogy with the
measurements for coaligned graphene layers Boschetto  _et al._ (2013); Tan
_et al._ (2012)), classical potentials can be refined to improve the
description of the PES of the twisted layers on the basis of the above
formalism. We also believe that the simple shape of the PES can be a universal
property for commensurate twisted bilayers consisting of diverse 2D materials.
This means that physical properties of moiré patterns of other 2D materials
can be estimated in a similar way.

The results of the geometrical analysis persented here are derived without
taking into account structural relaxation. Our calculations show that the
account of structural relaxation does not lead to changes in the shape of PES
and it is still described by the first Fourier harmonics (whereas some
increase of the amplitude $\Delta U$ of PES corrugations occurs). Thus, all
the equations presented here are still valid. The influence of the structural
relaxation on structural and tribological properties of commensurate moiré
patterns will be considered elsewhere.

##  IV Discussion and conclusions

PESs of interlayer interaction have been calculated for a set of commensurate
moiré patterns of twisted graphene bilayer using the registry-dependent
Kolmogorov–Crespi potential. The amplitude of PES corrugations is found to
exceed the calculation accuracy only for 5 moiré patterns with the smaller
unit cell sizes. All calculated PESs have the same simple shape which
corresponds to the symmetry of commensurate moiré patterns and with the size
of the unit cell of PES which is inversely related to the unit cell size of
the moiré pattern. The amplitude of PES corrugations exponentially decreases
with increasing the unit cell size of the moiré pattern. An analytical
expression which is based on the first Fourier harmonics describing the PES
has been derived. The calculated PESs can be approximated by the derived
expression with the accuracy within 1% relative to the amplitude of PES
corrugations. Since the derived expression contains a single energetic
parameter, it has been used to estimate a set of physical quantities
determined by the PES such as shear mode frequency, shear modulus, shear
strength and barrier for relative rotation of the commensurate twisted layers
to a fully incommensurate state. We propose that the latter barrier might
prevent the rotation of the layers from the twisted commensurate state to the
ground commensurate state through a fully incommensurate state and, therefore,
can possibly lead to the macroscopic robust superlubricity for a sufficiently
large contact area.

The approximation by the first Fourier harmonics can be applied not only for
consideration of the interlayer interaction energy. For example, such an
approximation for coaligned layers was used for the anaysis of electronic
properties of twisted graphene Jung  _et al._ (2014) and graphene/h-BN
heterostructure Jung  _et al._ (2014, 2015). Moreover, since the approximation
of the PES of interlayer interaction by the first Fourier harmonics is a
universal property for coaligned layers of diverse 2D materials Lebedev  _et
al._ (2020), we believe that the simple shape of the PES obtained here for
twisted commensurate graphene bilayer can be also universal for any
commensurate moiré patterns consisting of layers of diverse 2D materials.

The raw data on calculated PES required to reproduce our findings are
available to download from Ref. Minkin  _et al._ , 2023.

## Acknowledgments

A.S.M., A.M.P. and Y.E.L. acknowledges the support by the Russian Science
Foundation grant No. 23-42-10010, https://rscf.ru/en/project/23-42-10010/, for
the results described in subsections IIIA “PES of twisted graphene bilayer”
and IIIB “Approximation of PES by the first Fourier harmonics”. I.V.L.
acknowledges Bikaintek grant “Transport” from the Basque Government. A.M.P.
and Y.E.L. acknowledges the support by project FFUU-2021-0003 of the Institute
of Spectroscopy of the Russian Academy of sciences for the results described
in subsection IIIC “Properties related with PES”. S.A.V. and N.A.P.
acknowledge support by the Belarusian Republican Foundation for Fundamental
Research (Grant No. F23RNF-049) and by the Belarusian National Research
Program “Convergence-2025”. This work has been particularly carried out using
computing resources of the federal collective usage center Complex for
Simulation and Data Processing for Mega-science Facilities at NRC “Kurchatov
Institute”, http://ckp.nrcki.ru.

The authors declare no conflict of interest.

## References

  * Hirano and Shinjo (1990) M. Hirano and K. Shinjo, Phys. Rev. B 41, 11837 (1990).
  * Hirano  _et al._ (1991) M. Hirano, K. Shinjo, R. Kaneko, and Y. Murata, Phys. Rev. Lett. 67, 2642 (1991).
  * Hod  _et al._ (2018) O. Hod, E. Meyer, Q. Zheng, and M. Urbakh, Nature 563, 485 (2018).
  * Verhoeven  _et al._ (2004) G. S. Verhoeven, M. Dienwiebel, and J. W. M. Frenken, Phys. Rev. B 70, 165418 (2004).
  * Dienwiebel  _et al._ (2005) M. Dienwiebel, N. Pradeep, G. S. Verhoeven, H. W. Zandbergen, and J. W. M. Frenken, Surf. Sci. 576, 197 (2005).
  * Filippov  _et al._ (2008) A. E. Filippov, M. Dienwiebel, J. W. M. Frenken, J. Klafter, and M. Urbakh, Phys. Rev. Lett. 100, 046102 (2008).
  * Xu  _et al._ (2013) Z. Xu, X. Li, B. I. Yakobson, and F. Ding, Nanoscale 5, 6736 (2013).
  * Koren and Duerig (2016) E. Koren and U. Duerig, Phys. Rev. B 94, 045401 (2016).
  * Bonelli  _et al._ (2009) F. Bonelli, N. Manini, E. Cadelano, and L. Colombo, Eur. Phys. J. B 70, 449 (2009).
  * van Wijk  _et al._ (2013) M. M. van Wijk, M. Dienwiebel, J. W. M. Frenken, and A. Fasolino, Phys. Rev. B 88, 235423 (2013).
  * Guo  _et al._ (2007) Y. Guo, W. Guo, and C. Chen, Phys. Rev. B 76, 155429 (2007).
  * Shibuta and Elliott (2011) Y. Shibuta and J. A. Elliott, Chem. Phys. Lett. 512, 146 (2011).
  * Wang  _et al._ (2019) K. Wang, W. Ouyang, W. Cao, M. Ma, and Q. Zheng, Nanoscale 11, 2186 (2019).
  * Zhang  _et al._ (2015) H. Zhang, Z. Guo, H. Gao, and T. Chang, Carbon 94, 60 (2015).
  * Zhang and Chang (2018) H. Zhang and T. Chang, Nanoscale 10, 2447 (2018).
  * Zhang  _et al._ (2022) H. Zhang, J. Qu, Z. Guo, L. Huang, and Q. Xie, AIP Adv. 12, 115312 (2022).
  * Androulidakis  _et al._ (2020) C. Androulidakis, E. N. Koukaras, G. Paterakis, G. Trakakis, and C. Galiotis, Nat. Commun. 11, 1595 (2020).
  * Liu  _et al._ (2012) Z. Liu, J. Yang, F. Grey, J. Z. Liu, Y. Liu, Y. Wang, Y. Yang, Y. Cheng, and Q. Zheng, Phys. Rev. Lett. 108, 205503 (2012).
  * Vu  _et al._ (2016) C. C. Vu, S. Zhang, M. Urbakh, Q. Li, Q.-C. He, and Q. Zheng, Phys. Rev. B 94, 081405(R) (2016).
  * Song  _et al._ (2018) Y. Song, D. Mandelli, O. Hod, M. Urbakh, M. Ma, and Q. Zheng, Nature Mater. 17, 894 (2018).
  * Mandelli  _et al._ (2017) D. Mandelli, I. Leven, O. Hod, and M. Urbakh, Sci. Rep. 7, 10851 (2017).
  * Minkin  _et al._ (2021) A. S. Minkin, I. V. Lebedeva, A. M. Popov, and A. A. Knizhnik, Phys. Rev. B 104, 075444 (2021).
  * Minkin  _et al._ (2022) A. S. Minkin, I. V. Lebedeva, A. M. Popov, and A. A. Knizhnik, Nanotechnol. Russia 17, 472 (2022).
  * Campanera  _et al._ (2007) J. M. Campanera, G. Savini, I. Suarez-Martinez, and M. I. Heggie, Phys. Rev. B 75, 235449 (2007).
  * Lebedeva and Popov (2020) I. V. Lebedeva and A. M. Popov, J. Phys. Chem. C 124, 2120 (2020).
  * Popov  _et al._ (2011) A. M. Popov, I. V. Lebedeva, A. A. Knizhnik, Y. E. Lozovik, and B. V. Potapkin, Phys. Rev. B 84, 045404 (2011).
  * Mele (2012) E. J. Mele, J. Phys. D: Appl. Phys. 45, 154004 (2012).
  * Kabengele and Johnson (2021) T. Kabengele and E. R. Johnson, Nanoscale 13, 14399 (2021).
  * Kolmogorov and Crespi (2000) A. N. Kolmogorov and V. H. Crespi, Phys. Rev. Lett. 85, 4727 (2000).
  * Belikov  _et al._ (2004) A. V. Belikov, Y. E. Lozovik, A. G. Nikolaev, and A. M. Popov, Chem. Phys. Lett. 385, 72 (2004).
  * Bichoutskaia  _et al._ (2006) E. Bichoutskaia, M. I. Heggie, A. M. Popov, and Y. E. Lozovik, Phys. Rev. B 73, 045435 (2006).
  * Lebedev  _et al._ (2020) A. V. Lebedev, I. V. Lebedeva, A. M. Popov, A. A. Knizhnik, N. A. Poklonski, and S. A. Vyrko, Phys. Rev. B 102, 045418 (2020).
  * Ershova  _et al._ (2010) O. V. Ershova, T. C. Lillestolen, and E. Bichoutskaia, Phys. Chem. Chem. Phys. 12, 6483 (2010).
  * Lebedeva  _et al._ (2011a) I. V. Lebedeva, A. A. Knizhnik, A. M. Popov, Y. E. Lozovik, and B. V. Potapkin, Phys. Chem. Chem. Phys. 13, 5687 (2011a).
  * Popov  _et al._ (2012a) A. M. Popov, I. V. Lebedeva, A. A. Knizhnik, Y. E. Lozovik, and B. V. Potapkin, Chem. Phys. Lett. 536, 82 (2012a).
  * Lebedeva  _et al._ (2012) I. V. Lebedeva, A. A. Knizhnik, A. M. Popov, Y. E. Lozovik, and B. V. Potapkin, Physica E 44, 949 (2012).
  * Zhou  _et al._ (2015) S. Zhou, J. Han, S. Dai, J. Sun, and D. J. Srolovitz, Phys. Rev. B 92, 155438 (2015).
  * Reguzzoni  _et al._ (2012) M. Reguzzoni, A. Fasolino, E. Molinari, and M. C. Righi, Phys. Rev. B 86, 245434 (2012).
  * Lebedev  _et al._ (2016) A. V. Lebedev, I. V. Lebedeva, A. A. Knizhnik, and A. M. Popov, RSC Adv. 6, 6423 (2016).
  * Popov  _et al._ (2013) A. M. Popov, I. V. Lebedeva, A. A. Knizhnik, Y. E. Lozovik, B. V. Potapkin, N. A. Poklonski, A. I. Siahlo, and S. A. Vyrko, J. Chem. Phys. 139, 154705 (2013).
  * Jung  _et al._ (2015) J. Jung, A. M. DaSilva, A. H. MacDonald, and S. Adam, Nat. Commun. 6, 6308 (2015).
  * Kumar  _et al._ (2015) H. Kumar, D. Er, L. Dong, J. Li, and V. B. Shenoy, Sci. Rep. 5, 10872 (2015).
  * Lebedev  _et al._ (2017) A. V. Lebedev, I. V. Lebedeva, A. M. Popov, and A. A. Knizhnik, Phys. Rev. B 96, 085432 (2017).
  * Vuković  _et al._ (2003) T. Vuković, M. Damnjanović, and I. Milošević, Physica E 16, 259 (2003).
  * Bichoutskaia  _et al._ (2005) E. Bichoutskaia, A. M. Popov, A. El-Barbary, M. I. Heggie, and Y. E. Lozovik, Phys. Rev. B 71, 113403 (2005).
  * Bichoutskaia  _et al._ (2009a) E. Bichoutskaia, A. M. Popov, Y. E. Lozovik, O. V. Ershova, I. V. Lebedeva, and A. A. Knizhnik, Phys. Rev. B 80, 165427 (2009a).
  * Popov  _et al._ (2009) A. M. Popov, Y. E. Lozovik, A. S. Sobennikov, and A. A. Knizhnik, JETP 108, 621 (2009).
  * Popov  _et al._ (2012b) A. M. Popov, I. V. Lebedeva, and A. A. Knizhnik, Appl. Phys. Lett. 100, 173101 (2012b).
  * Kolmogorov and Crespi (2005) A. N. Kolmogorov and V. H. Crespi, Phys. Rev. B 71, 235415 (2005).
  * Jung  _et al._ (2014) J. Jung, A. Raoux, Z. Qiao, and A. H. MacDonald, Phys. Rev. B 89, 205414 (2014).
  * Lebedeva  _et al._ (2017) I. V. Lebedeva, A. V. Lebedev, A. M. Popov, and A. A. Knizhnik, Comput. Mater. Sci. 128, 45 (2017).
  * Boschetto  _et al._ (2013) D. Boschetto, L. Malard, C. H. Lui, K. F. Mak, Z. Li, H. Yan, and T. F. Heinz, Nano Lett. 13, 4620 (2013).
  * 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, Nature Mater. 11, 294 (2012).
  * Bichoutskaia  _et al._ (2009b) E. Bichoutskaia, O. V. Ershova, Y. E. Lozovik, and A. M. Popov, Tech. Phys. Lett. 35, 666 (2009b).
  * Lebedeva  _et al._ (2010) I. V. Lebedeva, A. A. Knizhnik, A. M. Popov, O. V. Ershova, Y. E. Lozovik, and B. V. Potapkin, Phys. Rev. B 82, 155460 (2010).
  * Lebedeva  _et al._ (2011b) I. V. Lebedeva, A. A. Knizhnik, A. M. Popov, O. V. Ershova, Y. E. Lozovik, and B. Potapkin, J. Chem. Phys. 134, 104505 (2011b).
  * Feng  _et al._ (2013) X. Feng, S. Kwon, J. Y. Park, and M. Salmeron, ACS Nano 7, 1718 (2013).
  * Minkin  _et al._ (2023) A. S. Minkin, I. V. Lebedeva, A. M. Popov, S. A. Vyrko, N. A. Poklonski, and Y. E. Lozovik, “Data for: Interlayer interaction, shear vibrational mode, and tribological properties of two-dimensional bilayers with a commensurate moiré pattern,” Mendeley Data (2023), doi: 10.17632/svnpzbjc3t.