Buckling in a rotationally invariant spin-elastic model

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.

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 nonzero curvature (buckled) state, under load.A prototypical example occurs when a long column, i.e. a onedimensional 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 [27].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 lowdimensional 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 [41].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 [33][34][35].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 twodimensional 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 → 0 + , the transition is first-order and there is metastability: again, two buckled phases B± 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 nonessential 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, . . ., 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 σ j = ±1.We introduce the following rotation-ally invariant energy (or Hamiltonian) for any configuration of the elastic modes u, their associated momenta p, and the pseudospins σ, 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 [42].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 electronphonon coupling [43], and the antiferromagnetic interaction between the pseudospins account for the Coulomb repulsion of electrons.Therefore, nearest-neighbour outof-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 ≡ 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 = 0 corresponds to a small rotation of angle A, sin(A) A.
Thermal equilibrium of the system at temperature T is described by the canonical distribution P eq (u, p, σ) = exp [−βH(u, p, σ)] /Z, where β = (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 p and σ to get P eq (u) ∝ e −βF (u) , where the free energy of the string F(u) is (2) Above, Z int (u) stands for the partial partition function associated with the internal degrees of freedom, which involves a sum over all the pseudospin configurations Now we introduce a continuum limit by assuming that the displacements u j vary slowly with j, i.e., u j → u(x) with x = ja.The total length of the system is L = N a.In this way, where we have defined the curvature of the profile χ(x) ≡ u (x), and 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 pseudospinelastic interaction have been introduced as Also, terms with higher powers of a in Eq. ( 4) have been omitted, since they involve higher-order derivatives of u [44].
In the continuum limit, the equilibrium probability of a certain string profile u(x) becomes a functional thereof, specifically where n ≡ N/L is the number density, and the free energy density (per particle) f is given by ζ 1d int (βh 0 χ, βJ) ≡ e −βJ cosh (βh 0 χ) Our notation explicitly tells us that f only depends on the profile through its curvature χ.The second term in f (χ) comes from the sum of the spin variables, ζ 1d int is the one-particle partition function of a one-dimensional Ising chain with coupling J and external field h 0 χ-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 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 sizethus decoupling the continuum and the thermodynamic limits.

III. EQUILIBRIUM PROFILES
The equilibrium profiles u eq (x) are those that maximise the probability distribution P eq [u] or, equivalently, those that minimise the free energy functional F[u].Therefore, we consider the first variation of the free energy functional upon the change u → u + δu, The expression for δ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].[48] Considering the vanishing of the integral term, one obtains the Euler-Lagrange equation 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 but the value of u at the boundaries is free-the so-called supported boundary conditions [18].Therefore, ∂f /∂χ must vanish at both boundaries for the equilibrium profile, which entails A = B = 0.The equilibrium profile is thus provided by ∂f ∂χ eq = 0, i.e.
with the boundary conditions (11).Equation ( 12) is a transcendental equation for the equilibrium curvature χ eq , the solution of which can be written as χ eq = χ eq (T, J, k 0 , h 0 ).[49] The equilibrium curvature is constant, independent of position, and so is f eq = f (χ eq ).Therefore, the equilibrium free energy is F eq = nLf eq = N f eq , i.e. it is extensive-note that the constant χ eq is an intensive quantity, independent of system size.If χ eq is a solution of Eq. ( 12), then −χ eq is also a solution.This is a consequence of the free energy density being an even function of χ, 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 χ 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 ( 12) tells us the characteristic length 0 , which in turn determines the characteristic temperature T 0 , Then we define dimensionless displacement u * , position x * , curvature χ * , temperature T * , and coupling constant J * by and the dimensionless free energy F * ≡ F/k B T 0 .Then, we can write where f * is the dimensionless free energy density and ζ 1d int is the same one-particle partition function for the pseudospins than that in Eq. (8b).[50] Our use of dimensionless variables simplifies the theoretical analysis, by reducing the number of parametersand 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. ( 12) is rewritten as where µ * is the local magnetisation of the pseudospins.Equation ( 16) implies that χ * eq is a certain function of J * and T * , χ * eq = χ * 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, the equilibrium free energy is extensive.
For the supported boundary conditions (11) we are considering, there is a unique smooth equilibrium profile given by 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 * eq in Eq. ( 18), with C and D being arbitrary constants-the boundary condition d( ∂f * /∂χ * )/dx * | 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 χ ZC eq = 0 of Eq. ( 16).Second, we consider the low-temperature limit T 1, where the equilibrium profiles can be exactly calculated.Lastly, we analyse the phase diagram of the model by numerically solving Eq. ( 16).

