# Buckling in a rotationally invariant spin-elastic model

Gregorio García-Valladares<EMAIL_ADDRESS>Carlos A. Plata<EMAIL_ADDRESS>Antonio Prados<EMAIL_ADDRESS>Física Teórica, Universidad de Sevilla, Apartado
de Correos 1065, E-41080 Sevilla, Spain

###### Abstract

Scanning tunneling microscopy experiments have revealed an spontaneous
rippled-to-buckled transition in heated graphene sheets, in absence of any
mechanical load. Several models relying on a simplified picture of the
interaction between elastic and internal, electronic, degrees of freedom have
been proposed to understand this phenomenon. Nevertheless, these models are
not fully consistent with the classical theory of elasticity, since they do
not preserve rotational invariance. Herein, we develop and analyse an
alternative classical spin-elastic model that preserves rotational invariance
while giving a qualitative account of the rippled-to-buckled transition. By
integrating over the internal degrees of freedom, an effective free energy for
the elastic modes is derived, which only depends on the curvature.
Minimisation of this free energy gives rise to the emergence of different
mechanical phases, whose thermodynamic stability is thoroughly analysed, both
analytically and numerically. All phases are characterised by a spatially
homogeneous curvature, which plays the role of the order parameter for the
rippled-to-buckled transition, in both the one- and two-dimensional cases. In
the latter, our focus is put on the honeycomb lattice, which is representative
of actual graphene.

††preprint: APS/123-QED

##  I Introduction

Behaviour of low-dimensional systems is a prolific research field [1, 2, 3, 4,
5, 6, 7, 8, 9]. Two-dimensional materials exhibit intriguing properties [10,
11, 12], such as flexibility, high thermal and electrical conductivity,
transparency, low Joule effect, to name just a few. These relevant features
can be useful for many applications, from realisation of microsystems for
healthcare [13, 14] up to energy harvesting [9, 15, 16]. For a recent review
on the effects of strains in graphene and other two dimensional materials, see
Ref. [17].

An outstanding phenomenon in the context of mechanical properties of elastic
systems is the buckling transition. In structural engineering, buckling is the
abrupt change of the considered structure, from a flat to a non-zero curvature
(buckled) state, under load. A prototypical example occurs when a long column,
i.e. a one-dimensional system, buckles under axial compression, which is known
as Euler buckling [18, 19]. Nevertheless, not only is buckling an engineering
subject but also an appealing physical behaviour to be better understood from
a theoretical standpoint. Moreover, phenomena similar to buckling emerge in a
wide range of systems, from low-dimensional systems [4, 20, 21, 22] to cells
[23] or in device development [24].

Graphene sheets are not perfectly two-dimensional. In fact, they display
transversal displacements (ripples) due to thermal fluctuations [25, 26] that
stabilise the membrane 111It is important to recall that the Mermin-Wagner
Theorem [66] forbids the existence of purely two-dimensional stable crystals.
Also, the instability of two-dimensional crystals can be understood within the
context of Landau’s theory of fluctuations [67].. Rippling has been studied
employing models where its origin and behaviour are linked with the coupling
between electronic and elastic degrees of freedom, i.e. the electron-phonon
coupling [1, 26, 28, 29, 30, 31, 32]. Rippling has also been investigated by
employing phenomenological models such as the Helfrich model [33], originally
developed for membranes [34, 35]. In this approach, the free energy of the
system is written as a functional of the configuration of the membrane,
specifically the free energy depends on its curvatures—both mean and Gaussian.

In slender structures such as graphene, there also appears a buckling-like
transition from a flat (or rippled, when taking into account thermal
fluctuations) state to a buckled state under mechanical load [4, 21]. In
addition, scanning tunneling microscopy (STM) experiments have shown that,
even in absence of mechanical load, local heating induces a transition from a
soft rippled sheet to a hard buckled graphene membrane [36, 37, 38, 39],
similar to the Euler buckling transition in structural engineering. This is
somehow unexpected, since the increase of temperature makes the membrane go
from a less ordered state (the flat or rippled one) to a more ordered
(buckled) state. Due to the strong resemblance between this transition and
Euler’s buckling under stress, and following the terminology usually employed
in previous literature, we also refer to this phenomenon as buckling.

The complexity of the interactions in real low-dimensional elastic systems
makes it difficult to obtain analytical results, either exactly or with
controlled approximations, that improve our understanding of the observed
transitions. Therefore, mesoscopic models are relevant, since they contain—in
a simplified way—the main ingredients of low-dimensional elastic systems and
allow for a detailed analytical approach. Specifically, models based on the
coupling between elastic degrees of freedom and Ising spins, termed spin-
string and spin-membrane in the one-dimensional and two-dimensional cases,
respectively, are able to reproduce in a qualitative way the aforementioned
transition. Despite their simplicity, these mesoscopic models present a rich
behaviour, with a complicated phase diagram that include buckled phases [31,
32, 40].

Nevertheless, the currently available spin-string and spin-membrane models
[31, 32, 40], although giving an appealing physical picture, are not
completely satisfactory from a fundamental point of view. First, some of the
microscopic parameters controlling the interactions have to depend on the
system size to ensure that the thermodynamic properties have the correct
behaviour in the large system size limit, e.g. an intensive transition
temperature and an extensive free energy. Second, they are not fully
consistent with the classical theory of elasticity [18], since rotational
symmetry is broken in absence of external forces.

In this work, we propose the simplest model that keeps the appealing physical
picture of the models in Refs. [31, 32, 40] and mends the above mentioned
fundamental flaws. We consider a classical elastic lattice (one or two-
dimensional) characterised by the out-of-plane deformations and internal
degrees of freedom described by a pseudospin variable 222We use the term
pseudospin to emphasise that it is not the real spin but a two-state variable
that mimics in a simple way the range of values taken by the internal degrees
of freedom.. Elastic deformations are modelled by a harmonic interaction
between the atoms that depends on the curvature of the membrane, the repulsive
interaction between electrons is modelled by an antiferromagnetic interaction
between pseudospins, and the electron-phonon coupling is modelled by the
coupling between pseudospins and the membrane curvature. No external stress is
acting on the system. In this framework, we will show that the equilibrium
state of the system can be completely characterised by a free energy that
depends on the membrane curvature—in close analogy with the phenomenological
approach based on the Helfrich model [34, 35, 33]. The equilibrium value of
the curvature only depends on two parameters: temperature $T$ and
antiferromagnetic coupling constant $J$.

In this work, we will prove that equilibrium curvature is spatially
homogeneous and can be interpreted as the order parameter of the rippled to
buckled transition. This is true both for the one-dimensional and the two-
dimensional lattices, regardless of the exact geometry of the two-dimensional
lattice—as long as it does not contain triangular loops. The phase diagram of
the system is built in the $(J,T)$ plane, with an approach that combines
analytical and numerical results. This allows us to characterise both the
existence and stability of all the possible phases.

In the spin-string model (one-dimensional case), the knowledge of the
partition function for the pseudospins make it possible to derive an
analytical expression for the free energy for all values of $(J,T)$.
Therefrom, many analytical results can be derived. In particular, the
transition lines are obtained analytically. In absence of antiferromagnetic
coupling, $J=0$, the transition from the rippled, zero-curvature, phase
(termed ZC) to the buckled phase (termed B) is second-order, at a certain
temperature $T_{0}$. As the antiferromagnetic coupling is increased, the
second-order transition from rippled to buckled moves to lower temperatures
and at a certain critical point $(J_{c},T_{c})$ the transition changes to
first order. Beyond this tricritical point, basically for low enough
temperatures in a given range of $J$’s, the behaviour of the system becomes
more complex: there appear two buckled phases, one of them locally stable
(termed B+) and the other unstable (termed B-). Moreover, there is
metastability: the rippled phase ZC is also locally stable—corresponding to a
local minimum of the free energy, and which phase, B+ or ZC, gives the
absolute minimum of the free energy depends on the value of the
antiferromagnetic coupling. This picture is supported by the analytical
results that can be worked out close to the second-order bifurcation curve and
in the limit of very low temperatures, and also by the numerical construction
of the complete phase diagram of the system in the $(J,T)$ plane.

In the spin-membrane system (two-dimensional case), there is not an analytical
closed form of the partition function for the pseudospins. Still, some
analytical results can be obtained, specifically in the limits of absence of
antiferromagnetic coupling among pseudospins or very low temperatures—for
lattices without triangular loops in the latter case. Despite the increase in
the dimensionality of the system, the qualitative picture is similar to the
spin-string system. For $J=0$, there appears a second-order transition from a
rippled phase (ZC) to a buckled phase (B). In the limit $T\to 0^{+}$, the
transition is first-order and there is metastability: again, two buckled
phases B$\pm$ emerge, one locally stable and the other unstable, and the zero-
curvature phase ZC is also locally stable. In complete analogy with the one-
dimensional situation, the relative stability of the B+ and ZC phases depends
on the intensity of the antiferromagnetic coupling.

The paper is organised as follows. In Sec. II, we define the one-dimensional
model and study its equilibrium behaviour, by deriving a free energy
functional from the probability of a given profile. In Sec. III, we introduce
the Euler-Lagrange equation, which provides us with the equilibrium profile,
and define dimensionless variables. In Sec. IV, we analyse, both analytically
and numerically, the emergence of buckled states and put forward the main
results of our study. In Sec. V, we generalise our one-dimensional model to
the two-dimensional case, focusing ourselves on the specific case of a
honeycomb lattice—that of actual graphene. We present the main conclusions of
our work and perspectives for future research in Sec. VI. The appendices deal
with some non-essential technical details that are omitted in the main text.

##  II Rotationally invariant spin-string model

We consider a string on a one-dimensional lattice, with lattice constant $a$.
At each site, $j=0,\ldots,N$, there is a particle of mass $m$ that is
characterised by its transversal displacement $u_{j}$, and, in addition, there
is a “pseudospin” variable $\sigma_{j}=\pm 1$. We introduce the following
rotationally invariant energy (or Hamiltonian) for any configuration of the
elastic modes $\bm{u}$, their associated momenta $\bm{p}$, and the pseudospins
$\bm{\sigma}$,

$\displaystyle\mathcal{H}(\bm{u},\bm{p},\bm{\sigma})=\sum_{j=0}^{N}$
$\displaystyle\bigg{[}\frac{p_{j}^{2}}{2m}+\frac{k}{2}\left(u_{j+1}-2u_{j}+u_{j-1}\right)^{2}$
$\displaystyle-$ $\displaystyle
h\left(u_{j+1}-2u_{j}+u_{j-1}\right)\sigma_{j}+J\sigma_{j+1}\sigma_{j}\bigg{]}.$
(1)

The first term on the rhs is the kinetic energy, the second one stands for the
elastic contribution to the energy—with elastic constant $k$, the third one
involves the interaction between the elastic displacements and the
pseudospins—the parameter $h$ tunes the strength thereof, and the last one
represents the nearest-neighbour interaction among the pseudospins—with
coupling constant $J$. Our Hamiltonian (1) can be analysed for both positive
and negative values of the coupling constant $J$. In the case of our concern,
the interaction between pseudospins models the repulsive interaction between
electrons, as discussed in detail below, and thus we focus on the
antiferromagnetic case $J>0$ 333Yet, at some specific points in the paper we
briefly consider the situation for ferromagnetic interaction, $J<0$.. The
pseudospins do not correspond to actual spin variables, they model in a simple
way internal degrees of freedom—as discussed below.

The above system is a mesoscopic toy model, which tries to resemble in a
simple manner some of the main features of low-dimensional elastic systems,
such as graphene. Therein, the interaction between the elastic modes and the
pseudospins would model the electron-phonon coupling 444More specifically, the
coupling between the transverse elastic modes and out-of-plane electrons,
which can be thought to be in two states corresponding to the two different
sides of the plane., and the antiferromagnetic interaction between the
pseudospins account for the Coulomb repulsion of electrons. Therefore,
nearest-neighbour out-of-plane electrons prefer to be at different sides of
the plane. There are two main differences with respect to other proposals in
the literature [31, 32, 40]: (i) the elastic term is proportional to the
square of the (discrete) curvature and not to the square of the (discrete)
gradient, and (ii) the coupling between the elastic modes and the pseudospins
also involves the curvature and not the transverse displacement.

The key role played by the curvature of the system, as shown below, makes our
model consistent with the classical theory of elasticity [18]—at variance with
previous spin-string (or spin-membrane in the two-dimensional case) models
[31, 32, 40]. In particular, a flat profile, $u_{j}\equiv 0$, and a rigid
rotation thereof, $u_{j}=Aj$ with constant $A$, have the same energy in Eq.
(1). In general, any two profiles differing just by a linear function $Aj+B$
has the same energy. The case $A=0$ corresponds to a rigid translation,
whereas the case $A\neq 0$ corresponds to a small rotation of angle $A$,
$\sin(A)\simeq A$.

Thermal equilibrium of the system at temperature $T$ is described by the
canonical distribution
$P_{\text{eq}}(\bm{u},\bm{p},\bm{\sigma})=\exp\left[-\beta\mathcal{H}(\bm{u},\bm{p},\bm{\sigma})\right]/Z$,
where $\beta=(k_{B}T)^{-1}$, with $k_{B}$ and $Z$ being the Boltzmann constant
and the partition function, respectively. We are interested in the equilibrium
profiles of the string, so we integrate over $\bm{p}$ and $\bm{\sigma}$ to get
$P_{\text{eq}}(\bm{u})\propto e^{-\beta\mathcal{F}(\bm{u})}$, where the free
energy of the string $\mathcal{F}(\bm{u})$ is

