Buckling in a Rotationally Invariant Spin-elastic Model

PHYSICAL REVIEW E 107, 014120 (2023) Buckling in a rotationally invariant spin-elastic model Gregorio García-Valladares ,*Carlos A. Plata ,†and Antonio Prados ‡ Física Teórica, Universidad de Sevilla, Apartado de Correos 1065, E-41080 Sevilla, Spain (Received 27 September 2022; accepted 22 December 2022; published 13 January 2023) Scanning tunneling microscopy experiments have revealed a 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 analyze 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. Minimization of this free energy gives rise to the emergence of different mechanical phases, whose thermodynamic stability is thoroughly analyzed, both analytically and numerically. All phases are characterized 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. DOI: 10.1103/PhysRevE.107.014120 I. INTRODUCTION Behavior of low-dimensional systems is a prolific research field [1–9]. Two-dimensional materials exhibit intriguing properties [10–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 realization 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 nonzero 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]. Not only is buckling an engineering subject but also an appealing physical behavior 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–22] to cells [23]orindevice development [24]. Graphene sheets are not perfectly two-dimensional. In fact, they display transversal displacements (ripples) due to thermal fluctuations [25,26] that stabilize the membrane.1Rippling *ggv[email protected] †[email protected] ‡[email protected] 1It is important to recall that the Mermin-Wagner Theorem [27] 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 [28]. has been studied employing models where its origin and behavior are linked with the coupling between electronic and elastic degrees of freedom, i.e., the electron-phonon coupling [1,26,29–33]. Rippling has also been investigated by employing phenomenological models such as the Helfrich model [34], originally developed for membranes [35,36]. 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 [37–40], 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 2470-0045/2023/107(1)/014120(17) 014120-1 ©2023 American Physical Society GARCÍA-VALLADARES, PLATA, AND PRADOS PHYSICAL REVIEW E 107, 014120 (2023) cases, respectively, are able to reproduce in a qualitative way the aforementioned transition. Despite their simplicity, these mesoscopic models present a rich behavior, with a complicated phase diagram that include buckled phases [32,33,41]. Nevertheless, the currently available spin-string and spinmembrane models [32,33,41], 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 behavior 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. [32,33,41] and mends the above-mentioned fundamental flaws. We consider a classical elastic lattice (one or two-dimensional) characterized by the out-of-plane deformations and internal degrees of freedom described by a pseudospin variable.2 Elastic deformations are modeled by a harmonic interaction between the atoms that depends on the curvature of the membrane, the repulsive interaction between electrons is modeled by an antiferromagnetic interaction between pseudospins, and the electron-phonon coupling is modeled 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 characterized by a free energy that depends on the membrane curvature—in close analogy with the phenomenological approach based on the Helfrich model [34–36]. The equilibrium value of the curvature only depends on two parameters: temperature Tand 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 characterize 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 makes 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 T0. As the antiferromagnetic 2We use the term pseudospin to emphasize 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. coupling is increased, the second-order transition from rippled to buckled moves to lower temperatures and at a certain critical point (Jc,Tc) the transition changes to first order. Beyond this tricritical point, basically for low enough temperatures in a given range of J’s, the behavior 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 spinstring system. On the one hand, for J=0, there appears a second-order transition from a rippled phase (ZC) to a buckled phase (B). On the other hand, 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 onedimensional situation, the relative stability of the B+and ZC phases depends on the intensity of the antiferromagnetic coupling. The paper is organized as follows. In Sec. II, we define the one-dimensional model and study its equilibrium behavior, 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 analyze, both analytically and numerically, the emergence of buckled states and put forward the main results of our study. In Sec. V, we generalize 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 mthat is characterized by its transversal displacement uj, and, in addition, there is a “pseudospin” variable σj=±1. We introduce the following rotationally invariant energy (or Hamiltonian) for any configuration of the elastic modes u, 014120-2 BUCKLING IN A ROTATIONALLY INVARIANT … PHYSICAL REVIEW E 107, 014120 (2023) their associated momenta p, and the pseudospins σ, H(u,p,σ)= N  j=0p2 j 2m+k 2(uj+1−2uj+uj−1)2 −h(uj+1−2uj+uj−1)σj+Jσj+1σj.(1) The first term on the right-hand side (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 htunes the strength thereof, and the last one represents the nearest-neighbor interaction among the pseudospins—with coupling constant J. Our Hamiltonian (1) can be analyzed 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 caseJ>0.3The 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,4and the antiferromagnetic interaction between the pseudospins account for the Coulomb repulsion of electrons. Therefore, nearestneighbor 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 [32,33,41]: (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 [32,33,41]. In particular, a flat profile, uj≡0, and a rigid rotation thereof, uj=Aj with constant A,havethesame energy in Eq. (1). In general, any two profiles differing just by a linear function Aj+Bhas 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 Peq(u,p,σ)= exp[−βH(u,p,σ)]/Z, where β=(kBT)−1, with kBand Z being the Boltzmann constant and the partition function, respectively. We are interested in the equilibrium profiles of the string, so we integrate over pand σto get Peq(u)∝e−βF(u), 3Yet, at some specific points in the paper we briefly consider the situation for ferromagnetic interaction, J<0. 4More 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. where the free energy of the string F(u)is F(u)= N  j=0k 2(uj+1−2uj+uj−1)2−kBTln Zint(u).(2) Above, Zint(u) stands for the partial partition function associated with the internal degrees of freedom, which involves a sum over all the pseudospin configurations Zint = σ e−βN j=0[−h(uj+1−2uj+uj−1)σj+Jσj+1σj].(3) Now we introduce a continuum limit by assuming that the displacements ujvary slowly with j, i.e., uj→u(x) with x= ja. The total length of the system is L=Na. In this way, uj+1−2uj+uj−1=a2χ(x)+O(a4),(4) where we have defined the curvature of the profile χ(x)≡ u(x), and k 2(uj+1−2uj+uj−1)2=k0 2χ2,(5a) h(uj+1−2uj+uj−1)σj=h0χσj.(5b) The continuum limit in the spin variables is skipped willingly, due to the marginalization over all spin configurations carried out in Eqs. (2) and (3). Scaled parameters for the elastic constant and the strength of the pseudospin-elastic interaction have been introduced as k0=ka4,h0=ha2.(6) Also, terms with higher powers of ain Eq. (4) have been omitted, since they involve higher-order derivatives of u.5 In the continuum limit, the equilibrium probability of a certain string profile u(x) becomes a functional thereof, specifically Peq[u]∝e−βF[u],F[u]=nL 0 dxf(χ),(7) where n≡N/Lis the number density, and the free-energy density (per particle) fis given by f(χ)=k0χ2 2−β−1ln ζ1D int (βh0χ,βJ),(8a) ζ1D int (βh0χ,βJ)≡e−βJcosh (βh0χ) +eβJ1+e−4βJsinh2(βh0χ).(8b) Our notation explicitly tells us that fonly depends on the profile through its curvature χ. The second term in f(χ) comes from the sum of the spin variables, ζ1D int is the oneparticle partition function of a one-dimensional Ising chain with coupling Jand external field h0χ—e.g., see Ref. [42]. In contrast to previous spin-string and spin-membrane models [32,33,41], the continuum limit consider here does not involve a large system size limit N1 (thermodynamic 5Therefore, these terms are expected to be subdominant when uj varies slowly with j. The (small) corrections introduced thereby could be incorporated into our theoretical framework in a perturbative way. 014120-3 GARCÍA-VALLADARES, PLATA, AND PRADOS PHYSICAL REVIEW E 107, 014120 (2023) 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 hand kdo not scale with the system size—thus decoupling the continuum and the thermodynamic limits. III. EQUILIBRIUM PROFILES The equilibrium profiles ueq(x) are those that maximize the probability distribution Peq[u] or, equivalently, those that minimize the free-energy functional F[u]. Therefore, we consider the first variation of the free-energy functional upon the change u→u+δu, δF[u]=L 0 dxδud2 dx2∂f ∂χ +δu∂f ∂χ −δud dx∂f ∂χ L 0 . (9) 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 [43,44].6Considering the vanishing of the integral term, one obtains the Euler-Lagrange equation ∂2 ∂x2∂f ∂χ eq =0⇒∂f ∂χ eq =Ax +B,(10) where Aand Bare 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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., k0χeq =kBT∂ln ζ1D int ∂χ χ=χeq =h0 e−2J kBTsinh h0χeq kBT 1+e−4J kBTsinh2h0χeq kBT ,(12) 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,k0,h0).7 The equilibrium curvature is constant, independent of position, and so is feq =f(χeq). Therefore, the equilibrium free energy is Feq =nL feq =Nf eq, i.e., it is extensive—note that 6For instance, for free boundary conditions, the variation δu(x)for x∈(0,L)and{δu(0),δu(L),δu(0),δu(L)}at the boundaries are independent. 7Zero curvature profiles, i.e., linear profiles with χ=0, always solve Eq. (12). 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 xaxis 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 characterize the equilibrium curvature χeq,itis 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 T0, 0=k0 h0=k ha2,T0=h2 0 kBk0=h2 kBk.(13) Then we define dimensionless displacement u∗, position x∗, curvature χ∗, temperature T∗, and coupling constant J∗by u∗=u 0 ,x∗=x 0 ,χ ∗=0χ, T∗=T T0 ,J∗=J kBT0 . (14) Consistently, we define β∗=kBT0β=1/T∗,L∗=L/0, n∗=n0=N/L∗, and the dimensionless free energy F∗≡ F/kBT0. Then, we can write F∗[u∗](J∗,T∗)=n∗L∗ 0 dx∗f∗(χ∗;J∗,T∗),(15a) f∗(χ∗;J∗,T∗)≡1 2χ∗2−T∗ln ζ1D int (β∗χ∗,β∗J∗), (15b) 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).8 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-like transition 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 (∂f∗/∂χ∗)eq =0, i.e., χ∗ eq =μ∗(β∗χ∗ eq,β∗J∗),(16a) μ∗(β∗χ∗,β∗J∗)≡∂ln ζ1D int ∂(β∗χ∗)=e−2β∗J∗sinh (β∗χ∗) 1+e−4β∗J∗sinh2(β∗χ∗), (16b) 8Note that βh0χ=β∗χ∗and βJ=β∗J∗. 014120-4 BUCKLING IN A ROTATIONALLY INVARIANT … PHYSICAL REVIEW E 107, 014120 (2023) where μ∗is the local magnetization 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 nonzero 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 T0, 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, F∗ eq(J∗,T∗)≡F∗[u∗ eq](J∗,T∗)=Nf∗(χ∗ 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∗ eq(x∗;J∗,T∗)=χ∗ eq(J∗,T∗) 2x∗(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∗+Dto u∗ eq in Eq. (18), with C and Dbeing 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∗+Dfor the free case. In the remainder of the paper, we employ dimensionless variables—dropping the asterisks in them not to clutter our formulas. IV. BUCKLED STATES In this section, we analyze the emergence of buckled states, i.e., with nonzero curvature, in our model system. First, we show how buckled profiles bifurcate from the zero-curvature solution χZC eq =0ofEq.(16). Second, we consider the lowtemperature limit T1, where the equilibrium profiles can be exactly calculated. Lastly, we analyze 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 fZC(J,T)≡f(χ=0; J,T)=−Tln [2 cosh(J/T)],(19) i.e., we consider the difference of the free-energy density from the zero-curvature value, f(χ;J,T)≡f(χ;J,T)−fZC(J,T) =1 2χ2−Tln ζ1D int (χ/T,J/T) 2 cosh(J/T).(20) The stability of the zero-curvature profile is determined by the sign of the second-derivative of the free-energy density at χ=0, f2(J,T)≡∂2f ∂χ2χ=0=1−e−2J/T T.(21) The change of stability takes place at the line over which f2 vanishes, i.e., J=Jb(T)=−1 2Tln T,(22) which is called the bifurcation curve. This line separates the plane (J,T) into two regions: (i) a region with f2>0, comprising two subregions I and III where the zero-curvature solution is stable and (ii) region II with f2<0, where the zero-curvature solution becomes unstable and other stable solution exist. Note that for J=0, we have that f2<0for T<1: the unit of temperature T0is 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(χ;J,T)=f2(J,T) 2χ2+f4(J,T) 4! χ4 +f6(J,T) 6! χ6+O(χ8),(23) where f4(J,T)=e−6J/T T3(3 −e4J/T),(24a) f6(J,T)=−e−10J/T T5(e8J/T−30e4J/T+45).(24b) 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 summarizes the main results of the equilibrium and stability analyses carried out in detail in Appendix A. Therefrom, we highlight here the tricritical point, Tc=1 √3,Jc=Jb(Tc)=ln 3 4√3,(25) and the asymptotic theoretical prediction for the curves JM(Jt) separating regions I and IIIb (IIIb and IIIa), JM(T)=Jb(T)+5 12√3(T−Tc)2,(26a) Jt(T)=Jb(T)+5 16√3(T−Tc)2.(26b) 014120-5 GARCÍA-VALLADARES, PLATA, AND PRADOS PHYSICAL REVIEW E 107, 014120 (2023) FIG. 1. Phase diagram of the rotationally invariant spin-string model. The plane (J,T) is divided into several regions, characterized 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 Jb(T), givenbyEq.(22), is a bifurcation curve. For T>Tc,itdefines 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) Mand J(0) 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)andJM(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)andJM(T). B. Low-temperature limit Now we turn our attention to the low-temperature limit T1. By defining e1D 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 lim T→0+−Tln ζ1D int =e1D gs .(27) The ground-state energy of the pseudospin system is e1D gs =− (|χ|−2J)H(|χ|−2J)−J,(28) 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 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 lowtemperature 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). to antiferromagnetic ordering, e1D gs =−J, whereas for large |χ|, the external field wins and the ground state corresponds to all spins aligned with the curvature, e1D gs =−|χ|+J. These two expressions for e1D gs become equal at |χ|=2J, which marks the crossover between them both. Of course, e1D 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 f(χ;J,T)∼χ2 2−(|χ|−2J)H(|χ|−2J)−J,T1. (29) Also, Eq. (16) for the equilibrium curvature converts to |χeq|=H(|χeq|−2J),T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 χeq →−χeq: we recall that we are restricting ourselves to positive curvatures, so |χeq|→χeq in the following. To study the solutions of Eq. (30), a graphical method is useful. In Fig. 2, we plot the functions on the left-hand side (lhs) and the rhs of Eq. (30) and look for the intersection points of both functions. It is clearly seen that χeq =0is 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 = 2Jand χeq =1. Both solutions coalesce at J(0) M: stronger antiferromagnetic interaction destroys the buckled phases, in which the magnetization of the pseudospins is different from zero—a signature of ferromagnetic ordering. 014120-6 BUCKLING IN A ROTATIONALLY INVARIANT … PHYSICAL REVIEW E 107, 014120 (2023) FIG. 3. Curvature of the most stable phase (numerically obtained). The labeling 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 Cto first order (dashed line) below, is clearly observed. Let us calculate the free-energy density of each of the phases, to identify them. First, the free-energy density of the ZC phase is given by fZC(J,T)∼−J,forT1, so that f(χ;J,T)=f(χ;J,T)−fZC(J,T) ∼χ2 2−(|χ|−2J)H(|χ|−2J),T1. (31) Therefore, one obtains f(χeq =2J;J,T→0) ∼2J2,(32a) f(χeq =1; J,T→0) ∼−1 2+2J.(32b) On the one hand, the phase with χeq =2Jhas 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(0) t(J(0) t<J⩽J(0) M). Therefore, the solution with χeq =1 is expected to correspond to the low-temperature limit of phase B+, whereas the solution with χeq =2Jshould 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 Jt(T) and the right border line of region III JM(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. 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 3shows 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 continuously from zero to nonzero values at the bifurcation curve Jb(T) (solid line). Below the tricritical point, χeq changes discontinuously from zero to nonzero values at the first-order line Jt(T) (dashed line). Recall that the curvature plays the role of the order parameter in our model. In fact, χeq equals the magnetization 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 FIG. 4. Free-energy density difference f(χeq;J,T) of the buckled states with respect to the ZC phase (numerically obtained). The labeling of the different regions and the code of the different lines is the same as in Fig. 1and 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. 014120-7 GARCÍA-VALLADARES, PLATA, AND PRADOS PHYSICAL REVIEW E 107, 014120 (2023) FIG. 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 neighbors below and one above, and o-sites (red), two nearest neighbors above and one below. The distance between sites of the same kind corresponding to consecutive indexes are x=3a/2, y=√3a/2. The spacing between e- and o-sites with the same value of the index iis xo=a/2. 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 JM(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 B−coexist in region III, which is divided into two subregions by the first-order line Jt(T) (dashed line): IIIa, where f(χB+ eq ;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 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 graphene—at each node (i,j) in a certain two-dimensional region . The transverse displacement with respect to the plane of each particle is denoted by uij. Analogously to the one-dimensional spin-string model analyzed 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 iand jare 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 neighbors are one above and two below; for the latter, the three nearest neighbors are two above and one below [41]. 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 neighbors 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, H(ˆu,ˆp,ˆσ)= |i−j|=even p2 ij 2m+k 2(ui−1,j+ui,j−1+ui,j+1−3uij)2−hσij(ui−1,j+ui,j−1+ui,j+1−3uij) +J 2σij(σi−1,j+σi,j−1+σi,j+1) + |i−j|=odd p2 ij 2m+k 2(ui+1,j+ui,j−1+ui,j+1−3uij)2−hσij(ui+1,j+ui,j−1+ui,j+1−3uij) +J 2σij(σi+1,j+σi,j−1+σi,j+1).(33) For the sake of clarity, the Hamiltonian is split into two sums because nearest neighbors 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, 014120-8 BUCKLING IN A ROTATIONALLY INVARIANT … PHYSICAL REVIEW E 107, 014120 (2023) and the fourth one provides the antiferromagnetic coupling among the pseudospins.9Again, we are interested in the equilibrium profiles of the membrane. Therefore, we integrate the canonical distribution over the momenta and the pseudospin variables, first 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 χ≡∇2u(x,y). Then, the sum over (i,j) goes to an integral over (x,y), which leads to Peq[u]∝e−βF[u],F[u]=n dxdy f(χ),(34) where n=(xy)−1is the number density and fis the freeenergy density per particle, given by f(χ)=k0χ2 2−kBTln ζ2D int (βh0χ,βJ),(35a) ζ2D int (βh0χ,βJ) stands for the partition function for the antiferromagnetic Ising chain with coupling Jand external field h0χ 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 minimizes 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 (x,y)≡∂f ∂χ =k0χ−h0 ∂ln ζ2D int (βh0χ,βJ) ∂(βh0χ).(36) Making use of this definition, the variation of the free energy can be written as δF[u]= dxdy∇2(δu) = dxdy∇·[∇(δu)−∇δu]+ dxdy∇2δu =δ ds ˆ n·(∇(δu)−∇δu)+ dxdy∇2δu. (37) In the contour integral along the boundary δ,ds is the length element and ˆ nis 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, ∇2=0,(38) 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 ∂ ∂nδu=0,δu∂ ∂n=0,(39) 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 9The coupling constant is J, the factor 1/2 is needed because each interaction is counted twice. boundary conditions: the value of the transversal displacements vanish at the boundaries, u(x,y)=0,(x,y)∈δ, (40) but the value of the normal derivative ∂ ∂nδu=0 at the boundary is free. Therefore, the coefficient of ∂ ∂nδuin Eq. (39) has to vanish, (x,y)=0,(x,y)∈δ. (41) The solution of the Laplace equation (38) with the homogeneous Dirichlet boundary condition (41) is identically zero: k0χeq −h0 ∂ln ζ2D int (βh0χeq,βJ) ∂(βh0χeq)=0,(42) 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, with 0=k0/h0.The dimensionless free energy and free-energy density are F∗[u∗](J∗,T∗)=n∗ dx∗dy∗f∗(χ∗;J∗,T∗),(43a) f∗(χ∗;J∗,T∗)≡1 2χ∗2−T∗ln ζ2D int (β∗χ∗,β∗J∗),(43b) where we have taken into account that βh0χ=β∗χ∗,βJ= β∗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 χeq =μ(βχeq,βJ),μ(βχeq,βJ)≡∂ln ζ2D int (βχeq,βJ) ∂(βχeq). (44) Note that μis the local magnetization of the two-dimensional pseudospin lattice. The equilibrium free energy Feq is given by Feq(J,T)≡F[ueq](J,T)=Nf(χ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 spinstring model, the two-dimensional spin-membrane model analyzed 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 ˜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 ∇2u(x,y). Similar to the one-dimensional case, the family of transformations defined by Ax +By +Ccontains 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 014120-9 GARCÍA-VALLADARES, PLATA, AND PRADOS PHYSICAL REVIEW E 107, 014120 (2023) [11] A. H. CastroNeto, F. Guinea, N. M. R. Peres, K. S. Novoselov, and A. K. Geim, The electronic properties of graphene, Rev. Mod. Phys. 81, 109 (2009). [12] 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). [13] H. Huang, S. Su, N. Wu, H. Wan, S. Wan, H. Bi, and L. Sun, Graphene-based sensors for human health monitoring, Front. Chem. 7, 399 (2019). [14] R. Zhang, J. Jiang, and W. Wu, Scalably nanomanufactured atomically thin materials-based wearable health sensors, Small Struct. 3, 2100120 (2022). [15] 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). [16] 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). [17] 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, Phys. Rep. 617, 1 (2016). [18] L. Landau, E. Lifshitz, A. Kosevich, J. Sykes, L. Pitaevskii, and W. Reid, Theory of Elasticity: Volume 7, Course of theoretical physics (Elsevier Science, Amsterdam, 1986). [19] L. Golubovic, D. Moldovan, and A. Peredera, Dynamics of the Euler Buckling Instability, Phys. Rev. Lett. 81, 3387 (1998). [20] S. K. Singh, M. Neek-Amal, S. Costamagna, and F. M. Peeters, Rippling, buckling, and melting of single- and multilayer MoS 2, Phys. Rev. B 91, 014101 (2015). [21] A. Plummer and D. R. Nelson, Buckling and metastability in membranes with dilation arrays, Phys. Rev. E 102, 033002 (2020). [22] S. Shankar and D. R. Nelson, Thermalized buckling of isotropically compressed thin sheets, Phys. Rev. E 104, 054141 (2021). [23] H. Luo and C. Pozrikidis, Buckling of a pre-compressed or prestretched membrane in shear flow, Int. J. Solids Struct. 44, 8074 (2007). [24] 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, Appl. Math. Modell. 52, 15 (2017). [25] 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 (London) 446, 60 (2007). [26] A. Fasolino, J. H. Los, and M. I. Katsnelson, Intrinsic ripples in graphene, Nat. Mater. 6, 858 (2007). [27] N. D. Mermin and H. Wagner, Absence of Ferromagnetism or Antiferromagnetism in One- Or Two-Dimensional Isotropic Heisenberg Models, Phys. Rev. Lett. 17, 1133 (1966). [28] L. Landau and E. Lifshitz, Statistical Physics: Volume 5, Course of theoretical physics (Elsevier Science, Amsterdam, 2013). [29] P. San-Jose, J. González, and F. Guinea, Electron-Induced Rippling in Graphene, Phys.Rev.Lett.106, 045502 (2011). [30] L. L. Bonilla and A. Carpio, Model of ripples in graphene, Phys. Rev. B 86, 195402 (2012). [31] L. L. Bonilla, A. Carpio, A. Prados, and R. R. Rosales, Ripples in a string coupled to glauber spins, Phys. Rev. E 85, 031125 (2012). [32] M. Ruiz-García, L. L. Bonilla, and A. Prados, STM-driven transition from rippled to buckled graphene in a spin-membrane model, Phys.Rev.B94, 205404 (2016). [33] M. Ruiz-Garcia, L. L. Bonilla, and A. Prados, Bifurcation analysis and phase diagram of a spin-string model with buckled states, Phys.Rev.E96, 062147 (2017). [34] Y. Wei, B. Wang, J. Wu, R. Yang, and M. L. Dunn, Bending rigidity and Gaussian bending stiffness of single-layered graphene, Nano Lett. 13, 26 (2013). [35] W. Helfrich, Elastic properties of lipid bilayers: Theory and possible experiments, Z. Naturforsch., C 28, 693 (1973). [36] R. Lipowsky, The conformation of membranes, Nature (London) 349, 475 (1991). [37] 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, Phys.Rev.B91, 045413 (2015). [38] 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, Nat. Commun. 5, 4962 (2014). [39] 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 Lett. 12, 3526 (2012). [40] 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 Lett. 13, 1934 (2013). [41] M. Ruiz-García, L. L. Bonilla, and A. Prados, Ripples in hexagonal lattices of atoms coupled to Glauber spins, J. Stat. Mech.: Theory Exp. (2015) P05015. [42] R. P. Feynman, Statistical Mechanics: A Set of Lectures, 19th ed., Frontiers in Physics (Addison-Wesley, Reading, MA, 1996). [43] I. M. Gelfand and S. V. Fomin, Calculus of Variations (Dover Publications, Mineola, NY, 2000). [44] C. Lanczos, The Variational Principles of Mechanics (Dover Publications, Mineola, NY, 1970). [45] J. Barré, F. Bouchet, T. Dauxois, and S. Ruffo, Large deviation techniques applied to systems with long-range interactions, J. Stat. Phys. 119, 677 (2005). [46] R. G. Winkler, Equivalence of statistical ensembles in stretching single flexible polymers, Soft Matter 6, 6183 (2010). [47] F. Manca, S. Giordano, P. L. Palla, and F. Cleri, On the equivalence of thermodynamics ensembles for flexible polymer chains, Physica A 395, 154 (2014). [48] A. Prados, A. Carpio, and L. L. Bonilla, Sawtooth patterns in force-extension curves of biomolecules: An equilibriumstatistical-mechanics theory, Phys.Rev.E88, 012704 (2013). 014120-16 BUCKLING IN A ROTATIONALLY INVARIANT … PHYSICAL REVIEW E 107, 014120 (2023) [49] L. L. Bonilla, A. Carpio, and A. Prados, Theory of forceextension curves for modular proteins and DNA hairpins, Phys. Rev. E 91, 052712 (2015). [50] M. Benedito and S. Giordano, Isotensional and isometric forceextension response of chains with bistable units and Ising interactions, Phys. Rev. E 98, 052146 (2018). [51] M. Fisher, Lattice statistics in a magnetic field, I. A twodimensional super-exchange antiferromagnet, Proc.R.Soc. London A 254, 66 (1960). [52] M. Peyrard, Biophysics: Melting the double helix, Nat. Phys. 2, 13 (2006). [53] A. Prados, A. Carpio, and L. L. Bonilla, Spin-oscillator model for the unzipping of biomolecules by mechanical force, Phys. Rev. E 86, 021919 (2012). [54] A. Cannizzo, L. Bellino, G. Florio, G. Puglisi, and S. Giordano, Thermal control of nucleation and propagation transition stresses in discrete lattices with nonlocal interactions and nonconvex energy, Eur. Phys. J. Plus 137, 569 (2022). [55] L. L. Bonilla, Chapman-Enskog method and synchronization of globally coupled oscillators, Phys. Rev. E 62, 4862 (2000). [56] 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, Phys. Rev. E 99, 060101(R) (2019). [57] M. Newman, Networks: An Introduction (Oxford University Press, Oxford, UK, 2010). [58] L. Onsager, Statistical hydrodynamics, Nuovo Cimento 6, 279 (1949). 014120-17