A. Bifurcation from the zero-curvature solution
To start with, we expand the free energy density around the zero-curvature profile, for which i.e. we consider the difference of the free energy density from the zero-curvature value, The stability of the zero-curvature profile is determined by the sign of the second-derivative of the free energy density at χ = 0, The change of stability takes place at the line over which f 2 vanishes, i.e.
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 zerocurvature 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 zerocurvature value is provided by the Taylor expansion ∆f (χ; where t .On the one hand, the predictions for the bifurcation curve and for J (0) t , J (0) M are exact.On the other hand, the asymptotic theoretical predictions for the curves Jt(T ) and JM (T ) show an excellent agreement with the numerical results, within their range of validity.See Sec.IV C for further details in the numerical determination of the curves Jt(T ) and JM (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, and the asymptotic theoretical prediction for the curves J M (J t ) separating regions I and IIIb (IIIb and IIIa), Now we turn our attention to the low-temperature limit T 1.By defining e 1d gs as the energy per site in the ground state of the one-dimensional pseudospin system with external field χ and antiferromagnetic coupling J, we can write The ground state energy of the pseudospin system is in which H(x) is Heaviside's step function.Equation ( 28) is readily understood on a physical basis: for small |χ|, the antiferromagnetic coupling wins and the ground state corresponds to antiferromagnetic ordering, e 1d gs = −J, whereas for large |χ|, the external field wins and the ground state corresponds to all spins aligned with the curvature, e 1d gs = − |χ| + J.These two expressions for e 1d gs become equal at |χ| = 2J, which marks the crossover between them both.Of course, e 1d gs can also be derived by taking the limit T → 0 + in Eq. ( 27), bringing to bear Eq. (8b) for ζ 1d int .The above discussion entails that the free energy density (15b) simplifies to (29) Also, Eq. ( 16) for the equilibrium curvature converts to 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 χ eq → −χ eq : we recall that we are restricting ourselves to positive curvatures, so |χ eq | → χ 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 FIG. 2. Graphical method for solving Eq. (30).For each J, the intersection points of the functions H(χeq − 2J) (solid lines) and χeq (dashed blue line) give the equilibrium curvature in the low-temperature limit.For J > J (0) M = 1/2, for instance J = 0.55 (green solid line), the only solution is χeq = 0, i.e. the ZC phase.For J ≤ J (0) M , for instance J = 0.25 (red solid line), there appear three solutions: in addition to the ZC phase corresponding to χeq = 0, we have two buckled phases corresponding to χeq = 2J (B-phase) and χeq = 1 (B+ phase).
lhs and the rhs of Eq. ( 30) and look for the intersection points of both functions.It is clearly seen that χ 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 ≤ 1, i.e.J ≤ J (0) M = 1/2: specifically, χ eq = 2J and χ eq = 1.Both solutions coalesce at J (0) M : stronger antiferromagnetic interaction destroys the buckled phases, in which the magnetisation of the pseudospins is different from zero-a signature of ferromagnetic ordering.
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 ZC (J, T ) ∼ −J, for T 1, so that ∆f (χ; Therefore, one obtains ∆f (χ eq =2J; On the one hand, the phase with χ 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, ∆f (χ eq = 1; J, T → 0) changes sign at J = J (0) t = 1/4: the phase with χ eq = 1 has a free energy that is smaller (larger) than that of the ZC phase for 0 ≤ J < J M .Therefore, the solution with χ eq = 1 is expected to correspond to the low-temperature limit of phase B+, whereas the solution with χ eq = 2J should correspond to the low-temperature limit of phase B-.Also, the points J (0) t , 0 and J (0) M , 0 should correspond to the endpoints of the first-order line J t (T ) and the right border line of region III J M (T ), respectivelywe 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.

C. Numerical analysis
Here, we perform a numerical analysis of the phase diagram of our model system.We construct a 200×200 mesh in the (J, T ) plane and find χ 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 χ 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, χ eq changes continously from zero to non-zero values at the bifurcation curve J b (T ) (solid line).Below the tricritical point, χ 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, χ eq equals the magnetisation of the pseudospins, which is given by the rhs of Eq. ( 16).
By inserting the numerical solution for χ 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, ∆f (χ 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 J (0) M , 0 , and the curve J M (T ) that marks the limit of existence of the phases B± (rightmost dotted line).In the right panel, ∆f (χ 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 Bcoexist in region III, which is divided into two subregions by the first-order line J t (T ) (dashed line): IIIa, where ∆f (χ B+ eq ; J, T ) < 0 and the most stable phase is B+, be- ing ZC metastable; and IIIb, where the roles are reversed, the most stable phase is ZC, being B+ metastable.
Now we move to a two-dimensional system, i.e. a spinmembrane model.Specifically, we consider the system on a honeycomb lattice, such as that of graphene.There is a particle-carbon atom if thinking of grapheneat each node (i, j) in a certain two-dimensional region Ω.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 σ ij = ±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 Ω 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 δΩ of the region Ω, 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, For the sake of clarity, the Hamiltonian is split into two sums because nearest neighbours are different for esites 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.[51] 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 ∇ 2 u(x, y).Then, the sum over (i, j) goes to an integral over (x, y), which leads to where n = (∆x∆y) −1 is the number density and f is the free energy density per particle, given by ζ 2d int (βh 0 χ, βJ) stands for the partition function for the antiferromagnetic Ising chain with coupling J and external field h 0 χ in the two-dimensional case.Unfortunately, at variance with the one-dimensional case, there is not an analytical expression for ζ 2d int when J = 0.The equilibrium profile minimises the free energy, so it is determined by the condition δF [u] = 0. Since the free energy density only depends on the curvature χ, it is convenient to define the field Making use of this definition, the variation of the free energy can be written as In the contour integral along the boundary δΩ, ds is the length element and 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, i.e. the Laplace equation for the field Φ(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 where ∂/∂n ≡ n • ∇ 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, but the value of the normal derivative ∂ ∂n δu = 0 at the boundary is free.Therefore, the coefficient of ∂ ∂n δu in Eq. ( 39) has to vanish, Φ(x, y) = 0, (x, y) ∈ δΩ. ( The solution of the Laplace equation ( 38) with the homogeneous Dirichlet boundary condition ( 41) is identically zero: where we have made use of the definition of Φ.Although we do not have an explicit expression for ζ 2d int , Eq. ( 42) neatly tells us that the equilibrium curvature is constant, independent of (x, y) over all the two-dimensional domain Ω.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 2 0 .The dimensionless free energy and free energy density are where we have taken into account that βh In what follows, we drop again the asterisks in the dimensionless variables to simplify our notation.Equation ( 42) is written in dimensionless variables as χ eq = µ(βχ eq , βJ), µ(βχ eq , βJ) ≡ ∂ ln ζ 2d int (βχ eq , βJ) ∂(βχ eq ) .
(44) Note that µ is the local magnetisation of the twodimensional pseudospin lattice.The equilibrium free energy F eq is given by F eq (J, T ) ≡ F[u eq ](J, T ) = N f (χ eq ; J, T ), (45) which is extensive, as expected on a physical basis.

A. 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 ũ(x, y) = u(x, y) + Ax + By + C, is another possible equilibrium profile with the same free energy, which only depends on the curvature ∇ 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: ∇ 2 u eq (x, y) = χ eq , (x, y) ∈ Ω; u(x, y) = 0, (x, y) ∈ δΩ.
(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 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 ζ 2d int can be worked out.Specifically, two situations are considered: (i) the case (J = 0, T = 0), in which the antiferromagnetic coupling among the pseudospins disappears, and (ii) the lowtemperature limit (J = 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 (0) M of the antiferromagnetic coupling.Moreover, metastability is present, signaling that the transition from zero curvature to buckled states is now first-order.

No antiferromagnetic coupling
Let us first consider the case J = 0.The one-particle partition ζ 2d int (J = 0, T ) corresponds to that of a non-interacting Ising system in an external field χ, i.e.
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 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 →0 + χ eq (J = 0, T ) = 1. (51)

Low-temperature limit
Now we turn to the low-temperature limit T → 0 + , with J = 0.In this regime, similarly to the situation in the one-dimensional spin-string system, we have lim where is the energy per site in the ground state of the twodimensional spin-membrane.Again, Eq. ( 53) for e 2d gs has a clear physical interpretation.On the one hand, the ground state energy corresponds to antiferromagnetic ordering for small χ, e 2d gs = −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 χ, e 2d gs = −χ + 3J/2.The crossing between these expressions for e 2d gs takes place at χ = 3J, which leads to Eq. ( 53)-see details in Appendix D.
The free energy density is now and the equilibrium curvature is fixed by We recall that, without loss of generality, we are restricting ourselves to positive curvatures in the paper, thus |χ eq | → χ eq henceforth.The analysis of the lowtemperature 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 ≤ 1, i.e.J ≤ J (0) M = 1/3.The two buckled phases have curvatures χ eq = 3J and χ eq = 1.For J > J (0) M , 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 ZC (J, ( and The phase with χ eq = 3J always has ∆f > 0. The phase with χ eq = 1 has ∆f < 0 (> 0) for 0 M .Then, the phase with χ eq = 3J should correspond to the (unstable) B-phase, whereas the phase with χ eq = 1 should correspond to the (at least locally stable) B+ phase.The buckled phases disappear discontinuously at J = J (0) M , which means that the transition is first-order.
The above identification of the phases is physically sensible, since in the limit J → 0 the B+ phase for T 1 is expected to converge to the T → 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.

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 onedimensional and the two-dimensional lattices.
FIG. 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 ≤ J (0) M = 1/3, corresponding to the most stable state for 0 ≤ J < J (0) t = 1/6 and to a metastable state for 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 (0) M .For J ≤ J (0) M , 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 ≤ J < J (0) t , with the situation being reversed for M .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 onedimensional case.Nevertheless, our analysis of limiting cases allows us to state that the same qualitative picture of phases applies for the two-dimensional spinmembrane model.Still, the two-dimensional values of J (0) t and J (0) M are smaller-in dimensionless variablesthan 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 onedimensional 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 [57,[59][60][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 (χ = 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: ∂∆f (χ;J,T ) This expansion is truncated at different orders, the sign of the last retained coefficient f n must be always positiveotherwise, the corresponding equilibrium distribution P eq would be non-normalisable.
a. f4 > 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 < ε 1 such that . This tells us that the separation of the point (J, T ) we are considering from the bifurcation curve is of order ε: either 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 ≤ 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 ϕ 2 > 0, the free energy density has only one minimum at χ ZC eq = 0. On the other hand, if ϕ 2 < 0, χ ZC eq = 0 turns out to be a maximum and there appear two symmetric minima at where f 4,b ≡ f 4 (J b , T b ) is the value of the coefficient f 4 over the bifurcation curve, The superindex B in Eq. (A2) 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. (A1), f 6 χ 5 eq /6! is in fact negligible: it is of order ε 5/2 , whereas the first two terms on the rhs are of order ε 3/2 -this shows the consistency of the approximations introduced.
b. f4 < 0 and f6 > 0: Tricritical point and first-order phase transition The coefficient f 4,b over the bifurcation curve vanishes at the tricritical point C ≡ (J c , T c ), given by Close to the tricritical point, the approximation given by Eq. (A2) 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. (A1) is needed.We still write f 2 = εϕ 2 : the terms involving f 2 and f 6 are of the same order if χ eq = O(ε 1/4 ).Then, the terms involving f The partition function of the Ising model with nearestneighbour 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 lowtemperature limit of the partition function for systems with specific properties, as shown below.σ i σ j , h > 0, J > 0, (D1) where σ i stands for the two-state spin variable at site i and sum over < i, j > means sum over all nearestneighbour 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 ≡ n h /N and x J = n J /N , which are particularly adequate to analyse the large system size limit N → ∞.Taking into account this parametrisation, we have In this way, a macrostate of the spin system is defined by a point (x h , x J ).The partition function of the model is ξ (x h , x J ) e −βH(x h ,x J ) , (D3) where we have introduced the multiplicity ξ (x h , x J ) of the macrostate-i.e. the number of microstates compatible with the values (x h , x J ).In the low-temperature limit, T → 0 + (β → +∞), the leading order of the partition function stems from the ground state energy, 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. (D2) is a linear function of the two variables (x h , x J ), the minimum energy must be found at some point belonging to the boundary of the physically available set (x h , x J ).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 (x h , x J ) ∈ [0, 1] × [0, c/2] are physically acceptable.The specific shape of the physically available region in the (x h , x J ) 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 (x h , x J ) space.If this condition holds, the number of anti-aligned couples verifies (D6) Hence, the three vertices (0, 0), (1/2, c/2), and (1, 0) define a triangular region in the (x h , x J ) plane in which all physically possible macrostates are contained.Moreover, in the large system size limit as N → ∞, 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 systemwere we not working in the limit N → ∞, only some points inside the blue region would represent acceptable macrostates.
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 absentotherwise, it would not be possible to have all nearest neighbours of any site i antiparallel to σ i , some of them would be nearest neighbours among themselves.If so, n h = N/2 and n J = cN/2, i.e.
Taking into account the above discussion, the energy per site in the ground state can be cast in a single equation,

f 6 (
FIG.1.Phase diagram of the rotationally invariant spinstring.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 = (Jc, Tc) 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 > Tc, it defines the already mentioned second-order phase transition (black solid line).For T < Tc, it demarcates the change from region II to region III (leftmost dotted line).The curve JM (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 Jt(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 (0) M and J (0)

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

FIG. 4 .
FIG.4.Free energy density difference ∆f (χ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.
2 and f 4,b are of the same order if f 4,b = O(ε 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

FIG. 7 .
FIG.7.Absolute value of the stable buckled state on a square domain.The graph corresponds to the low-temperature limit, in which χeq = 1, for a unit area square with supported boundary conditions.

FIG. 8 .
FIG. 8. Antiferromagnetic Ising model parameter space (x h , xJ ).The number of spins aligned with the external field is n h = N x h , whereas the number of antiferromagnetically aligned pairs in nJ = N xJ .The filled blue region indicates the physically available macrostates in the large system size limit as N → ∞, for an arbitrary d-dimensional lattice without triangular loops.

e
gs = − (h − cJ) H(h − cJ) − c 2 J,(D7)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 |χ| 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. (D7) holds for an arbitrary d-dimensional lattice-as long as it does not contain triangular loops.