$\mathcal{F}(\bm{u})=\sum_{j=0}^{N}\left[\frac{k}{2}\left(u_{j+1}-2u_{j}+u_{j-1}\right)^{2}\right]-k_{B}T\ln
Z_{\text{int}}(\bm{u}).$ (2)

Above, $Z_{\text{int}}(\bm{u})$ stands for the partial partition function
associated with the internal degrees of freedom, which involves a sum over all
the pseudospin configurations

$Z_{\text{int}}=\sum_{\bm{\sigma}}e^{-\beta\sum_{j=0}^{N}\left[-h\left(u_{j+1}-2u_{j}+u_{j-1}\right)\sigma_{j}+J\sigma_{j+1}\sigma_{j}\right]}.$
(3)

Now we introduce a continuum limit by assuming that the displacements $u_{j}$
vary slowly with $j$, i.e., $u_{j}\to u(x)$ with $x=ja$. The total length of
the system is $L=Na$. In this way,

$u_{j+1}-2u_{j}+u_{j-1}=a^{2}\chi(x)+O(a^{4}),$ (4)

where we have defined the curvature of the profile $\chi(x)\equiv
u^{\prime\prime}(x)$, and

$\displaystyle\frac{k}{2}(u_{j+1}-2u_{j}+u_{j-1})^{2}=\frac{k_{0}}{2}\chi^{2},$
(5a) $\displaystyle h(u_{j+1}-2u_{j}+u_{j-1})\sigma_{j}=h_{0}\chi\sigma_{j}.$
(5b)

The continuum limit in the spin variables is skipped willingly, since a
marginal sum over all spin configurations will be carried out in brief. Scaled
parameters for the elastic constant and the strength of the pseudospin-elastic
interaction have been introduced as

$k_{0}=ka^{4},\quad h_{0}=ha^{2}.$ (6)

Also, terms with higher powers of $a$ in Eq. (4) have been omitted, since they
involve higher-order derivatives of $u$ 555Therefore, these terms are expected
to be subdominant when $u_{j}$ varies slowly with $j$. The (small) corrections
introduced thereby could be incorporated into our theoretical framework in a
perturbative way..

In the continuum limit, the equilibrium probability of a certain string
profile $u(x)$ becomes a functional thereof, specifically

$P_{\text{eq}}[u]\propto
e^{-\beta\mathcal{F}[u]},\quad\mathcal{F}[u]=n\int_{0}^{L}dxf(\chi),$ (7)

where $n\equiv N/L$ is the number density, and the free energy density (per
particle) $f$ is given by

$\displaystyle f(\chi)=$
$\displaystyle\,\frac{k_{0}\chi^{2}}{2}-\beta^{-1}\ln\zeta_{\text{int}}^{\text{1d}}(\beta
h_{0}\chi,\beta J),$ (8a) $\displaystyle\zeta_{\text{int}}^{\text{1d}}(\beta
h_{0}\chi,\beta J)\equiv$ $\displaystyle\,e^{-\beta J}\cosh\left(\beta
h_{0}\chi\right)$ $\displaystyle+e^{\beta J}\sqrt{1+e^{-4\beta
J}\sinh^{2}\left(\beta h_{0}\chi\right)},$ (8b)

Our notation explicitly tells us that $f$ only depends on the profile through
its curvature $\chi$. The second term in $f(\chi)$ comes from the sum of the
spin variables, $\zeta_{\text{int}}^{\text{1d}}$ is the one-particle partition
function of a one-dimensional Ising chain with coupling $J$ and external field
$h_{0}\chi$—e.g. see Ref. [45].

In contrast to previous spin-string and spin-membrane models [31, 32, 40], the
continuum limit consider here does not involve a large system size limit $N\gg
1$ (thermodynamic limit). This approach to the continuum limit is more
consistent from a physical point of view. More specifically, the microscopic
parameters in the Hamiltonian, such as $h$ and $k$ do not scale with the
system size—thus decoupling the continuum and the thermodynamic limits.

##  III Equilibrium profiles

The equilibrium profiles $u_{\text{eq}}(x)$ are those that maximise the
probability distribution $P_{\text{eq}}[u]$ or, equivalently, those that
minimise the free energy functional $\mathcal{F}[u]$. Therefore, we consider
the first variation of the free energy functional upon the change $u\to
u+\delta u$,

$\displaystyle\delta F[u]=$ $\displaystyle\int_{0}^{L}\\!\\!\\!d{x}\,\delta
u\,\frac{d^{2}}{dx^{2}}\\!\left(\partialderivative{f}{\chi}\right)\\!\\!+\\!\left[\delta
u^{\prime}\,\partialderivative{f}{\chi}-\delta
u\,\derivative{x}\\!\left(\partialderivative{f}{\chi}\right)\right]_{0}^{L}.$
(9)

The expression for $\delta F[u]$ is more involved than in the usual case where
the integrand only depends on the first derivative, but the integral and the
boundary terms must vanish separately nonetheless [46, 47].666For instance,
for free boundary conditions, the variation $\delta u(x)$ for $x\in(0,1)$ and
$\\{\delta u(0),\delta u(1),\delta u^{\prime}(0),\delta u^{\prime}(1)\\}$ at
the boundaries are independent. Considering the vanishing of the integral
term, one obtains the Euler-Lagrange equation

$\left.\derivative[2]{x}\\!\left(\partialderivative{f}{\chi}\right)\right|_{\text{eq}}=0\implies\left.\partialderivative{f}{\chi}\right|_{\text{eq}}=Ax+B,$
(10)

where $A$ and $B$ are arbitrary constants, to be determined by imposing the
boundary conditions. The boundary conditions depend on the physical situation:
here, for the sake of simplicity we are going to consider that the ends of the
string are fixed, so that

$u(0)=u(L)=0,$ (11)

but the value of $u^{\prime}$ at the boundaries is free—the so-called
supported boundary conditions [18]. Therefore, $\partialderivative*{f}{\chi}$
must vanish at both boundaries for the equilibrium profile, which entails
$A=B=0$.

The equilibrium profile is thus provided by
$\left.\partialderivative{f}{\chi}\right|_{\text{eq}}\\!\\!\\!=0$, i.e.

$\displaystyle k_{0}\chi_{\text{eq}}$
$\displaystyle=\left.k_{B}T\partialderivative{\ln\zeta_{\text{int}}^{\text{1d}}}{\chi}\right|_{\chi=\chi_{\text{eq}}}\\!\\!\\!$
$\displaystyle=h_{0}\frac{e^{-\frac{2J}{k_{B}T}}\sinh{\left(\frac{h_{0}\chi_{\text{eq}}}{k_{B}T}\right)}}{\sqrt{1+e^{-\frac{4J}{k_{B}T}}\sinh^{2}{\left(\frac{h_{0}\chi_{\text{eq}}}{k_{B}T}\right)}}},$
(12)

with the boundary conditions (11). Equation (III) is a transcendental equation
for the equilibrium curvature $\chi_{\text{eq}}$, the solution of which can be
written as $\chi_{\text{eq}}=\chi_{\text{eq}}(T,J,k_{0},h_{0})$.777Zero
curvature profiles, i.e. linear profiles with $\chi=0$, always solve Eq.
(III). The equilibrium curvature is constant, independent of position, and so
is $f_{\text{eq}}=f(\chi_{\text{eq}})$. Therefore, the equilibrium free energy
is $\mathcal{F}_{\text{eq}}=nLf_{\text{eq}}=Nf_{\text{eq}}$, i.e. it is
extensive—note that the constant $\chi_{\text{eq}}$ is an intensive quantity,
independent of system size.

If $\chi_{\text{eq}}$ is a solution of Eq. (III), then $-\chi_{\text{eq}}$ is
also a solution. This is a consequence of the free energy density being an
even function of $\chi$, and it is clear on a physical basis: the specular
reflection of an equilibrium profile with respect to the $x$ axis is also an
equilibrium profile, since there is no external field. Therefore, without loss
of generality, we restrict ourselves to solutions with positive curvature in
the remainder of this work.

To better characterise the equilibrium curvature $\chi_{\text{eq}}$, it is
adequate to introduce dimensionless variables: this allows us to identify the
natural units of length and energy—or temperature—in our model system, and
elucidate the dimensionless combinations of the system parameters on which the
curvature depends. Equation (III) tells us the characteristic length
$\ell_{0}$, which in turn determines the characteristic temperature $T_{0}$,

$\ell_{0}=\frac{k_{0}}{h_{0}}=\frac{k}{h}a^{2},\qquad
T_{0}=\frac{h_{0}^{2}}{k_{B}k_{0}}=\frac{h^{2}}{k_{B}k}.$ (13)

Then we define dimensionless displacement $u^{*}$, position $x^{*}$, curvature
$\chi^{*}$, temperature $T^{*}$, and coupling constant $J^{*}$ by

$u^{*}=\frac{u}{\ell_{0}},\,\,x^{*}=\frac{x}{\ell_{0}},\,\,\chi^{*}=\ell_{0}\chi,\,\,T^{*}=\frac{T}{T_{0}},\,\,J^{*}=\frac{J}{k_{B}T_{0}}.$
(14)

Consistently, we define $\beta^{*}=k_{B}T_{0}\beta=1/T^{*}$,
$L^{*}=L/\ell_{0}$, $n^{*}=n\ell_{0}=N/L^{*}$, and the dimensionless free
energy $\mathcal{F}^{*}\equiv\mathcal{F}/k_{B}T_{0}$. Then, we can write

$\displaystyle\mathcal{F}^{*}[u^{*}](J^{*},T^{*})$
$\displaystyle=n^{*}\int_{0}^{L^{*}}dx^{*}f^{*}(\chi^{*};J^{*},T^{*}),$ (15a)
$\displaystyle f^{*}(\chi^{*};J^{*},T^{*})$
$\displaystyle\equiv\frac{1}{2}{\chi^{*}}^{2}-T^{*}\ln\zeta_{\text{int}}^{\text{1d}}(\beta^{*}\chi^{*},\beta^{*}\\!J^{*}),$
(15b)

where $f^{*}$ is the dimensionless free energy density and
$\zeta_{\text{int}}^{\text{1d}}$ is the same one-particle partition function
for the pseudospins than that in Eq. (8b).888Note that $\beta
h_{0}\chi=\beta^{*}\chi^{*}$ and $\beta J=\beta^{*}J^{*}$.

Our use of dimensionless variables simplifies the theoretical analysis, by
reducing the number of parameters—and identifying the relevant combinations
thereof. This course of action aligns well with our main goal, which is the
explanation of the experimentally observed buckling phenomenon with a simple
model. A quantitative comparison of the results predicted by our model with
experimental data of a specific system is beyond the scope of this work, and
in any case would involve fitting the model parameters to the experimental
data.

In dimensionless variables, Eq. (III) is rewritten as
$\left.\partialderivative{f^{*}}{\chi^{*}}\right|_{\text{eq}}=0$, i.e.

$\displaystyle\chi_{\text{eq}}^{*}$
$\displaystyle=\mu^{*}(\beta^{*}\chi_{\text{eq}}^{*},\beta^{*}\\!J^{*}),$
(16a) $\displaystyle\mu^{*}(\beta^{*}\chi^{*},\beta^{*}\\!J^{*})$
$\displaystyle\equiv\partialderivative{\ln\zeta_{\text{int}}^{\text{1d}}}{(\beta^{*}\chi^{*})}=\\!\frac{e^{-{2\beta^{*}\\!J^{*}}}\sinh{\left(\beta^{*}\chi^{*}\right)}}{\sqrt{1\\!+\\!e^{-4\beta^{*}\\!J^{*}}\sinh^{2}{\left(\beta^{*}\chi^{*}\right)}}},$
(16b)

where $\mu^{*}$ is the local magnetisation of the pseudospins. Equation (16)
implies that $\chi_{\text{eq}}^{*}$ is a certain function of $J^{*}$ and
$T^{*}$, $\chi_{\text{eq}}^{*}=\chi_{\text{eq}}^{*}(J^{*},T^{*})$. Therefore,
it is the dimensionless pseudospin-elastic coupling $J^{*}$ and the
dimensionless temperature $T^{*}$ that control the emergence of buckled
profiles with non-zero curvature. To keep our notation simple, we do not
explicitly state that both $J^{*}$ and $T^{*}$ depend on the microscopic
parameters of the model through the characteristic temperature $T_{0}$,
defined in Eq. (13), employed to make energy dimensionless. We recall that the
equilibrium curvature is constant, independent of position, and an intensive
quantity—and so is the free energy density. Therefore,

$\mathcal{F}_{\text{eq}}^{*}(J^{*},T^{*})\equiv\mathcal{F}^{*}[u_{\text{eq}}^{*}](J^{*},T^{*})=Nf^{*}(\chi_{\text{eq}}^{*};J^{*},T^{*}),$
(17)

the equilibrium free energy is extensive.

For the supported boundary conditions (11) we are considering, there is a
unique smooth equilibrium profile given by

$u_{\text{eq}}^{*}(x^{*};J^{*},T^{*})=\frac{\chi_{\text{eq}}^{*}(J^{*},T^{*})}{2}x^{*}(x^{*}-L^{*}),$
(18)

For completely free boundary conditions, i.e. when neither the displacement
nor its derivative is fixed at the boundary, there are multiple equilibrium
profiles that are obtained by adding $Cx^{*}+D$ to $u_{\text{eq}}^{*}$ in Eq.
(18), with $C$ and $D$ being arbitrary constants—the boundary condition
$\derivative*{(\partialderivative*{f^{*}}{\chi^{*}})}{x^{*}}|_{x^{*}=0,L^{*}}=0$
is fulfilled for all $(C,D)$. In particular, the flat profile, $u^{*}=0$, and
any transversal shift plus rotation thereof, $u^{*}=Cx^{*}+D$, are both
possible equilibrium profiles for free boundary conditions.

Since the free energy density only depends on the curvature, and the
equilibrium curvature is the same for both supported and completely free
boundary conditions, the analysis of buckled states that follows is valid for
both supported and free boundary conditions. Then, for the sake of
concreteness and without loss of generality, we will employ the supported
boundary condition expression (18) for the equilibrium profile—omitting the
additional linear contribution $Cx^{*}+D$ for the free case.

In the remainder of the paper, we employ dimensionless variables—dropping the
asterisks in them not to clutter our formulae.

##  IV Buckled states

In this section, we analyse the emergence of buckled states, i.e. with non-
zero curvature, in our model system. First, we show how buckled profiles
bifurcate from the zero-curvature solution $\chi_{\text{eq}}^{\text{ZC}}=0$ of
Eq. (16). Second, we consider the low-temperature limit $T\ll 1$, where the
equilibrium profiles can be exactly calculated. Lastly, we analyse the phase
diagram of the model by numerically solving Eq. (16).

###  IV.1 Bifurcation from the zero-curvature solution

To start with, we expand the free energy density around the zero-curvature
profile, for which

$f_{\text{ZC}}(J,T)\equiv f(\chi=0;J,T)=-T\ln\left[2\cosh(J/T)\right],$ (19)

i.e. we consider the difference of the free energy density from the zero-
curvature value,

$\displaystyle\Delta f(\chi;J,T)$ $\displaystyle\equiv
f(\chi;J,T)-f_{\text{ZC}}(J,T)$
$\displaystyle=\frac{1}{2}\chi^{2}-T\ln\left[\frac{\zeta_{\text{int}}^{\text{1d}}(\chi/T,J/T)}{2\cosh(J/T)}\right].$
(20)

The stability of the zero-curvature profile is determined by the sign of the
second-derivative of the free energy density at $\chi=0$,

$f_{2}(J,T)\equiv\left.\partialderivative[2]{f}{\chi}\right|_{\chi=0}=1-\frac{e^{-2J/T}}{T}.$
(21)

The change of stability takes place at the line over which $f_{2}$ vanishes,
i.e.

$J=J_{b}(T)=-\frac{1}{2}T\ln T,$ (22)

which is called the bifurcation curve. This line separates the plane $(J,T)$
into two regions: (i) a region with $f_{2}>0$, comprising two subregions I and
III where the zero-curvature solution is stable and (ii) region II with
$f_{2}<0$, where the zero-curvature solution becomes unstable and other stable
solution exist. Note that for $J=0$, we have that $f_{2}<0$ for $T<1$: the
unit of temperature $T_{0}$ is thus the transition temperature in the absence
of antiferromagnetic coupling among the pseudospins.

The deviation of the free energy density from the zero-curvature value is
provided by the Taylor expansion

$\displaystyle\Delta f(\chi;J,T)=$
$\displaystyle\frac{f_{2}(J,T)}{2}\chi^{2}+\frac{f_{4}(J,T)}{4!}\chi^{4}$
$\displaystyle+\frac{f_{6}(J,T)}{6!}\chi^{6}+O\left(\chi^{8}\right),$ (23)

where

$\displaystyle f_{4}(J,T)=$
$\displaystyle\frac{e^{-6J/T}}{T^{3}}\\!\left(3-e^{4J/T}\right),$ (24a)
$\displaystyle f_{6}(J,T)=$
$\displaystyle-\frac{e^{-10J/T}}{T^{5}}\\!\left(e^{8J/T}-30e^{4J/T}+45\right).$
(24b)

Figure 1: Phase diagram of the rotationally invariant spin-string. The plane
$(J,T)$ is divided into several regions, characterised by the existing phases
and their stability. In region I, only the ZC phase exists. In region II,
phase ZC is unstable, and there appears a stable buckled phase B. In region
III, three phases coexist: ZC, and two buckled phases, B+ and B-; the latter
being unstable, whereas ZC and B+ are locally stable. Phase B+ is the most
stable one in IIIa, with ZC being metastable, whereas the roles are reversed
in IIIb. Point $C=(J_{c},T_{c})$ is the tricritical point. Above it, the phase
transition is second-order; below, it is first-order. The curve $J_{b}(T)$,
given by Eq. (22), is a bifurcation curve. For $T>T_{c}$, it defines the
already mentioned second-order phase transition (black solid line). For
$T<T_{c}$, it demarcates the change from region II to region III (leftmost
dotted line). The curve $J_{M}(T)$ (rightmost dotted line), theoretically
predicted by Eq. (26a) close to the tricritical point (blue dotted line),
demarcates the change from region III to region I. The curve $J_{t}(T)$
(dashed line), theoretically predicted by Eq. (26b) close to the tricritical
point (red dashed line), is the first-order transition line: over it, phases
B+ and ZC have the same free energy. Also plotted are the theoretical low-
temperature limiting values $J_{M}^{(0)}$ and $J_{t}^{(0)}$. On the one hand,
the predictions for the bifurcation curve and for
$\left(J_{t}^{(0)},J_{M}^{(0)}\right)$ are exact. On the other hand, the
asymptotic theoretical predictions for the curves $J_{t}(T)$ and $J_{M}(T)$
show an excellent agreement with the numerical results, within their range of
validity. See Sec. IV.3 for further details in the numerical determination of
the curves $J_{t}(T)$ and $J_{M}(T)$.

The phase diagram of the system is presented in Fig. 1, paying special
attention to the demarcation of the different regions and subregions where
both the number of equilibrium solutions and its stability change. This figure
summarises the main results of the equilibrium and stability analyses carried
out in detail in Appendix A. Therefrom, we highlight here the tricritical
point,

$T_{c}=\frac{1}{\sqrt{3}},\quad J_{c}=J_{b}(T_{c})=\frac{\ln 3}{4\sqrt{3}}$
(25)

and the asymptotic theoretical prediction for the curves $J_{M}$ ($J_{t}$)
separating regions I and IIIb (IIIb and IIIa),

$\displaystyle J_{M}(T)=J_{b}(T)+\frac{5}{12}\sqrt{3}(T-T_{c})^{2}$ (26a)
$\displaystyle J_{t}(T)=J_{b}(T)+\frac{5}{16}\sqrt{3}(T-T_{c})^{2}.$ (26b)

###  IV.2 Low-temperature limit

Now we turn our attention to the low-temperature limit $T\ll 1$. By defining
$e_{\text{gs}}^{\text{1d}}$ as the energy per site in the ground state of the
one-dimensional pseudospin system with external field $\chi$ and
antiferromagnetic coupling $J$, we can write

$\lim_{T\to
0^{+}}-T\ln\zeta_{\text{int}}^{\text{1d}}=e_{\text{gs}}^{\text{1d}}.$ (27)

The ground state energy of the pseudospin system is

$e_{\text{gs}}^{\text{1d}}=-\left(\left|\chi\right|-2J\right)H\left(\left|\chi\right|-2J\right)-J,$
(28)

in which $H(x)$ is Heaviside’s step function. Equation (28) is readily
understood on a physical basis: for small $\left|\chi\right|$, the
antiferromagnetic coupling wins and the ground state corresponds to
antiferromagnetic ordering, $e_{\text{gs}}^{\text{1d}}=-J$, whereas for large
$\left|\chi\right|$, the external field wins and the ground state corresponds
to all spins aligned with the curvature,
$e_{\text{gs}}^{\text{1d}}=-\left|\chi\right|+J$. These two expressions for
$e_{\text{gs}}^{\text{1d}}$ become equal at $\left|\chi\right|=2J$, which
marks the crossover between them both. Of course, $e_{\text{gs}}^{\text{1d}}$
can also be derived by taking the limit $T\to 0^{+}$ in Eq. (27), bringing to
bear Eq. (8b) for $\zeta_{\text{int}}^{\text{1d}}$.

The above discussion entails that the free energy density (15b) simplifies to

$f(\chi;J,T)\sim\frac{\chi^{2}}{2}-\left(\left|\chi\right|-2J\right)H\left(\left|\chi\right|-2J\right)-J,\quad
T\ll 1.$ (29)

Also, Eq. (16) for the equilibrium curvature converts to

$\displaystyle\absolutevalue{\chi_{\text{eq}}}=H(\absolutevalue{\chi_{\text{eq}}}-2J),\quad
T\ll 1,$ (30)

Both in Eqs. (29) and (30), Heaviside’s step function must be understood in a
“physical way”, since it appears here as the low-temperature limit of the rhs
of Eq. (16). Therefore, $H(x)$ comprises three strokes: two flat strokes,
equal to $0$ and $1$ for $x<0$ and $x>0$, respectively, and a vertical stroke
at $x=0$, which joins the previous two. Logically, Eq. (30) preserves the
symmetry $\chi_{\text{eq}}\to-\chi_{\text{eq}}$: we recall that we are
restricting ourselves to positive curvatures, so
$\absolutevalue{\chi_{\text{eq}}}\to\chi_{\text{eq}}$ in the following.

In order to study the solutions of Eq. (30), a graphical method is useful. In
Fig. 2, we plot the functions on the lhs and the rhs of Eq. (30) and look for
the intersection points of both functions. It is clearly seen that
$\chi_{\text{eq}}=0$ is always a solution, for all $J$: the phase ZC survives
in the low-temperature limit. In addition, there appear two buckled solutions
for $2J\leq 1$, i.e. $J\leq J_{M}^{(0)}=1/2$: specifically,
$\chi_{\text{eq}}=2J$ and $\chi_{\text{eq}}=1$. Both solutions coalesce at
$J_{M}^{(0)}$: stronger antiferromagnetic interaction destroys the buckled
phases, in which the magnetisation of the pseudospins is different from zero—a
signature of ferromagnetic ordering.

Figure 2: Graphical method for solving Eq. (30). For each $J$, the
intersection points of the functions $H(\chi_{\text{eq}}-2J)$ (solid lines)
and $\chi_{\text{eq}}$ (dashed blue line) give the equilibrium curvature in
the low-temperature limit. For $J>J_{M}^{(0)}=1/2$, for instance $J=0.55$
(green solid line), the only solution is $\chi_{\text{eq}}=0$, i.e. the ZC
phase. For $J\leq J_{M}^{(0)}$, for instance $J=0.25$ (red solid line), there
appear three solutions: in addition to the ZC phase corresponding to
$\chi_{\text{eq}}=0$, we have two buckled phases corresponding to
$\chi_{\text{eq}}=2J$ (B- phase) and $\chi_{\text{eq}}=1$ (B+ phase).

Let us calculate the free energy density of each of the phases, in order to
identify them. First, the free energy density of the ZC phase is given by
$f_{\text{ZC}}(J,T)\sim-J$, for $T\ll 1$, so that

$\displaystyle\Delta f(\chi;J,T)$
$\displaystyle=f(\chi;J,T)-f_{\text{ZC}}(J,T)$
$\displaystyle\sim\frac{\chi^{2}}{2}-\left(\left|\chi\right|-2J\right)H\left(\left|\chi\right|-2J\right),\quad
T\ll 1.$ (31)

Therefore, one obtains

$\displaystyle\Delta f(\chi_{\text{eq}}=$ $\displaystyle 2J;J,T\to 0)\sim
2J^{2},$ (32a) $\displaystyle\Delta f(\chi_{\text{eq}}=$ $\displaystyle
1;J,T\to 0)\sim-\frac{1}{2}+2J.$ (32b)

On the one hand, the phase with $\chi_{\text{eq}}=2J$ has always a free energy
that is larger than that of the ZC phase—recall that the total free energy is
proportional to its density as stated by Eq. (17). On the other hand, $\Delta
f(\chi_{\text{eq}}=1;J,T\to 0)$ changes sign at $J=J_{t}^{(0)}=1/4$: the phase
with $\chi_{\text{eq}}=1$ has a free energy that is smaller (larger) than that
of the ZC phase for $0\leq J<J_{t}^{(0)}$ $\left(J_{t}^{(0)}<J\leq
J_{M}^{(0)}\right)$. Therefore, the solution with $\chi_{\text{eq}}=1$ is
expected to correspond to the low-temperature limit of phase B+, whereas the
solution with $\chi_{\text{eq}}=2J$ should correspond to the low-temperature
limit of phase B-. Also, the points $\left(J_{t}^{(0)},0\right)$ and
$\left(J_{M}^{(0)},0\right)$ should correspond to the endpoints of the first-
order line $J_{t}(T)$ and the right border line of region III $J_{M}(T)$,
respectively—we have analytical predictions for these lines, them being
rigorously valid just in the vicinity of the tricritical point, Eqs. (26b) and
(26a). These expectations are confirmed by the numerical analysis below.

###  IV.3 Numerical analysis

Here, we perform a numerical analysis of the phase diagram of our model
system. We construct a 200$\times$200 mesh in the $(J,T)$ plane and find
$\chi_{\text{eq}}$, by numerically solving Eq. (16) at each point of this
mesh. At some points $(J,T)$, more than one solution is found: note that this
is expected, since in certain regions of the plane $(J,T)$ several phases
coexist. This analysis allows us to build a numerical phase diagram for our
system and compare it with our analytical predictions above—depicted in Fig.
1.

Figure 3 shows the values of $\chi_{\text{eq}}$ for the thermodynamically
stable phase—the most stable one when there is coexistence. Therein, it is
clearly observed that the transition changes from second-order to first-order
at the tricritical point. Above the tricritical point $C$, $\chi_{\text{eq}}$
changes continously from zero to non-zero values at the bifurcation curve
$J_{b}(T)$ (solid line). Below the tricritical point, $\chi_{\text{eq}}$
changes discontinuously from zero to non-zero values at the first-order line
$J_{t}(T)$ (dashed line). Recall that the curvature plays the role of the
order parameter in our model. In fact, $\chi_{\text{eq}}$ equals the
magnetisation of the pseudospins, which is given by the rhs of Eq. (16).

Figure 3: Curvature of the most stable phase (numerically obtained). The
labelling of the different regions and the code of the different lines is the
same as in Fig. 1. The change of order of the transition, from second-order
(solid line) above $C$ to first-order (dashed line) below, is clearly
observed.  Figure 4: Free energy density difference $\Delta
f(\chi_{\text{eq}};J,T)$ of the buckled states with respect to the ZC phase
(numerically obtained). The labelling of the different regions and the code of
the different lines is the same as in Figs. 1 and Fig. 3. On the left panel,
we show a density plot for the (at least locally) stable buckled phase: B in
region II and B+ in region III. On the right panel, we show the density plot
for the unstable buckled phase B-, which only exists in region III.

By inserting the numerical solution for $\chi_{\text{eq}}$ into Eq. (15b), we
get the values of the free energy density over the considered mesh. We present
our results in Fig. 4, where we plot the difference of the free energy density
with respect to the ZC profile, $\Delta f(\chi_{\text{eq}};J,T)$, as defined
in Eq. (20). In the left panel, the value for the stable buckled phase,
wherever it exists, is shown—in region I, within there is no buckled phase, a
constant value zero is plotted. Therefore, phase B is shown in region II,
inside the bifurcation curve (22) (solid line above the tricritical point $C$,
leftmost dotted line below it), and phase B+ is plotted in region III. Region
III is demarcated by the bifurcation curve, the segment of the $J$-axis
between the origin and the point $\left(J_{M}^{(0)},0\right)$, and the curve
$J_{M}(T)$ that marks the limit of existence of the phases B$\pm$ (rightmost
dotted line). In the right panel, $\Delta f(\chi_{\text{eq}};J,T)$ is plotted
for the unstable phase B-, which only exists in region III—consistently, in
regions I and II, a constant value zero is plotted. The three phases, ZC, B+,
and B- coexist in region III, which is divided into two subregions by the
first-order line $J_{t}(T)$ (dashed line): IIIa, where $\Delta
f(\chi_{\text{eq}}^{\text{B+}};J,T)<0$ and the most stable phase is B+, being
ZC metastable; and IIIb, where the roles are reversed, the most stable phase
is ZC, being B+ metastable.

##  V Spin-membrane model on a honeycomb lattice

Now we move to a two-dimensional system, i.e. a spin-membrane model.
Specifically, we consider the system on a honeycomb lattice, such as that of
graphene. There is a particle—carbon atom if thinking of graphene—at each node
$(i,j)$ in a certain two-dimensional region $\Omega$. The transverse
displacement with respect to the plane of each particle is denoted by
$u_{ij}$. Analogously to the one-dimensional spin-string analysed before,
displacements parallel to the plane are not taken into consideration.
Moreover, at each site $(i,j)$ we have a pseudospin $\sigma_{ij}=\pm 1$, which
models in a simple way other degrees of freedom—for example, the out-of-plane
electron that is not bonded.

A sketch of the lattice honeycomb lattice, with lattice constant $a$, is shown
in Fig. 5. Indices $i$ and $j$ are employed for rows and columns,
respectively. Note that atoms at each row are distributed in zigzag. There are
two types of sites, e-sites and o-sites. For the former, the three nearest
neighbours are one above and two below; for the latter, the three nearest
neighbours are two above and one below [40]. Note that the two-dimensional
domain $\Omega$ can have an arbitrary shape, not necessarily rectangular—for
the sake of simplicity, we assume that it has no holes. This means that all
the nearest neighbours are present for the bulk sites shown in Fig. 5—boundary
sites will be taken into account by introducing appropriate boundary
conditions over the contour $\delta\Omega$ of the region $\Omega$, as
discussed later on.

The rotationally invariant Hamiltonian is the extension of the spin-string
one, as given by Eq. (1), to the honeycomb lattice,

$\displaystyle\mathcal{H}(\hat{u},\hat{p},\hat{\sigma})$
$\displaystyle=\sum_{\absolutevalue{i-j}=\text{even}}\bigg{[}\frac{p_{ij}^{2}}{2m}+\frac{k}{2}\left(u_{i-1,j}+u_{i,j-1}+u_{i,j+1}-3u_{ij}\right)^{2}-h\sigma_{ij}\left(u_{i-1,j}+u_{i,j-1}+u_{i,j+1}-3u_{ij}\right)$
$\displaystyle\qquad\qquad\qquad+\frac{J}{2}\sigma_{ij}\left(\sigma_{i-1,j}+\sigma_{i,j-1}+\sigma_{i,j+1}\right)\bigg{]}$
$\displaystyle+\sum_{\absolutevalue{i-j}=\text{odd}}\bigg{[}\frac{p_{ij}^{2}}{2m}+\frac{k}{2}\left(u_{i+1,j}+u_{i,j-1}+u_{i,j+1}-3u_{ij}\right)^{2}-h\sigma_{ij}\left(u_{i+1,j}+u_{i,j-1}+u_{i,j+1}-3u_{ij}\right)$
$\displaystyle\qquad\qquad\qquad+\frac{J}{2}\sigma_{ij}\left(\sigma_{i+1,j}+\sigma_{i,j-1}+\sigma_{i,j+1}\right)\bigg{]}.$
(33)

For the sake of clarity, the Hamiltonian is split into two sums because
nearest neighbours are different for e-sites and o-sites. Within each sum on
the rhs, the first term stands for the kinetic energy, the second one
corresponds to the elastic contribution to the Hamiltonian, the third one
describes the interaction between the transverse displacements and the
pseudospins, and the fourth one provides the antiferromagnetic coupling among
the pseudospins.999The coupling constant is $J$, the factor $1/2$ is needed
because each interaction is counted twice.

Figure 5: Sketch of the honeycomb lattice with parameter $a$. The lattice is
bipartite and can be thus split into two sublattices: e-sites (blue), two
nearest neighbours below and one above, and o-sites (red), two nearest
neighbours above and one below. The distance between sites of the same kind
corresponding to consecutive indexes are $\Delta x=3a/2$, $\Delta
y=\sqrt{3}a/2$. The spacing between e- and o-sites with the same value of the
index $i$ is $x_{o}=a/2$.

Again, we are interested in the equilibrium profiles of the membrane.
Therefore, we integrate the canonical distribution over the momenta and the
pseudospin variables, firstly introducing the continuum limit—see details in
Appendix B. Now, the curvature is the extension of the one-dimensional case,
i.e., the Laplacian of the displacements $\nabla^{2}u(x,y)$. Then, the sum
over $(i,j)$ goes to an integral over $(x,y)$, which leads to

$P_{\text{eq}}[u]\propto
e^{-\beta\mathcal{F}[u]},\quad\mathcal{F}[u]=n\int_{\Omega}dx\,dy\,f(\chi),$
(34)

where $n=(\Delta x\Delta y)^{-1}$ is the number density and $f$ is the free
energy density per particle, given by

$\displaystyle f(\chi)$
$\displaystyle=\frac{k_{0}\chi^{2}}{2}-k_{B}T\ln\zeta_{\text{int}}^{\text{2d}}(\beta
h_{0}\chi,\beta J),$ (35a)

$\zeta_{\text{int}}^{\text{2d}}(\beta h_{0}\chi,\beta J)$ stands for the
partition function for the antiferromagnetic Ising chain with coupling $J$ and
external field $h_{0}\chi$ in the two-dimensional case. Unfortunately, at
variance with the one-dimensional case, there is not an analytical expression
for $\zeta_{\text{int}}^{\text{2d}}$ when $J\neq 0$.

The equilibrium profile minimises the free energy, so it is determined by the
condition $\delta F[u]=0$. Since the free energy density only depends on the
curvature $\chi$, it is convenient to define the field

$\Phi(x,y)\equiv\partialderivative{f}{\chi}=k_{0}\chi-
h_{0}\partialderivative{\ln\zeta_{\text{int}}^{\text{2d}}(\beta
h_{0}\chi,\beta J)}{(\beta h_{0}\chi)}.$ (36)

Making use of this definition, the variation of the free energy can be written
as

$\displaystyle\delta F[u]=$
$\displaystyle\int_{\Omega}d{x}\,d{y}\,\Phi\,\nabla^{2}(\delta u)$
$\displaystyle=$
$\displaystyle\int_{\Omega}d{x}\,d{y}\,\nabla\\!\cdot\\!\left[\Phi\,\nabla(\delta
u)-\nabla\Phi\delta u\right]+\\!\int_{\Omega}d{x}\,d{y}\,\nabla^{2}\Phi\delta
u$ $\displaystyle=$
$\displaystyle\int_{\delta\Omega}ds\,\hat{\bm{n}}\\!\cdot\\!\left(\Phi\,\nabla(\delta
u)-\nabla\Phi\,\delta u\right)+\int_{\Omega}d{x}\,d{y}\,\nabla^{2}\Phi\,\delta
u.$ (37)

In the contour integral along the boundary $\delta\Omega$, $ds$ is the length
element and $\hat{\bm{n}}$ is the outward pointing unit normal. On the one
hand, the surface integral in the second term on the rhs of Eq. (37) provides
us with the Euler-Lagrange equation,

$\nabla^{2}\Phi=0,$ (38)

i.e. the Laplace equation for the field $\Phi(x,y)$. On the other hand, the
contour integral in the first term gives us the boundary conditions: since
each term must vanish separately, we have

$\Phi\,\partialderivative{n}\delta u=0,\quad\delta
u\,\partialderivative{n}\Phi=0,$ (39)

where $\partial/\partial n\equiv\hat{\bm{n}}\cdot\nabla$ is the normal
derivative. For the sake of concreteness, and consistently with our analysis
of the one-dimensional spin-string model, we consider supported boundary
conditions: the value of the transversal displacements vanish at the
boundaries,

$u(x,y)=0,\quad(x,y)\in\delta\Omega,$ (40)

but the value of the normal derivative $\partialderivative{n}\delta u=0$ at
the boundary is free. Therefore, the coefficient of
$\partialderivative{n}\delta u$ in Eq. (39) has to vanish,

$\Phi(x,y)=0,\quad(x,y)\in\delta\Omega.$ (41)

The solution of the Laplace equation (38) with the homogeneous Dirichlet
boundary condition (41) is identically zero:

$\displaystyle
k_{0}\chi_{\text{eq}}-h_{0}\partialderivative{\ln\zeta_{\text{int}}^{\text{2d}}(\beta
h_{0}\chi_{\text{eq}},\beta J)}{(\beta h_{0}\chi_{\text{eq}})}=0,$ (42)

where we have made use of the definition of $\Phi$. Although we do not have an
explicit expression for $\zeta_{\text{int}}^{\text{2d}}$, Eq. (42) neatly
tells us that the equilibrium curvature is constant, independent of $(x,y)$
over all the two-dimensional domain $\Omega$. Interestingly, the
nondimensional variables introduced in the one-dimensional spin-string model,
as given by Eq. (14), also work in the two-dimensional case—only now the
dimensionless number density is $n^{*}=n\ell_{0}^{2}$. The dimensionless free
energy and free energy density are

$\displaystyle\mathcal{F}^{*}[u^{*}](J^{*},T^{*})$
$\displaystyle=n^{*}\int_{\Omega}dx^{*}\,dy^{*}\,f^{*}(\chi^{*};J^{*},T^{*}),$
(43a) $\displaystyle f^{*}(\chi^{*};J^{*},T^{*})$
$\displaystyle\equiv\frac{1}{2}{\chi^{*}}^{2}-T^{*}\ln\zeta_{\text{int}}^{\text{2d}}(\beta^{*}\chi^{*},\beta^{*}\\!J^{*}),$
(43b)

where we have taken into account that $\beta h_{0}\chi=\beta^{*}\chi^{*}$,
$\beta J=\beta^{*}J^{*}$. In what follows, we drop again the asterisks in the
dimensionless variables to simplify our notation. Equation (42) is written in
dimensionless variables as

$\chi_{\text{eq}}=\mu(\beta\chi_{\text{eq}},\beta
J),\quad\mu(\beta\chi_{\text{eq}},\beta
J)\equiv\partialderivative{\ln\zeta_{\text{int}}^{\text{2d}}(\beta\chi_{\text{eq}},\beta
J)}{(\beta\chi_{\text{eq}})}.$ (44)

Note that $\mu$ is the local magnetisation of the two-dimensional pseudospin
lattice. The equilibrium free energy $\mathcal{F}_{\text{eq}}$ is given by

$\mathcal{F}_{\text{eq}}(J,T)\equiv\mathcal{F}[u_{\text{eq}}](J,T)=Nf(\chi_{\text{eq}};J,T),$
(45)

which is extensive, as expected on a physical basis.

###  V.1 Buckled states

An important novel feature of our model is its rotational invariance. Indeed,
as was the case of the one-dimensional spin-string model, the two-dimensional
spin-membrane model analysed here is invariant under rotations. Let us
consider a certain equilibrium profile $u(x,y)$, which thus solves the Euler-
Lagrange equation (44). Note that $\tilde{u}(x,y)=u(x,y)+Ax+By+C$, is another
possible equilibrium profile with the same free energy, which only depends on
the curvature $\nabla^{2}u(x,y)$. Similarly to the one-dimensional case, the
family of transformations defined by $Ax+By+C$ contains any small two-
dimensional rotation.

To univocally determine the equilibrium profile for the membrane, we have to
bring to bear the remainder of the supported boundary conditions, i.e. Eq.
(40). Therefore, we have to solve Poisson’s equation, with the uniform
curvature playing the role of the density and homogeneous boundary conditions:

$\nabla^{2}u_{\text{eq}}(x,y)=\chi_{\text{eq}},\;(x,y)\in\Omega;\;\;u(x,y)=0,\;(x,y)\in\delta\Omega.$
(46)

The solution for the profile depends on the geometry of the problem. The
simplest situation appears for a circular membrane of radius $R$: therein, the
equilibrium profile only depends on the distance $r$ to the centre of the
membrane and straightforward integration of the Poisson equation gives

$u_{\text{eq}}(r)=\frac{\chi_{\text{eq}}}{4}(r^{2}-R^{2}).$ (47)

For a rectangle of sides $L_{x}$ and $L_{y}$, the solution of the Poisson
equation involves an expansion in eigenfunctions—for details, see Appendix C.

Below, we analytically investigate the behaviour of the two-dimensional spin-
membrane model in some limits of interest, for which the expression for
$\zeta_{\text{int}}^{\text{2d}}$ can be worked out. Specifically, two
situations are considered: (i) the case $(J=0,T\neq 0)$, in which the
antiferromagnetic coupling among the pseudospins disappears, and (ii) the low-
temperature limit $(J\neq 0,T=0)$. In the former, we show that there appears a
second-order phase transition at $T=1$ (in dimensionless variables), below
which the zero curvature membrane is unstable and stable buckled profiles
bifurcate. In the latter, we show that the zero curvature membrane is always
locally stable but membrane buckled states (one locally stable, another
unstable) are present below a certain value $J_{M}^{(0)}$ of the
antiferromagnetic coupling. Moreover, metastability is present, signaling that
the transition from zero curvature to buckled states is now first-order.

####  V.1.1 No antiferromagnetic coupling

Let us first consider the case $J=0$. The one-particle partition
$\zeta_{\text{int}}^{\text{2d}}(J=0,T)$ corresponds to that of a non-
interacting Ising system in an external field $\chi$, i.e.

$\displaystyle
f(\chi;J=0,T)=\frac{1}{2}\chi^{2}-T\ln\left[2\cosh\left(\frac{\chi}{T}\right)\right].$
(48)

and the Euler-Lagrange for the curvature reduces to

$\chi_{\text{eq}}(J=0,T)=\tanh\left[\frac{\chi_{\text{eq}}(J=0,T)}{T}\right].$
(49)

This coincides with the particularisation of Eq. (16) for $J=0$—this is
logical, without antiferromagnetic interaction the transition is of mean field
type for any dimension. Therefore, the ZC phase is the only one for $T>1$ and
is stable, whereas it becomes unstable for $T<1$. At $T=1$, a stable buckled
phase B continuously bifurcates from the ZC phase, namely

$\chi_{\text{eq}}\sim\sqrt{3(1-T)},\quad\text{for}\ T\to 1^{-}.$ (50)

This phase B is the stable one for $T<1$ and continues to exist up to $T=0$.
In fact, Eq. (49) predicts that

$\lim_{T\to 0^{+}}\chi_{\text{eq}}(J=0,T)=1.$ (51)

####  V.1.2 Low-temperature limit

Now we turn to the low-temperature limit $T\to 0^{+}$, with $J\neq 0$. In this
regime, similarly to the situation in the one-dimensional spin-string system,
we have

$\lim_{T\to
0^{+}}-T\ln\zeta_{\text{int}}^{\text{2d}}=e_{\text{gs}}^{\text{2d}}$ (52)

where

$e_{\text{gs}}^{\text{2d}}=-\left(\left|\chi\right|-3J\right)H\left(\left|\chi\right|-3J\right)-\frac{3J}{2}.$
(53)

is the energy per site in the ground state of the two-dimensional spin-
membrane. Again, Eq. (53) for $e_{\text{gs}}^{\text{2d}}$ has a clear physical
interpretation. On the one hand, the ground state energy corresponds to
antiferromagnetic ordering for small $\chi$,
$e_{\text{gs}}^{\text{2d}}=-3J/2$—in the honeycomb lattice, each pseudospin
has three nearest neighbours. On the other hand, the ground state corresponds
to all pseudospins aligned with the curvature for large $\chi$,
$e_{\text{gs}}^{\text{2d}}=-\chi+3J/2$. The crossing between these expressions
for $e_{\text{gs}}^{\text{2d}}$ takes place at $\chi=3J$, which leads to Eq.
(53)—see details in Appendix D.

The free energy density is now

$f(\chi;J,T)\sim\frac{\chi^{2}}{2}-\left(\left|\chi\right|-3J\right)H\left(\left|\chi\right|-3J\right)-\frac{3J}{2},\quad
T\ll 1,$ (54)

and the equilibrium curvature is fixed by

$\displaystyle\absolutevalue{\chi_{\text{eq}}}=H(\absolutevalue{\chi_{\text{eq}}}-3J),\quad
T\ll 1.$ (55)

We recall that, without loss of generality, we are restricting ourselves to
positive curvatures in the paper, thus
$\absolutevalue{\chi_{\text{eq}}}\to\chi_{\text{eq}}$ henceforth. The analysis
of the low-temperature limit follows along the same lines as in the one-
dimensional case, so we present the results in a concise way. The ZC phase is
always a solution of Eq. (55), and two buckled solutions emerge for $3J\leq
1$, i.e. $J\leq J_{M}^{(0)}=1/3$. The two buckled phases have curvatures
$\chi_{\text{eq}}=3J$ and $\chi_{\text{eq}}=1$. For $J>J_{M}^{(0)}$, the
antiferromagnetic interaction is so strong that prevents the emergence of
buckled solutions.

To identify the phases, we evaluate the free energy density over them, the ZC
phase verifies $f_{\text{ZC}}(J,T)\sim-3J/2$, for $T\ll 1$. Therefore

$\displaystyle\Delta f(\chi;J,T)$
$\displaystyle=f(\chi;J,T)-f_{\text{ZC}}(J,T)$
$\displaystyle\sim\frac{\chi^{2}}{2}-\left(\left|\chi\right|-3J\right)H\left(\left|\chi\right|-3J\right),\quad
T\ll 1.$ (56)

and

$\displaystyle\Delta f(\chi_{\text{eq}}=3J;J,T)$
$\displaystyle\sim\frac{9}{2}J^{2},$ (57a) $\displaystyle\Delta
f(\chi_{\text{eq}}=1;J,T)$ $\displaystyle\sim-\frac{1}{2}+3J.$ (57b)

The phase with $\chi_{\text{eq}}=3J$ always has $\Delta f>0$. The phase with
$\chi_{\text{eq}}=1$ has $\Delta f<0$ ($>0$) for $0\leq J<J_{t}^{(0)}=1/6$
$\left(J_{t}^{(0)}<J\leq J_{M}^{(0)}\right)$. Then, the phase with
$\chi_{\text{eq}}=3J$ should correspond to the (unstable) B- phase, whereas
the phase with $\chi_{\text{eq}}=1$ should correspond to the (at least locally
stable) B+ phase. The buckled phases disappear discontinuously at
$J=J_{M}^{(0)}$, which means that the transition is first-order.

The above identification of the phases is physically sensible, since in the
limit $J\to 0$ the B+ phase for $T\ll 1$ is expected to converge to the $T\to
0^{+}$ limit of the B phase found for $J=0$—and, in fact, the curvature of the
latter has the right behaviour, as expressed by Eq. (51). The profile of the
B+ phase for a circular membrane of unit area is shown in Fig. 6; the
anagolous profile for a rectangular membrane can be found in Appendix C.

Figure 6: Low-temperature stable buckled configuration for the spin-membrane
model. The plotted surface corresponds to the absolute value of Eq. (47), for
a unit area circle. The curvature equals unity for all $J\leq
J_{M}^{(0)}=1/3$, corresponding to the most stable state for $0\leq
J<J_{t}^{(0)}=1/6$ and to a metastable state for $J_{t}^{(0)}<J\leq
J_{M}^{(0)}$.

##  VI Conclusions

In this paper, we have introduced novel, rotationally invariant, spin-string
and spin-membrane models. The main role is played by the curvature associated
with the transversal displacements, which not only rules the elastic energy
but also the coupling with internal degrees of freedom. These models give an
appealing physical picture of the buckling transition while keeping
consistency with both thermodynamics and classical theory of elasticity, which
is important from a theoretical perspective. Moreover, the equilibrium
profiles are very simple, with a position-independent curvature, both for the
one-dimensional and the two-dimensional lattices.

The phase diagram in the $(J,T)$ plane—$J$ being the antiferromagnetic
coupling among the pseudospins and $T$ the system temperature—of the one-
dimensional rotationally invariant spin-string model has been characterised
analytically in great detail. It is complex and qualitatively similar to that
found in Ref. [32] and, in particular, buckled phases emerge for low enough
values of the temperature and the antiferromagnetic coupling. In absence of
the latter, there appears a second-order phase transition at a certain
temperature, below which the zero curvature phase becomes unstable and a
buckled B phase is the stable one. In the buckled phase, the pseudospins
magnetisation equals the non-zero curvature and, thus, they predominantly
align with the sign of the curvature. For low temperatures, the ferromagnetic
ordering of the pseudospins is frustrated by the antiferromagnetic coupling,
which destroys the buckled states for high enough coupling constants $J$,
$J>J_{M}^{(0)}$. For $J\leq J_{M}^{(0)}$, there appears coexistence between
the ZC phase and a stable B+ buckled phase, phase B+ is the most stable one
and phase ZC is metastable for small $J$, $0\leq J<J_{t}^{(0)}$, with the
situation being reversed for $J_{t}^{(0)}<J\leq J_{M}^{(0)}$. This qualitative
map of phases is the key ingredient that makes this kind of models
qualitatively explain the transition from a rippled to a buckled state when
the temperature is increased [31, 36].

Despite our qualitative explaining of the experimentally observed rippled to
buckled transition, the model presented here does have some limitations. The
main one stems from our starting from a mesoscopic description. Our focus is
thus the characterisation of the kind of physical interactions that are able
to cause the observed phenomenology, and not a quantitative account of the
mechanical properties of the graphene. In order to achieve this latter goal,
one would have to start from a microscopic description and go to a mesoscopic
model by taking an adequate physical limit, which is outside the scope of the
present work. Nevertheless, it is precisely the mesoscopic nature of our
modelling that allows us to obtain analytical results and a complete
characterisation of the different phases.

To be concrete, and motivated by the geometry of graphene, we have studied the
spin-membrane model in a honeycomb lattice. The analytical expression of the
partition function for the pseudospins is not known in the two-dimensional
case and thus the phase diagram cannot be characterised as completely as in
the one-dimensional case. Nevertheless, our analysis of limiting cases allows
us to state that the same qualitative picture of phases applies for the two-
dimensional spin-membrane model. Still, the two-dimensional values of
$J_{t}^{(0)}$ and $J_{M}^{(0)}$ are smaller—in dimensionless variables—than
the one-dimensional ones. This can be understood as stemming from the
antiferromagnetic interaction being stronger, since the number of nearest
neighbours is larger for $d=2$ than for $d=1$. It should be stressed that the
main conclusions of our study of the two-dimensional case hold for a quite
general lattice, as long as it does not contain triangular loops—such as a
rectangular one.

Previous spin-string and spin-membrane models [31, 32, 40] had several
weaknesses. First, it was necessary to assume the parameters in the system
Hamiltonian to have certain scalings with system size to get consistency with
thermodynamics, specifically extensive free energy and transition temperature
independent of system size. Second, even assuming size-dependent microscopic
parameters, the models were not consistent with classical elasticity theory:
e.g. a flat but rotated string or membrane was not an acceptable equilibrium
profile.

The models developed in this work mend the aforementioned weaknesses. The
microscopic parameters $k$ and $h$ do not depend on system size and determine
natural units for length and temperature, which we have used to introduce
dimensionless variables. From them, we have naturally obtained a system-size-
independent transition temperature and an extensive free energy. Moreover, the
Euler-Lagrange equation is consistent with the classical theory of elasticity,
with the term stemming from the internal interactions involving a biLaplacian
instead of a Laplacian. These two properties hold for both the one-dimensional
and two-dimensional lattices.

Note that there are no external fields in our Hamiltonian, and therefore
buckling here is a purely thermal phenomenon—there is no way of making the
system buckle as a consequence of mechanical actions. For instance, it would
be interesting to consider a prescribed force applied to one extremity of the
spin-string chain (or to the border of the spin-membrane). From the point of
view of statistical mechanics, two different statistical ensembles emerge
(Gibbs-like, controlled force, and Helmholtz-like, controlled length or area),
which are relevant both from a theoretical standpoint and for practical
applications. The continuum limit considered here does not need—as already
discussed, and at variance with previous works—the system to be large.
Therefore, the problem of ensemble equivalence, which has been analysed in
different contexts [52, 53, 54, 55, 56, 57], is also worth studying here—both
ensembles are expected to be equivalent only in the thermodynamic (large
system size) limit.

Also, perspectives for future work are opened. In the context of elasticity of
low-dimensional systems, our work may trigger new approaches to the analysis
of buckled phases in more realistic models. Also, it would be interesting to
consider variants of the two-dimensional model that allow for a more detailed
analytical treatment. In this regard, effective antiferromagnetic interactions
for which the pseudospins partition function is exactly known [58] deserve to
be studied. Moreover, the framework developed here may be useful to look into
the elastic behaviour of systems in which elastic modes are coupled to other,
internal, degrees of freedom—such as biomolecules [59, 60, 57, 61]. A
different but also interesting prospect is the investigation of dynamics,
since the present study only concerns equilibrium configurations. On a
physical basis, one expects the dynamics of the internal degrees of freedom to
be much faster than that of the elastic modes. This separation of time scales
should make it possible—e.g. making use of a Chapman-Enskog expansion [62,
63]—to derive an effective stochastic dynamics, either at the Langevin or
Fokker-Planck level of description, for the elastic modes.

##  Appendix A Bifurcation theory details

In this Appendix, the different regions of the system phase diagram are
derived in detail using the expansion of the free energy, through bifurcation
theory.

Equation (23) allows us to analyse the emergence of buckled ($\chi\neq 0$)
profiles from the zero-curvature solution, using the Landau theory of phase
transitions—with the curvature playing the role of the order parameter. Within
this approximation, the transcendental Euler-Lagrange equation (16) for the
curvature becomes algebraic: $\left.\partialderivative{\Delta
f(\chi;J,T)}{\chi}\right|_{\text{eq}}=0$ implies

$\displaystyle 0=$ $\displaystyle
f_{2}(J,T)\chi_{\text{eq}}+\frac{f_{4}(J,T)}{3!}\chi_{\text{eq}}^{3}+\frac{f_{6}(J,T)}{5!}\chi_{\text{eq}}^{5}.$
(58)

This expansion is truncated at different orders, the sign of the last retained
coefficient $f_{n}$ must be always positive—otherwise, the corresponding
equilibrium distribution $P_{\text{eq}}$ would be non-normalisable.

####  A.0.1  $f_{4}>0$: Second-order phase transition

When $f_{4}>0$, the sign of $f_{2}$ determines the possible phases of the
system. Close to any point $(J_{b},T_{b})$ over the bifurcation curve (22),
$f_{2}$ is small: we explicitly indicate this by introducing a parameter
$0<\varepsilon\ll 1$ such that $f_{2}=\varepsilon\varphi_{2}$, with
$\varphi_{2}=O(1)$. This tells us that the separation of the point $(J,T)$ we
are considering from the bifurcation curve is of order $\varepsilon$: either
$T-T_{b}=O(\varepsilon)$ or $J-J_{b}=O(\varepsilon)$. We can take advantage of
these results to analyse our model in the ferromagnetic case $J<0$. Therein,
no tricritical point stemming from the competition of the interactions emerge,
because $f_{4}(J,T)>0$ for any $J\leq 0$. Therefore, the bifurcation curve
inside the second quadrant $(J<0,T>0)$ defines a typical second-order phase
transition line. In order to clarify how the regions are extended to the
ferromagnetic situation, a small portion of the second quadrant is included in
Fig. 1.

On the one hand, if $\varphi_{2}>0$, the free energy density has only one
minimum at $\chi_{\text{eq}}^{\text{ZC}}=0$. On the other hand, if
$\varphi_{2}<0$, $\chi_{\text{eq}}^{\text{ZC}}=0$ turns out to be a maximum
and there appear two symmetric minima at

$\chi_{\text{eq}}^{\text{B}}=\pm\varepsilon^{1/2}\sqrt{-6\varphi_{2}/f_{4,b}},$
(59)

where $f_{4,b}\equiv f_{4}(J_{b},T_{b})$ is the value of the coefficient
$f_{4}$ over the bifurcation curve,

$f_{4,b}=\frac{3T_{b}^{2}-1}{T_{b}^{2}}.$ (60)

The superindex B in Eq. (59) indicates that this solution corresponds to a
buckled phase, with non-zero curvature. Inside region II of Fig. 1, the stable
profile for the string is thus buckled.

Moreover, the term we have neglected in Eq. (58),
$f_{6}\chi_{\text{eq}}^{5}/6!$ is in fact negligible: it is of order
$\varepsilon^{5/2}$, whereas the first two terms on the rhs are of order
$\varepsilon^{3/2}$—this shows the consistency of the approximations
introduced.

####  A.0.2 $f_{4}<0$ and $f_{6}>0$: Tricritical point and first-order phase
transition

The coefficient $f_{4,b}$ over the bifurcation curve vanishes at the
tricritical point $C\equiv(J_{c},T_{c})$, given by

$T_{c}=\frac{1}{\sqrt{3}},\quad
J_{c}=J_{b}(T_{c})=\frac{\ln{3}}{4\sqrt{{3}}}.$ (61)

Close to the tricritical point, the approximation given by Eq. (59) thus
breaks down. In fact, below the tricritical point—over the leftmost dotted
line in Fig. 1—we have that $f_{4,b}<0$ and the term involving $f_{6}$ in the
expansion of the free energy density must be retained.

Close to the tricritical point, a different scaling for the solutions of Eq.
(58) is needed. We still write $f_{2}=\varepsilon\varphi_{2}$: the terms
involving $f_{2}$ and $f_{6}$ are of the same order if
$\chi_{\text{eq}}=O(\varepsilon^{1/4})$. Then, the terms involving $f_{2}$ and
$f_{4,b}$ are of the same order if $f_{4,b}=O(\varepsilon^{1/2})$—this tells
us the order of the distance of the point $(J_{b},T_{b})$ over the bifurcation
curve to the tricritical point, i.e. either $T_{b}-T_{c}=O(\varepsilon^{1/2})$
or $J_{b}-J_{c}=O(\varepsilon^{1/2})$. Then, we write
$f_{4,b}=\varepsilon^{1/2}\varphi_{4}$, $f_{6,c}=f_{6}(J_{c},T_{c})=36$, and
the dominant balance for the curvature equation in this region is

$\varphi_{2}\Xi+\frac{1}{6}\varphi_{4}\Xi^{3}+\frac{1}{120}f_{6,c}\Xi^{5}=0,\qquad\chi_{\text{eq}}\equiv\varepsilon^{1/4}\Xi.$
(62)

We find the zero-curvature solution $\chi_{\text{eq}}^{\text{ZC}}=0$, and two
buckled solutions

$\displaystyle\chi_{\text{eq}}^{\text{B+}}=$
$\displaystyle\pm\varepsilon^{1/4}\sqrt{\frac{-5\varphi_{4}+\sqrt{25\varphi_{4}^{2}-30\varphi_{2}f_{6,c}}}{f_{6,c}/2}},$
(63a) $\displaystyle\chi_{\text{eq}}^{\text{B-}}=$
$\displaystyle\pm\varepsilon^{1/4}\sqrt{\frac{-5\varphi_{4}-\sqrt{25\varphi_{4}^{2}-30\varphi_{2}f_{6,c}}}{f_{6,c}/2}}.$
(63b)

The zero-curvature phase corresponds to a local minimum of the free energy as
long as $f_{2}>0$, i.e. outside the bifurcation curve (22)—it becomes unstable
inside it, as already discussed. The curvatures
$\chi_{\text{eq}}^{\text{B}\pm}$ correspond to two buckled phases B$\pm$, the
domain of existence and stability of which is discussed below.

The buckled phase B+ exists as long as $f_{4,b}<0$ (or $\varphi_{4}<0$), i.e.
below the tricritical point, provided that $5f_{4,b}^{2}-6f_{2}f_{6,c}>0$ (or
$5\varphi_{4}^{2}-6\varphi_{2}f_{6,c}>0$). The line

$5f_{4,b}^{2}-6f_{2}f_{6,c}=0$ (64)

marks the end of the existence of this phase (rightmost dotted curve in Fig.
1). Moreover, when traversing from right to left this line, the buckled phase
pops up discontinuously with a finite curvature: specifically, from
$\chi_{\text{eq}}^{\text{ZC}}=0$ to
$\chi_{\text{eq}}^{\text{B+}}=\pm\sqrt{-5f_{4,b}/18}\neq 0$. It can be easily
checked—by evaluating the second derivative of the free energy density with
respect to $\chi$—that this buckled phase is locally stable in its domain of
existence.

The buckled phase B- also exists as long as $f_{4,b}<0$ (or
$\varphi_{4,b}<0$), i.e. below the tricritical point. However, not only does
it need $5f_{4}^{2}-6f_{2}f_{6,c}>0$ but also $f_{2}>0$. This phase B- only
exists between the line (64) and the part of the bifurcation curve (22) below
the tricritical point (dotted lines in Fig. 1). Phase B- is always unstable:
the second derivative of the free energy density for
$\chi=\chi_{\text{eq}}^{\text{B-}}$ is negative.

The curve demarcating the limit of existence of the buckled phases, Eq. (64),
can be written in a more transparent way. Let us consider, for any $T$ close
to $T_{c}$, a point $(J,T)$ close to the bifurcation curve $(J_{b}(T),T)$ by
taking into account that

$f_{4,b}\sim 6\sqrt{3}(T-T_{c}),\quad f_{2}\sim 2\sqrt{3}(J-J_{b}(T))$ (65)

in the approximation we are employing. Therefore, Eq. (64) is equivalent to

$J=J_{M}(T)\equiv J_{b}(T)+\frac{5}{12}\sqrt{3}(T-T_{c})^{2}.$ (66)

The analysis above entails the emergence of the region III in the plane
$(J,T)$, also depicted in Fig. 1. This region extends over the zone of the
plane $(J,T)$ below the tricritical point $C$, with its left border at the
bifurcation line (22) and its right border at the line (66) (both dotted lines
in Fig. 1). Inside region III, the three phases ZC, B+, and B- coexist. Phases
B+ and ZC correspond to local minima of the free energy, whereas B- always
corresponds to a local maximum. The relative stability of phases B+ and ZC is
elucidated in the following.

In region III, the phases ZC and B+ are both locally stable: one of them
corresponds to the deepest minima, being the most favourable thermodynamic
state, whereas the other one corresponds to a metastable state. This means
that the transition changes from second-order above the tricritical point to
first-order below it, since the curvature (order parameter) of the most stable
phase changes abruptly below the tricitical point. The change of stability
takes place at the first-order transition line determined by the condition
$\Delta\mathcal{F}(\text{B+})=\Delta\mathcal{F}(\text{ZC})=0$. Close to the
tricritical point, this is equivalent to

$f_{2}+\frac{1}{12}f_{4,b}(\chi_{\text{eq}}^{\text{B+}})^{2}+\frac{1}{360}f_{6,c}(\chi_{\text{eq}}^{\text{B+}})^{4}=0.$
(67)

Bringing to bear Eqs. (62) and (63a), over the first-order line one has
$8f_{2}f_{6,c}=5f_{4,b}^{2}$ that translates into

$J=J_{t}(T)\equiv J_{b}(T)+\frac{5}{16}\sqrt{3}(T-T_{c})^{2},$ (68)

making use of Eq. (65). The buckled phase B+ is the most stable one in region
IIIa of Fig. 1, i.e. between the branch of the bifurcation line below the
tricritical point (leftmost dotted line) and the first-order line (68) (dashed
line), with the ZC phase being metastable. The situation is just reversed in
region IIIb between the first-order line (68) and the curve demarcating the
limit of existence of phases B$\pm$, Eq. (66) (rightmost dotted line):
therein, the ZC phase is most stable and B+ is metastable.

##  Appendix B Continuum limit in honeycomb lattice

We go to a continuum limit, by assuming that $u_{ij}$ varies slowly with
$(i,j)$. We introduce continuous spatial variables $x$ and $y$ as depicted in
Fig. 5,

$\displaystyle x=$ $\displaystyle i\Delta x,$ $\displaystyle y=$
$\displaystyle j\Delta y,$ $\displaystyle\text{e-sites sublattice},$ (69a)
$\displaystyle x=$ $\displaystyle x_{o}+i\Delta x,$ $\displaystyle y=$
$\displaystyle j\Delta y,$ $\displaystyle\text{o-sites sublattice}.$ (69b)

where we have defined $x_{o}=a/2$, $\Delta x=3a/2$, and $\Delta
y=\sqrt{3}a/2$. All nearest neighbours of o-sites are e-sites, and vice
versa—the honeycomb lattice is bipartite [64]. Therefore, the term in
parentheses of the elastic term for e-sites is, neglecting $O(a^{3})$ terms,

$\displaystyle u_{i+1,j}+$ $\displaystyle u_{i,j-1}+u_{i,j+1}-3u_{ij}$
$\displaystyle=u(x+l,y)+u(x-l/2,y-a/2)$
$\displaystyle\quad+u(x-l/2,y+a/2)-3u(x,y)$
$\displaystyle\sim\left(\sqrt{3}\frac{a}{2}\right)^{2}\chi(x,y),\qquad\chi(x,y)\equiv\nabla^{2}u(x,y).$
(70)

where, in analogy with the one-dimensional case, we have introduced the
notation $\chi$ for the curvature of the membrane. The same result holds for
o-sites,
$u_{i-1,j}+u_{i,j-1}+u_{i,j+1}-3u_{ij}\sim\left(\sqrt{3}a/2\right)^{2}\chi(x,y)$.
Then, we get

$\displaystyle\frac{k}{2}(u_{i+1,j}+u_{i,j-1}+u_{i,j+1}-3u_{i,j})^{2}$
$\displaystyle=\frac{k_{0}}{2}\chi^{2},$ (71a) $\displaystyle
h(u_{i+1,j}+u_{i,j-1}+u_{i,j+1}-3u_{i,j})\sigma_{ij}$
$\displaystyle=h_{0}\chi\sigma_{ij},$ (71b)

in which we have defined $k_{0}=k(\sqrt{3}a/2)^{4}$ and
$h_{0}=h(\sqrt{3}a/2)^{2}$.

##  Appendix C Buckled profile on a rectangular domain

Here, we derive the analytical expression for the solution of Poisson’s
equation (46) on a rectangular domain of sides $L_{x}$ and $L_{y}$. We propose
the solution

$u(x,y)=\chi_{\text{eq}}\left[\frac{x(x-L_{x})}{2}+w(x,y)+z(x,y)\right].$ (72)

The first term on the rhs verifies Poisson’s equation and the homogeneous
Dirichlet boundary conditions at the sides $x=0$ and $x=L_{x}$, but not at the
sides $y=0$ and $y=L_{y}$. That is the reason why we introduce the functions
$w(x,y)$ and $z(x,y)$; both satisfy Laplace’s equation but with different
boundary conditions, specifically

$\displaystyle w(0,y)$ $\displaystyle=w(L_{x},y)=w(x,0)=0,$ (73a)
$\displaystyle w(x,L_{y})$ $\displaystyle=-x(x-L_{x})/2,$ (73b) $\displaystyle
z(0,y)$ $\displaystyle=z(L_{x},y)=z(x,L_{y})=0,$ (73c) $\displaystyle z(x,0)$
$\displaystyle=-x(x-L_{x})/2,$ (73d)

so that $u(x,0)=u(x,L_{y})=0$. Applying separation of variables and imposing
the boundary conditions, we obtain

$\displaystyle w(x,y)$
$\displaystyle=\frac{4L_{x}^{2}}{\pi^{3}}\sum_{n\scriptsize{\text{
odd}}}\frac{\sinh(\frac{n\pi y}{L_{x}})}{\sinh(\frac{n\pi
L_{y}}{L_{x}})}\frac{\sin\left(\frac{n\pi x}{L_{x}}\right)}{n^{3}},$ (74)
$\displaystyle z(x,y)$
$\displaystyle=\frac{4L_{x}^{2}}{\pi^{3}}\sum_{n\scriptsize{\text{
odd}}}\frac{\sinh(\frac{n\pi(L_{y}-y)}{L_{x}})}{\sinh(\frac{n\pi
L_{y}}{L_{x}})}\frac{\sin\left(\frac{n\pi x}{L_{x}}\right)}{n^{3}}.$ (75)

This result is presented in Fig. 7.

Figure 7: Absolute value of the stable buckled state on a square domain. The
graph corresponds to the low-temperature limit, in which $\chi_{\text{eq}}=1$,
for a unit area square with supported boundary conditions.

##  Appendix D Antiferromagnetic ground state for lattices without triangular
loops

The partition function of the Ising model with nearest-neighbour interactions
has only been analytically solved for the case of one-dimensional and two-
dimensional lattices [65]—in absence of external field for the latter.
Nevertheless, it is possible to derive an expression of the low-temperature
limit of the partition function for systems with specific properties, as shown
below.

For a general $d$-dimensional lattice, let the general antiferromagnetic Ising
Hamiltonian be

$\mathcal{H}(\bm{\sigma})=-h\sum_{i}\sigma_{i}+J\sum_{<i,j>}\sigma_{i}\sigma_{j},\qquad
h>0,J>0,$ (76)

where $\sigma_{i}$ stands for the two-state spin variable at site $i$ and sum
over $<i,j>$ means sum over all nearest-neighbour pairs. It is handy to
introduce the parameters defining the basic topological properties of the
lattice, which are the number of spins $N$ and the coordination number $c$,
the latter being the number of nearest neighbours of any site within the
lattice. The energy of the system can be fully characterised through the
number of spins aligned with the external field $n_{h}$, and the number of
nearest-neighbour pairs with anti-aligned spins $n_{J}$. Equivalently, we may
use the fractions $x_{h}\equiv n_{h}/N$ and $x_{J}=n_{J}/N$, which are
particularly adequate to analyse the large system size limit $N\to\infty$.
Taking into account this parametrisation, we have

$\displaystyle\mathcal{H}\left(x_{h},x_{J}\right)$
$\displaystyle=-h[n_{h}-(N-n_{h})]+\\!J\\!\left[\left(\dfrac{c}{2}N-n_{J}\right)-n_{J}\right]$
$\displaystyle=-N\left[h(2x_{h}-1)+J\left(2x_{J}-\dfrac{c}{2}\right)\right].$
(77)

In this way, a macrostate of the spin system is defined by a point
$\left(x_{h},x_{J}\right)$.

The partition function of the model is

$Z=\sum_{\hat{\sigma}}e^{-\beta\mathcal{H}(\bm{\sigma})}=\sum_{\left(x_{h},x_{J}\right)}\xi\left(x_{h},x_{J}\right)e^{-\beta\mathcal{H}\left(x_{h},x_{J}\right)},$
(78)

where we have introduced the multiplicity $\xi\left(x_{h},x_{J}\right)$ of the
macrostate—i.e. the number of microstates compatible with the values
$(x_{h},x_{J})$. In the low-temperature limit, $T\to 0^{+}$
($\beta\to+\infty$), the leading order of the partition function stems from
the ground state energy,

$-\frac{1}{N\beta}\ln Z\simeq e_{\text{gs}},\quad\beta\to\infty,$ (79)

where

$\displaystyle e_{\text{gs}}$
$\displaystyle\equiv\frac{1}{N}\min_{\left(x_{h},x_{J}\right)}\mathcal{H}\left(x_{h},x_{J}\right)$
$\displaystyle=\min_{\left(x_{h},x_{J}\right)}\left\\{-h\left(2x_{h}-1\right)-J\left(2x_{J}-\dfrac{c}{2}\right)\right\\}$
(80)

is the minimum energy of the system per lattice site. Therefore, the energy of
the ground state plays the role of the free energy of the system in the low-
temperature limit. Since Eq. (77) is a linear function of the two variables
$\left(x_{h},x_{J}\right)$, the minimum energy must be found at some point
belonging to the boundary of the physically available set
$\left(x_{h},x_{J}\right)$.

Note that $x_{h}$ is bounded both from below and from above by $0$ and $1$
respectively. Instead, $x_{J}$ is bounded by $0$ and $c/2$. Nevertheless, not
all points inside the rectangle
$\left(x_{h},x_{J}\right)\in[0,1]\times[0,c/2]$ are physically acceptable. The
specific shape of the physically available region in the
$\left(x_{h},x_{J}\right)$ space is not straightforward to derive for
arbitrary lattices. Nevertheless, let us restrict ourselves to the family of
lattices with no triangular loops that contains, for instance, all bipartite
networks—such as the honeycomb lattice in Fig. 5 or the rectangular lattice.
The assumption of the absence of triangular loops allows to characterise
exactly the available control set in the $\left(x_{h},x_{J}\right)$ space. If
this condition holds, the number of anti-aligned couples verifies

$x_{J}\leq c\ \text{min}(x_{h},1-x_{h}).$ (81)

Hence, the three vertices $(0,0)$, $(1/2,c/2)$, and $(1,0)$ define a
triangular region in the $\left(x_{h},x_{J}\right)$ plane in which all
physically possible macrostates are contained. Moreover, in the large system
size limit as $N\to\infty$, this region is densely filled by such
macrostates—both $x_{h}$ and $x_{J}$ become continuous variables in this
limit. The plane $(x_{h},x_{J})$ is shown in Fig. 8, the blue region
representing the physically available macrostates in the system—were we not
working in the limit $N\to\infty$, only some points inside the blue region
would represent acceptable macrostates.

Figure 8: Antiferromagnetic Ising model parameter space $(x_{h},x_{J})$. The
number of spins aligned with the external field is $n_{h}=Nx_{h}$, whereas the
number of antiferromagnetically aligned pairs in $n_{J}=Nx_{J}$. The filled
blue region indicates the physically available macrostates in the large system
size limit as $N\to\infty$, for an arbitrary $d$-dimensional lattice without
triangular loops.

After evaluating the energy at the boundaries, it is possible to obtain
finally what is the energy of the ground state—and thus the low-temperature
limit of the free energy of the system. Depending on the parameters $h$, $J$,
there appear only two possibilities:

  * •

All spins are aligned with the external field for $h>cJ$. In this situation,
the external field prevails over the ferromagnetic interaction. If so,
$n_{h}=N$ and $n_{J}=0$, i.e. $x_{h}=1$ and $x_{J}=0$.

  * •

The state is purely antiferromagnetic, i.e. all spins are antiparallel to
their nearest neighbours, for $cJ>h$. In this regime, it is the
antiferromagnetic coupling that dominates. Such a state exists due to our
assuming that triangular loops are absent—otherwise, it would not be possible
to have all nearest neighbours of any site $i$ antiparallel to $\sigma_{i}$,
some of them would be nearest neighbours among themselves. If so, $n_{h}=N/2$
and $n_{J}=cN/2$, i.e. $x_{h}=1/2$ and $x_{J}=c/2$.

Taking into account the above discussion, the energy per site in the ground
state can be cast in a single equation,

$e_{\text{gs}}=-\left(h-cJ\right)H(h-cJ)-\frac{c}{2}J,$ (82)

with $H$ denoting the Heaviside step function. This expression agrees with
Eqs. (28) and (53) in the main text for the one-dimensional and two-
dimensional cases, where (i) the absolute value of the curvature $|\chi|$
plays the role of the external field, and (ii) the coordination number is
$c=2$ for $d=1$ and $c=3$ on the honeycomb lattice for $d=2$. Still, it must
be remarked that Eq. (82) holds for an arbitrary $d$-dimensional lattice—as
long as it does not contain triangular loops.

## Acknowledgments

We acknowledge financial support from Grant No. PID2021-122588NB-I00 funded by
MCIN/AEI/10.13039/501100011033 and by “ERDF A way of making Europe.” C.A.P.
acknowledges the funding from European Union for the support from Horizon-
Marie Sklodowska-Curie 2021 programme through the Postdoctoral Fellowship Ref.
101065902 (ORION). We also thank Miguel Ruiz-García and Abigail Plummer for
useful discussions.

## References

  * Cea  _et al._ [2020] T. Cea, M. Ruiz-García, L. L. Bonilla, and F. Guinea, Numerical study of the rippling instability driven by electron-phonon coupling in graphene, Physical Review B 101, 235428 (2020).
  * Amnuanpol [2021] S. Amnuanpol, Buckling instability on a rotating plane, EPL (Europhysics Letters) 135, 50003 (2021).
  * Chen  _et al._ [2022] Z. Chen, D. Wan, and M. J. Bowick, Spontaneous Tilt of Single-Clamped Thermal Elastic Sheets, Physical Review Letters 128, 028006 (2022).
  * Hanakata  _et al._ [2021] P. Z. Hanakata, S. S. Bhabesh, M. J. Bowick, D. R. Nelson, and D. Yllanes, Thermal buckling and symmetry breaking in thin ribbons under compression, Extreme Mechanics Letters 44, 101270 (2021).
  * Hanakata  _et al._ [2022] P. Z. Hanakata, A. Plummer, and D. R. Nelson, Anomalous Thermal Expansion in Ising-like Puckered Sheets, Physical Review Letters 128, 075902 (2022).
  * Jain  _et al._ [2021] H. Jain, S. Ghosh, and K. C. Sahu, Compression-controlled dynamic buckling in thin soft sheets, Physical Review E 104, L033001 (2021).
  * Poincloux  _et al._ [2021] S. Poincloux, T. Chen, B. Audoly, and P. M. Reis, Bending Response of a Book with Internal Friction, Physical Review Letters 126, 218004 (2021).
  * Le Doussal and Radzihovsky [2021] P. Le Doussal and L. Radzihovsky, Thermal Buckling Transition of Crystalline Membranes in a Field, Physical Review Letters 127, 015702 (2021).
  * Thibado  _et al._ [2020] P. M. Thibado, P. Kumar, S. Singh, M. Ruiz-Garcia, A. Lasanta, and L. L. Bonilla, Fluctuation-induced current from freestanding graphene, Physical Review E 102, 042101 (2020).
  * Novoselov  _et al._ [2004] K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, Y. Zhang, S. V. Dubonos, I. V. Grigorieva, and A. A. Firsov, Electric Field Effect in Atomically Thin Carbon Films, Science 306, 666 (2004).
  * Neto  _et al._ [2009] A. C. Neto, F. Guinea, N. M. R. Peres, K. S. Novoselov, and A. K. Geim, The electronic properties of graphene, Reviews of modern physics 81, 109 (2009).
  * Samy  _et al._ [2021] O. Samy, S. Zeng, M. D. Birowosuto, and A. El Moutaouakil, A Review on MoS2 Properties, Synthesis, Sensing Applications and Challenges, Crystals 11, 355 (2021).
  * Huang  _et al._ [2019] H. Huang, S. Su, N. Wu, H. Wan, S. Wan, H. Bi, and L. Sun, Graphene-Based Sensors for Human Health Monitoring, Frontiers in Chemistry 7, 399 (2019).
  * Zhang  _et al._ [2022] R. Zhang, J. Jiang, and W. Wu, Scalably Nanomanufactured Atomically Thin Materials‐Based Wearable Health Sensors, Small Structures 3, 2100120 (2022).
  * Mangum  _et al._ [2021] J. M. Mangum, F. Harerimana, M. N. Gikunda, and P. M. Thibado, Mechanisms of Spontaneous Curvature Inversion in Compressed Graphene Ripples for Energy Harvesting Applications via Molecular Dynamics Simulations, Membranes 11, 516 (2021).
  * Gikunda  _et al._ [2022] M. N. Gikunda, F. Harerimana, J. M. Mangum, S. Rahman, J. P. Thompson, C. T. Harris, H. O. H. Churchill, and P. M. Thibado, Array of Graphene Variable Capacitors on 100 mm Silicon Wafers for Vibration-Based Applications, Membranes 12, 533 (2022).
  * Amorim  _et al._ [2016] B. Amorim, A. Cortijo, F. de Juan, A. Grushin, F. Guinea, A. Gutiérrez-Rubio, H. Ochoa, V. Parente, R. Roldán, P. San-Jose, J. Schiefele, M. Sturla, and M. Vozmediano, Novel effects of strains in graphene and other two dimensional materials, Physics Reports 617, 1 (2016).
  * Landau  _et al._ [1986] L. Landau, E. Lifshitz, A. Kosevich, J. Sykes, L. Pitaevskii, and W. Reid,  _Theory of Elasticity: Volume 7_ , Course of theoretical physics (Elsevier Science, 1986).
  * Golubović  _et al._ [1998] L. Golubović, D. Moldovan, and A. Peredera, Dynamics of the euler buckling instability, Phys. Rev. Lett. 81, 3387 (1998).
  * Singh  _et al._ [2015] S. K. Singh, M. Neek-Amal, S. Costamagna, and F. M. Peeters, Rippling, buckling, and melting of single- and multilayer MoS 2, Physical Review B 91, 014101 (2015).
  * Plummer and Nelson [2020] A. Plummer and D. R. Nelson, Buckling and metastability in membranes with dilation arrays, Physical Review E 102, 033002 (2020).
  * Shankar and Nelson [2021] S. Shankar and D. R. Nelson, Thermalized buckling of isotropically compressed thin sheets, Physical Review E 104, 054141 (2021).
  * Luo and Pozrikidis [2007] H. Luo and C. Pozrikidis, Buckling of a pre-compressed or pre-stretched membrane in shear flow, International Journal of Solids and Structures 44, 8074 (2007).
  * Hu  _et al._ [2017] W. Hu, M. Song, Z. Deng, T. Yin, and B. Wei, Axial dynamic buckling analysis of embedded single-walled carbon nanotube by complex structure-preserving method, Applied Mathematical Modelling 52, 15 (2017).
  * Meyer  _et al._ [2007] J. C. Meyer, A. K. Geim, M. I. Katsnelson, K. S. Novoselov, T. J. Booth, and S. Roth, The structure of suspended graphene sheets, Nature 446, 60 (2007).
  * Fasolino  _et al._ [2007] A. Fasolino, J. H. Los, and M. I. Katsnelson, Intrinsic ripples in graphene, Nature Materials 6, 858 (2007).
  * Note [1] It is important to recall that the Mermin-Wagner Theorem [66] forbids the existence of purely two-dimensional stable crystals. Also, the instability of two-dimensional crystals can be understood within the context of Landau’s theory of fluctuations [67].
  * San-Jose  _et al._ [2011] P. San-Jose, J. González, and F. Guinea, Electron-induced rippling in graphene, Physical Review Letters 106, 045502 (2011).
  * Bonilla and Carpio [2012] L. L. Bonilla and A. Carpio, Model of ripples in graphene, Physical Review B 86, 195402 (2012).
  * Bonilla  _et al._ [2012] L. L. Bonilla, A. Carpio, A. Prados, and R. R. Rosales, Ripples in a string coupled to Glauber spins, Physical Review E 85, 031125 (2012).
  * Ruiz-García  _et al._ [2016] M. Ruiz-García, L. L. Bonilla, and A. Prados, STM-driven transition from rippled to buckled graphene in a spin-membrane model, Physical Review B 94, 205404 (2016).
  * Ruiz-Garcia  _et al._ [2017] M. Ruiz-Garcia, L. L. Bonilla, and A. Prados, Bifurcation analysis and phase diagram of a spin-string model with buckled states, Physical Review E 96, 062147 (2017).
  * Wei  _et al._ [2013] Y. Wei, B. Wang, J. Wu, R. Yang, and M. L. Dunn, Bending Rigidity and Gaussian Bending Stiffness of Single-Layered Graphene, Nano Letters 13, 26 (2013).
  * Helfrich [1973] W. Helfrich, Elastic Properties of Lipid Bilayers: Theory and Possible Experiments, Zeitschrift für Naturforschung C 28, 693 (1973).
  * Lipowsky [1991] R. Lipowsky, The conformation of membranes, Nature 349, 475 (1991).
  * Schoelz  _et al._ [2015] J. K. Schoelz, P. Xu, V. Meunier, P. Kumar, M. Neek-Amal, P. M. Thibado, and F. M. Peeters, Graphene ripples as a realization of a two-dimensional Ising model: A scanning tunneling microscope study, Physical Review B 91, 045413 (2015).
  * Neek-Amal  _et al._ [2014] M. Neek-Amal, P. Xu, J. Schoelz, M. Ackerman, S. Barber, P. Thibado, A. Sadeghi, and F. Peeters, Thermal mirror buckling in freestanding graphene locally controlled by scanning tunnelling microscopy, Nature Communications 5, 4962 (2014).
  * Lindahl  _et al._ [2012] N. Lindahl, D. Midtvedt, J. Svensson, O. A. Nerushev, N. Lindvall, A. Isacsson, and E. E. B. Campbell, Determination of the Bending Rigidity of Graphene via Electrostatic Actuation of Buckled Membranes, Nano Letters 12, 3526 (2012).
  * Eder  _et al._ [2013] F. R. Eder, J. Kotakoski, K. Holzweber, C. Mangler, V. Skakalova, and J. C. Meyer, Probing from Both Sides: Reshaping the Graphene Landscape via Face-to-Face Dual-Probe Microscopy, Nano Letters 13, 1934 (2013).
  * Ruiz-García  _et al._ [2015] M. Ruiz-García, L. L. Bonilla, and A. Prados, Ripples in hexagonal lattices of atoms coupled to glauber spins, Journal of Statistical Mechanics: Theory and Experiment 2015, P05015 (2015).
  * Note [2] We use the term pseudospin to emphasise that it is not the real spin but a two-state variable that mimics in a simple way the range of values taken by the internal degrees of freedom.
  * Note [3] Yet, at some specific points in the paper we briefly consider the situation for ferromagnetic interaction, $J<0$.
  * Note [4] More specifically, the coupling between the transverse elastic modes and out-of-plane electrons, which can be thought to be in two states corresponding to the two different sides of the plane.
  * Note [5] Therefore, these terms are expected to be subdominant when $u_{j}$ varies slowly with $j$. The (small) corrections introduced thereby could be incorporated into our theoretical framework in a perturbative way.
  * Feynman [1996] R. P. Feynman,  _Statistical mechanics: a set of lectures_ , 19th ed., Frontiers in physics (Addison-Wesley, Reading, Mass., 1996).
  * Gelfand and Fomin [2000] I. M. Gelfand and S. V. Fomin,  _Calculus of Variations_ (Dover Publications, 2000).
  * Lanczos [1970] C. Lanczos,  _The variational principles of mechanics_ (Dover Publications, 1970).
  * Note [6] For instance, for free boundary conditions, the variation $\delta u(x)$ for $x\in(0,1)$ and $\\{\delta u(0),\delta u(1),\delta u^{\prime}(0),\delta u^{\prime}(1)\\}$ at the boundaries are independent.
  * Note [7] Zero curvature profiles, i.e. linear profiles with $\chi=0$, always solve Eq. (III).
  * Note [8] Note that $\beta h_{0}\chi=\beta^{*}\chi^{*}$ and $\beta J=\beta^{*}J^{*}$.
  * Note [9] The coupling constant is $J$, the factor $1/2$ is needed because each interaction is counted twice.
  * Barré  _et al._ [2005] J. Barré, F. Bouchet, T. Dauxois, and S. Ruffo, Large Deviation Techniques Applied to Systems with Long-Range Interactions, Journal of Statistical Physics 119, 677 (2005).
  * Winkler [2010] R. G. Winkler, Equivalence of statistical ensembles in stretching single flexible polymers, Soft Matter 6, 6183 (2010).
  * Manca  _et al._ [2014] F. Manca, S. Giordano, P. L. Palla, and F. Cleri, On the equivalence of thermodynamics ensembles for flexible polymer chains, Physica A: Statistical Mechanics and its Applications 395, 154 (2014).
  * Prados  _et al._ [2013] A. Prados, A. Carpio, and L. L. Bonilla, Sawtooth patterns in force-extension curves of biomolecules: An equilibrium-statistical-mechanics theory, Physical Review E 88, 012704 (2013).
  * Bonilla  _et al._ [2015] L. L. Bonilla, A. Carpio, and A. Prados, Theory of force-extension curves for modular proteins and DNA hairpins, Physical Review E 91, 052712 (2015).
  * Benedito and Giordano [2018] M. Benedito and S. Giordano, Isotensional and isometric force-extension response of chains with bistable units and Ising interactions, Physical Review E 98, 052146 (2018).
  * Fisher [1960] M. Fisher, Lattice statistics in a magnetic field, I. A two-dimensional super-exchange antiferromagnet, Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences 254, 66 (1960).
  * Peyrard [2006] M. Peyrard, Biophysics: Melting the double helix, Nature Physics 2, 13 (2006).
  * Prados  _et al._ [2012] A. Prados, A. Carpio, and L. L. Bonilla, Spin-oscillator model for the unzipping of biomolecules by mechanical force, Physical Review E 86, 021919 (2012).
  * Cannizzo  _et al._ [2022] A. Cannizzo, L. Bellino, G. Florio, G. Puglisi, and S. Giordano, Thermal control of nucleation and propagation transition stresses in discrete lattices with non-local interactions and non-convex energy, The European Physical Journal Plus 137, 569 (2022).
  * Bonilla [2000] L. L. Bonilla, Chapman-Enskog method and synchronization of globally coupled oscillators, Physical Review E 62, 4862 (2000).
  * Baldovin  _et al._ [2019] M. Baldovin, A. Vulpiani, A. Puglisi, and A. Prados, Derivation of a Langevin equation in a system with multiple scales: The case of negative temperatures, Physical Review E 99, 060101 (2019).
  * Newman [2010] M. Newman,  _Networks: An Introduction_ (OUP Oxford, 2010).
  * Onsager [1949] L. Onsager, Statistical hydrodynamics, Nuovo Cimento 6, 279 (1949).
  * Mermin and Wagner [1966] N. D. Mermin and H. Wagner, Absence of Ferromagnetism or Antiferromagnetism in One- or Two-Dimensional Isotropic Heisenberg Models, Physical Review Letters 17, 1133 (1966).
  * Landau and Lifshitz [2013] L. Landau and E. Lifshitz,  _Statistical Physics: Volume 5_ , Course of theoretical physics (Elsevier Science, 2013).