# Thermalized buckling of isotropically compressed thin sheets

Suraj Shankar1 and David R. Nelson1,2 1Department of Physics, Harvard
University, Cambridge, MA 02138, USA
2Department of Molecular and Cellular Biology and School of Engineering and
Applied Sciences, Harvard University, Cambridge, Massachusetts 02138, USA

###### Abstract

The buckling of thin elastic sheets is a classic mechanical instability that
occurs over a wide range of scales. In the extreme limit of atomically thin
membranes like graphene, thermal fluctuations can dramatically modify such
mechanical instabilities. We investigate here the delicate interplay of
boundary conditions, nonlinear mechanics and thermal fluctuations in
controlling buckling of confined thin sheets under isotropic compression. We
identify two inequivalent mechanical ensembles based on the boundaries at
constant strain (isometric) or at constant stress (isotensional) conditions.
Remarkably, in the isometric ensemble, boundary conditions induce a novel
long-ranged nonlinear interaction between the local tilt of the surface at
distant points. This interaction combined with a spontaneously generated
thermal tension leads to a renormalization group description of two distinct
universality classes for thermalized buckling, realizing a mechanical variant
of Fisher-renormalized critical exponents. We formulate a complete scaling
theory of buckling as an unusual phase transition with a size dependent
critical point, and discuss experimental ramifications for the mechanical
manipulation of ultrathin nanomaterials.

##  1 Introduction

Thin sheets with a resistance to shear can accommodate compressive stresses
through an array of mechanical instabilities, including buckling Landau and
Lifshitz (1986); Koiter (1967), wrinkling Cerda  _et al._ (2002); *Cerda2003;
Davidovitch  _et al._ (2011); Vandeparre  _et al._ (2011), folding Pocivavsek
_et al._ (2008); Diamant and Witten (2011) and crumpling Ben Amar and Pomeau
(1997); *Cerda1998; *Cerda1999; Witten (2007), all controlled essentially by
geometry. Although once disregarded as undesirable modes of failure,
instabilities now play a central role in the design of mechanical
metamaterials Krieger (2012); Bertoldi  _et al._ (2017) as they combine
complex morphologies with mechanical functionality. In recent years, rapid
miniaturization has driven intense research efforts in developing similar
metamaterials on a much smaller scale Blees  _et al._ (2015); Rogers  _et al._
(2016); *xu2017ultrathin; *miskin2018graphene; *reynolds2019capillary;
*grosso2020graphene. In this regard, atomically thin two dimensional (2D)
materials such as graphene, MoS2 or BN Novoselov  _et al._ (2005); Katsnelson
(2012) are particularly promising and offer unprecedented opportunities to
study classical elasticity and mechanics in the ultimate limit in thin sheets,
where thermal fluctuations can play a dominant role Nelson  _et al._ (2004);
Katsnelson (2012).

In such ultrathin flexible materials thermal fluctuations can dramatically
renormalize the mechanical properties in a scale-dependent fashion Nelson and
Peliti (1987). Out of plane (flexural) deformations allow tensionless solid
membranes to exhibit a remarkable thermally wrinkled, yet flat phase with a
scale-dependent bending rigidity and strongly softened elastic moduli Bowick
and Travesset (2001); Nelson  _et al._ (2004). While thin sheets favour
bending over energetically expensive stretching, geometry links the two as any
bending-induced Gaussian curvature inevitably causes stretching as well. This
basic feature underlies many of the impressive finite temperature properties.
Nanoindentation measurements in graphene Lee  _et al._ (2008) and MoS2
Bertolazzi  _et al._ (2011) monolayers yield exceptionally high Young’s moduli
on the nanoscale as expected from strong covalent bonding. Yet on larger
scales $\sim 10~{}\mu$m, recent experiments with freely suspended graphene
have demonstrated a $\sim 4000$ fold enhancement of the bending rigidity Blees
_et al._ (2015) and a factor $\sim 20$ reduction in the in-plane stiffness
Nicholl  _et al._ (2015); *nicholl2017hidden, due to a combination of
thermally generated and static ripples Meyer  _et al._ (2007); Košmrlj and
Nelson (2013); Xu  _et al._ (2014), highlighting the importance of flexural
fluctuations.

While the anomalous mechanics of thermalized membranes has been extensively
explored, the role of confinement and boundaries is much less appreciated.
Supported or clamped edges are one of the most commonly encountered boundary
conditions, in electromechanical resonators Bunch  _et al._ (2007); Chen  _et
al._ (2009), multistable switches Loh and Espinosa (2012) and in
nanomechanical devices Blees  _et al._ (2015); Nicholl  _et al._ (2015, 2017).
Geometric confinement at the boundary can induce prestrains in the sample that
can cause large scale instabilities such as wrinkling Bao  _et al._ (2009). As
a result, in recent years, exploring the influence of external stresses on the
mechanics of fluctuating membranes has been a topic of prime interest Amorim
_et al._ (2016). While there has been some theoretical work, both old Guitter
_et al._ (1988); *guitter1989thermodynamical and new Roldán  _et al._ (2011);
Košmrlj and Nelson (2016); Bonilla and Ruiz-Garcia (2016); Burmistrov  _et
al._ (2018a); *burmistrov2018differential; *saykin2020absolute, complemented
by more recent large scale numerical simulations Jiang (2014); Bowick  _et
al._ (2017); Wan  _et al._ (2017); Sgouros  _et al._ (2018); Morshedifard  _et
al._ (2021); Hanakata  _et al._ (2021), elucidating the role of boundaries in
controlling the nonlinear mechanics and buckling of thermalized sheets,
particularly for compressions which attempt to impose a nonzero Gaussian
curvature, remains a challenging problem.

Motivated by the above, in this paper, we pose and answer the following
question - what is the finite temperature version of the buckling transition
in an isotropically compressed thin sheet? Euler buckling represents the
simplest mechanical instability a thin elastic body can undergo and it
provides an attractive setting to investigate the interplay of thermal
fluctuations and boundary conditions along with the geometric nonlinearities
inherent to thin plate mechanics. In particular, we focus on the universal
aspects of the transition such as critical scaling exponents that are
independent of microscopic details. By combining a full renormalization group
analysis along with a general scaling theory, we provide a complete
description of thermalized buckling as a genuine phase transition that
exhibits critical scaling along with more unusual features such as a sensitive
dependence on system size and the choice of boundary conditions. Apart from
new statistical mechanical results, the scaling framework we propose also
yields key predictions that have important ramifications for experiments that
we highlight below.

In the rest of the introduction, we summarize our main results and outline the
structure of the paper.

###  1.1 Results and outline

A key outcome of our work is a renormalization group analysis, augmented by a
scaling description of buckling in isotropically compressed thermalized thin
sheets. There are two main reasons why finite temperature buckling, when
viewed as a phase transition, is distinct from conventional critical
phenomena. The first is that buckling is strongly size dependent, even at zero
temperature by virtue of it being a long-wavelength instability Landau and
Lifshitz (1986). The second is the remarkable fact that freely fluctuating
elastic sheets exhibit a flat phase Nelson and Peliti (1987) with critical
fluctuations over an extended range of temperatures. Both these features are
characteristic of thin sheets, arising from an interplay of geometry and
mechanics, and form the basis of our results below.

###  1 _Ensemble inequivalence & Fisher renormalization_

Figure 1: A sketch showing a possible realization of the two mechanical
ensembles, for example, in an atomically thin sheet of graphene. (a) In the
isotensional ensemble, a constant inward external stress $\sigma_{0}$ is
applied to the membrane, while the boundary displacement fluctuates. This set-
up might be realized in the same fashion as in single molecule experiments, by
using feedback-controlled multiplexed optical tweezers to actuate under
constant force conditions, similar to experiments recently used to probe the
mechanical response of red blood cells Turlier  _et al._ (2016). (b) The
isometric ensemble instead corresponds to a clamped boundary with the external
load imposed via a global strain $\epsilon$. While, current experiments with
graphene typically suspend monolayers across fixed size holes Nicholl  _et
al._ (2015), a variable aperture size tuned by a camera shutter mechanism
could be used to tune the strain isotropically. An external symmetry breaking
field $\mathcal{E}$ perpendicular to average plane of the sheet can also be
applied in either ensemble to bias the direction of buckling.

Thin sheets can be loaded in-plane either by prescribing the external strain
(isometric) or the external stress (isotensional) at the boundary (see Fig. 1
for an illustration). These constitute dual mechanical ensembles, in analogy
with thermodynamic ensembles Gibbs (1902). While it is well known that
boundary conditions can modify nonuniversal quantities such as the buckling
threshold Landau and Lifshitz (1986), we discover that universal scaling
exponents (defined in Sec. 4) can also exhibit a remarkable sensitivity to
boundary conditions! We demonstrate this fact explicitly within a systematic
$\varepsilon=4-D$ expansion for a general $D$-dimensional solid embedded in
$d>D$ dimensional space (discussed in Secs. 7 and 8 with details relegated to
Appendix C) along with a simpler, but uncontrolled, one-loop calculation
performed directly in the physical dimensions of $D=2$ and $d=3$ in Appendix
D. Our calculations show that buckling in the two mechanical ensembles is in
fact controlled by two distinct fixed points, with different scaling exponents
that are summarized in Table 1. This remarkable departure from conventional
wisdom demonstrates the nonequivalence of mechanical ensembles in thermally
fluctuating thin sheets and highlights the subtle nature of membrane
statistical mechanics. Thus, in the simplest setting of isotropic compression,
we find thermalized plate buckling is described by two distinct critical
points characterizing the isometric and isotensional universality classes that
are distinguished simply by the boundary conditions imposed.

Although surprising, the inequivalence of ensembles has precedence in critical
phenomena. A critical point engenders fluctuations on all scales which under
appropriate conditions can result in scaling exponents that change upon
switching to a dual or constrained ensemble. This is known as Fisher
renormalization Fisher (1968); *lipa1968critical; *fisher1970visibility. In
our case, however, by tuning to the buckling threshold, we approach the flat
phase of a free-standing membrane which characterizes an entire low
temperature _critical phase_ with scale invariant fluctuations! Furthermore,
we find that the imposition of a fixed strain (in contrast to a fixed stress)
boundary condition induces a novel long-ranged nonlinearity that couples the
local _tilt_ of the surface at far away points, which we derive in Sec. 3 and
Appendix A through a careful consideration of zero modes and boundary
conditions. This nonlocal term, which can also be important far from the
buckling transition, softly enforces the geometric constraint of global
inextensibilitiy, which simultaneously shifts the buckling threshold by a
spontaneously generated thermal tension and modifies the critical exponents
via a mechanical variant of Fisher renormalization.

###  2 _Size-dependent scaling theory_

The long-wavelength nature of the buckling instability endows it with both a
system size dependent threshold and a macroscopic mechanical response Landau
and Lifshitz (1986), features that are retained even at finite temperature.
This size dependence is unusual though, from the point of view of critical
phenomena, and behaves as a _dangerously irrelevant_ variable Amit and Peliti
(1982); *Gunton1973RenormalizationGI; *nelson1976coexistence that modifies
scaling exponents in nontrivial ways. As a result, we derive new exponent
identites in Sec. 9 that relate different scaling exponents in both ensembles,
mirroring classic results from conventional critical phenomena Goldenfeld
(2018). Many of these relations are also summarized in Table 1. By combining
scaling with general thermodynamic arguments, we also explicitly demonstrate
how buckling physics in both ensembles is a mechanical variant of Fisher
renormalization. Note that, the construction of a consistent scaling theory
for thermalized buckling is a significant achievement as it not only clarifies
previous confounding results Guitter  _et al._ (1988, 1989); Košmrlj and
Nelson (2016) by correctly accounting for nontrivial system size dependence
and ensemble inequivalence, but it also yields a unified framework that
incorporates experimentally relevant boundary conditions and symmetry breaking
fields.

###  3 _Experimental consequences_

Our work illustrates the spectacular ways in which geometry, boundary effects
and thermal fluctuations can conspire to produce unexpected phenomena, and
suggests that extending thin body mechanics to finite temperature is a rich
and challenging enterprise, requiring great care. Given the ubiquity and ease
of manipulating strain rather than stress in experiments, our results have
important implications for the rational design of strain engineered
nanodevices and interpretation of mechanical measurements in ultrathin
materials in the presence of thermal fluctuations.

Figure 2: Experimentally measured height response curve of a clamped graphene
sheet suspended over a circular hole of radius $6.2~{}\mu$m in a cryostat at
two different temperatures, $T=78~{}$K (blue dots) and $T=297~{}$K (red dots).
The data is reproduced from Ref. Storch  _et al._ (2018). The electrically
integrated graphene devices are capacitively actuated out of the plane using
an electrostatic force $\mathcal{E}\propto V_{g}^{2}$ ($V_{g}$ is the gate
voltage) and the average deflection $\langle h\rangle$ of the sheet is
measured using laser interferometry, see inset for a sketch of the setup and
Ref. Storch  _et al._ (2018) for further details. The red and blue lines are
guides to the eye showing the exponent of the nonlinear force response. While
the $T=297~{}$K data (red dots) shows a strong cubic dependence on the height,
the $T=78~{}$K data (blue dots) shows a smaller slope that matches well with
our theoretical prediction using the _isometric_ ensemble exponent
$3-1/\beta^{\prime}\approx 1.607$. Note that, this slope is significantly
different from a slope of unity (black line), as would be predicted in both
the _isotensional_ ensemble and within mean field theory.

A common setup to probe the mechanical properties of graphene involves
measuring the force-displacement curve for a clamped monolayers that is
deformed by the application of stresses or external fields. As a simple
example, we shall focus here on a simple circular geometry with sheets draped
across holes and deflected by a normal electric field ($\mathcal{E}$) as
employed in recent experiments Nicholl  _et al._ (2015); Storch  _et al._
(2018) (see inset in Fig. 2 for a sketch). Other geometries including ribbons
cantilevered at an edge Blees  _et al._ (2015) or suspended across trenches
Bunch  _et al._ (2007); Chen  _et al._ (2009); Bao  _et al._ (2009) are also
possible, but we don’t address them here. A crucial ingredient in the
interpretation of force response measurements in these devices is a mechanical
“equation of state” that relates the externally applied field $\mathcal{E}$
and the in-plane tension $\Delta\sigma$ to the magnitude of the average
deflection of the sheet $\langle h\rangle$. A zero-temperature mean-field
description following classical elasticity (detailed in Sec. 5) gives

$\mathcal{E}=c_{1}\dfrac{\Delta\sigma}{R^{2}}\langle
h\rangle+c_{2}\dfrac{Y}{R^{4}}\langle h\rangle^{3}\;,$ (1.1)

where $c_{1,2}$ are calculable numerical constants, $Y$ is the sought after 2D
Young’s modulus and $R$ is the size of the sheet. Importantly, when the sheet
is constrained at fixed stress $\Delta\sigma$, i.e., the isotensional
ensemble, Eq. 1.1 applies exactly even at finite temperature upon simply
replacing $Y$ by its scale dependent renormalized value (see Sec. 9 for
details).

However, nearly all the force measurement experiments are instead conducted
with fixed strain and clamped boundary conditions, i.e., in the _isometric_
ensemble. Our general scaling theory and renormalization group analysis
provides a different equation

$\mathcal{E}=c_{1}\dfrac{\Delta\sigma}{R^{2}}\left(\dfrac{Y}{R^{2}}\right)^{1-1/2\beta^{\prime}}\langle
h\rangle^{3-1/\beta^{\prime}}+c_{2}\dfrac{Y}{R^{4}}\langle h\rangle^{3}\;,$
(1.2)

where the tension $\Delta\sigma=B\Delta\epsilon$ is now given by the 2D bulk
modulus $B$ and the imposed in-plane strain $\Delta\epsilon$. Strikingly, Eq.
1.2 involves a new order parameter exponent $\beta^{\prime}\approx 0.718$ that
controls the asymptotic nonlinear force response. The difference between Eq.
1.1 and Eq. 1.2 makes it clear that using Eq. 1.1 for clamped sheets, as is
conventionally done, can lead to significantly wrong results. In fact, recent
experimental measurements Storch  _et al._ (2018) on graphene drumheads match
well with our theoretical predictions, see Fig. 2. The strain in the sample
increases upon lowering temperature Storch  _et al._ (2018); Chen  _et al._
(2009), notwithstanding the theoretically expected negative thermal expansion
coefficient of graphene Košmrlj and Nelson (2016); Bao  _et al._ (2009),
presumably due to surface contaminants. As a result, while the classical cubic
response dominates at higher temperature with weak tension (red dots in Fig.
2), the anomalous nonlinear response with $\mathcal{E}\propto\langle
h\rangle^{1.607}$ (blue dots in Fig. 2) emerges for higher strains at lower
temperature. A systematic analysis exploring how this nonlinear response
affects the extraction of the Young’s modulus is left for future work.

This result highlights the direct relevance of our work to the correct
interpretation of mechanical measurements in graphene devices, not only for
the circular geometries studied here, but also cantilevers and doubly clamped
ribbons. We believe that recognizing the fundamental distinction between the
isotensional and the isometric ensembles is essential to such endeavours.
While Fig. 2 depicts a static example, a dynamical extension of Eq. 1.2
including dissipation and inertia along with a time-varying field
$\mathcal{E}(t)$ provides a simple description of periodicaly driven
electromechanical oscillators Lifshitz and Cross (2008). Although a full
dynamical analysis is beyond the scope of this current work, we can already
appreciate the presence of a strong nonlinear response (Eq. 1.2) in the small
deflection limit, which allows for higher quality factors and bistability
Lifshitz and Cross (2008) even for a weak drive. Note that, such an anomalous
response is only elicited in the isometric ensemble, emphasizing once again
the importance of boundary conditions.

Although we focus on isotropic compression, we expect suitable extensions of
this framework to be applicable to recent numerical simulations Morshedifard
_et al._ (2021); Hanakata  _et al._ (2021) and experiments Blees  _et al._
(2015); Bao  _et al._ (2009); Xu  _et al._ (2014) on compressed ribbons that
have begun addressing the anomalous mechanics of _anisotropic_ buckling in
constrained sheets. While much of our discussion has been based on atomic
crystals, we emphasize our results are generic and also apply to organic 2D
polymers, such as naturally occuring DNA kinetoplasts Klotz  _et al._ (2020),
the spectrin cytoskeleton Schmidt  _et al._ (1993) and reconstituted spider
silk Hermanson  _et al._ (2007), or synthetic molecular crystals Servalli  _et
al._ (2018); *payamyar2016two and possibly polymerosomes Shum  _et al._ (2008)
with very large radii. In Sec. 10, we conclude with a brief discussion of the
broader significance of our results to experiments on atomically thin
crystalline membranes and possible future directions.

##  2 Thin plate elasticity

Here, we focus on the physically relevant continuum elastic description of a
thin 2D sheet fluctuating in 3D space (a generalization for a $D$-dimensional
solid embedded in $d$-dimensional ($d>D$) Euclidean space, useful for certain
calculations, is provided in Appendix C). The deformation of a thin flat sheet
is parametrized in the Monge gauge using a 2D in-plane displacement vector
${\bf u}$ and a height field $h$. The total elastic energy of the sheet
involves both stretching and bending contributions and is given by Landau and
Lifshitz (1986)

$\mathcal{H}=\int\mathrm{d}^{2}r\left[\dfrac{\kappa}{2}(\nabla^{2}h)^{2}+\mu
u_{ij}^{2}+\dfrac{\lambda}{2}u_{kk}^{2}-\mathcal{E}h\right]-\oint_{\mathcal{C}}\mathrm{d}\ell\;\hat{\nu}_{i}\sigma^{\rm
ext}_{ij}u_{j}\;.$ (2.1)

The Lamé parameters are $\mu$ and $\lambda$, and $\kappa$ is the bending
rigidity. The final boundary integral is the work done by an external stress
$\bm{\sigma}^{\rm{ext}}$ with $\hat{\bm{\nu}}$ being the outward normal
(within the plane) to the boundary curve $\mathcal{C}$. The penultimate term
corresponds to the potential energy due to an external out of plane field
$\mathcal{E}$ which couples directly to the height of the membrane. Such a
perturbation can be realized by an electric field $E$, with
$\mathcal{E}=\rho_{q}E$, where $\rho_{q}$ is the electric charge density on
the surface, while in the presence of gravity, we have
$\mathcal{E}=\rho_{m}g$, where $\rho_{m}$ is the mass density and $g$ is the
gravitational acceleration. The strain tensor

$u_{ij}=\dfrac{1}{2}\left(\partial_{i}u_{j}+\partial_{j}u_{i}+\partial_{i}h\partial_{j}h\right)\;,$
(2.2)

encodes the geometric nonlinearity inherent to thin sheets. We neglect higher
order terms in the in-plane displacements which are small and irrelevant on
large scales for a thin sheet Landau and Lifshitz (1986). As an aside, note
that, for a $D$-dimensional manifold embedded in $d$-dimensional space, the
displacement field ${\bf u}$ has $D$ components and the height field is no
longer a scalar, but instead a vector ${\bf h}$ with codimension $d_{c}=d-D>0$
components. The nonlinear strain tensor in this case is then
$u_{ij}=(\partial_{i}u_{j}+\partial_{j}u_{i}+\partial_{i}{\bf
h}\cdot\partial_{j}{\bf h})/2$. The relative importance of stretching versus
bending energies is captured by a dimensionless Föppl-von Kármán number, which
in 2D is given by

${\rm vK}=\dfrac{YR^{2}}{\kappa}\;,$ (2.3)

where $Y=4\mu(\mu+\lambda)/(2\mu+\lambda)$ is the 2D Young’s modulus and $R$
is a characteristic linear dimension of the sheet. When viewing the sheet as
thin slice of a bulk elastic material, its bending modulus and stiffness are
related as $\kappa=Y_{\rm 3D}t^{3}/[12(1-\nu_{\rm 3D}^{2})]$ and $Y=Y_{\rm
3D}t$, where $Y_{\rm 3D}$ is the 3D Young’s modulus, $\nu_{\rm 3D}$ the 3D
Poisson’s ratio and $t$ is the thickness of the sheet Audoly and Pomeau
(2010). As a result, ${\rm vK}\sim(R/t)^{2}$ is essentially controlled by
geometry with ${\rm vK}\gg 1$ for a thin sheet ($t/R\ll 1$), reflecting the
dominance of geometrically nonlinear isometric deformations, i.e., bending
without stretching. An ordinary sheet of paper has a ${\rm vK}\approx 10^{6}$
while a $1~{}\mu$m size graphene monolayer has a microscopic ${\rm vK}\approx
10^{9}$ (using the atomic scale values for $\kappa\approx 1.1~{}$eV Katsnelson
(2012) and $Y\approx 340~{}$N/m Lee  _et al._ (2008)).

##  3 Mechanical ensembles

The properties of a thermalized elatic membrane at temperature $T$ are
computed through the equilibrium partition function,

$\mathcal{Z}=\int\mathcal{D}h\mathcal{D}{\bf u}\;e^{-\mathcal{H}/k_{B}T}\;,$
(3.1)

where $k_{B}$ is the Boltzmann constant. As the in-plane phonons (${\bf u}$)
only appear quadratically in $\mathcal{H}$ they can be integrated out exactly
to give an effective free energy $\mathcal{F}=-k_{B}T\ln\int\mathcal{D}{\bf
u}\;e^{-\mathcal{H}/k_{B}T}$. To do this, we separate out the average strain
and Fourier transform the nonzero wavelength deformations. While the
calculation for the wavevector ${\bf q}\neq{\bf 0}$ modes is standard Nelson
_et al._ (2004), the homogeneous ${\bf q}={\bf 0}$ strain mode needs to be
handled with care in the presence of various boundary conditions. We shall
focus on isotropic loading at the boundary and neglect external shear or
compressional loading of ribbons for simplicity. This leaves us with two
possibilities, which are the

  1. (i)

Fixed stress or _isotensional_ ensemble, and the

  2. (ii)

Fixed strain or _isometric_ ensemble.

Note the latter could be realized with clamped circular boundary conditions.
Stress and strain (equivalently, force and displacement) are thermodynamically
conjugate variables and the elliptic nature of elasticity prohibits specifying
both at a boundary simultaneously. Hence, we have two natural mechanical
ensembles akin to the isobaric $(N,P,T)$ and isochoric $(N,V,T)$ ensembles of
statistical mechanics respectively, that are dual to each other. In Fig. 1, we
sketch a possible realization of the two mechanical ensembles in an atomically
or molecularly thin suspended sheet.

In the isotensional ensemble, the sheet is driven by an external isotropic
stress

$\sigma^{\rm ext}_{ij}=\sigma_{0}\delta_{ij}\;,$ (3.2)

with no further constraints on the zero modes of the displacement or strain.
As a result, the boundary can freely displace in the plane under the action of
$\sigma_{0}\neq 0$. Note that $\sigma_{0}>0$ corresponds to a tensile stress
while $\sigma_{0}<0$ is compressive stress, a situation studied in Ref.
Košmrlj and Nelson (2016). In the isometric ensemble, on the other hand, we
clamp the boundary with a fixed displacement and allow the stress to fluctuate
freely instead. If the constant displacement on the boundary is ${\bf
u}_{\mathcal{C}}=\Delta_{\mathcal{C}}\hat{\bm{\nu}}$, we have

$\oint_{\mathcal{C}}\mathrm{d}\ell\;\hat{\bm{\nu}}\cdot{\bf
u}=L_{\mathcal{C}}\Delta_{\mathcal{C}}\;,$ (3.3)

where $L_{\mathcal{C}}$ is the length of the boundary. By using Stokes’
theorem, we can rewrite this as a bulk integral,

$\dfrac{1}{A}\int\mathrm{d}^{2}r\;\bm{\nabla}\cdot{\bf u}=\epsilon\;.$ (3.4)

Here we have defined the average strain induced by the boundary as
$\epsilon=L_{\mathcal{C}}\Delta_{\mathcal{C}}/A$, where $A=\int\mathrm{d}{\bf
r}$ is the area of the sheet. To leading order this strain results in an
isotropic change in the projected area as $\delta A/A=\epsilon/2$. As before,
$\epsilon>0$ ($\Delta_{\mathcal{C}}>0$) corresponds to a uniform dilational
strain, while $\epsilon<0$ ($\Delta_{\mathcal{C}}<0$) is an isotropic
compressive strain. We can now integrate out the phonons in either ensemble to
get the free energy solely in terms of the flexural modes. In the isotensional
ensemble, we obtain

$\displaystyle\mathcal{F}_{\sigma}$
$\displaystyle=\int\mathrm{d}^{2}r\left[\vphantom{\left(\dfrac{1}{2}\right)^{2}}\dfrac{\kappa}{2}(\nabla^{2}h)^{2}+\dfrac{\sigma_{0}}{2}|\bm{\nabla}h|^{2}-\mathcal{E}h\right.$
$\displaystyle\qquad\qquad\left.+\dfrac{Y}{2}\left(\dfrac{1}{2}\mathcal{P}_{ij}^{T}\partial_{i}h\partial_{j}h\right)^{2}\right]\;.$
(3.5)

The subscript $\sigma$ denotes the fixed stress boundary condition imposed and
$\mathcal{F}_{\sigma}$ is the analogous to the Gibbs free energy. The
externally applied stress $\sigma_{0}$ enters in a quadratic term that has
been identified previously Guitter  _et al._ (1989); Košmrlj and Nelson
(2016). A tensile stress ($\sigma_{0}>0$) suppresses height fluctuations
Roldán  _et al._ (2011), while a compressive stress ($\sigma_{0}<0$) signals
the onset of the buckling instability. The Young’s modulus controls the now
well known nonlinear stretching term Nelson and Peliti (1987); Bowick and
Travesset (2001); Nelson  _et al._ (2004) via the transverse projection
operator
($\mathcal{P}_{ij}^{T}=\delta_{ij}-\partial_{i}\partial_{j}/\nabla^{2}$). In
the nonlinear stretching term (in Eq. 3.5 and below in Eq. 3.7), it is
understood that the ${\bf q}={\bf 0}$ Fourier mode has been projected out
Nelson and Peliti (1987); Bowick and Travesset (2001); Nelson  _et al._
(2004), which ensures that $\mathcal{P}^{T}_{ij}$ is well-defined. This fact
also means that the free energy is rotationally invariant when $\sigma_{0}=0$,
i.e., $h({\bf r})\to h({\bf r})+\bm{\alpha}\cdot{\bf r}$ where $\alpha_{i}$
are rotation angles, is a symmetry of the system. The average areal strain
($\epsilon=2\delta A/A$) conjugate to the imposed stress in this ensemble can
be computed from the partition function
$\mathcal{Z}_{\sigma}=\int\mathcal{D}h\;e^{-\mathcal{F}_{\sigma}/k_{B}T}$ to
give

$\langle\epsilon(\sigma_{0})\rangle=\dfrac{k_{B}T}{A}\dfrac{\partial\ln\mathcal{Z}_{\sigma}}{\partial\sigma_{0}}=\dfrac{\sigma_{0}}{B}-\dfrac{1}{2A}\int\mathrm{d}^{2}r\langle|\bm{\nabla}h|^{2}\rangle\;.$
(3.6)

The thermal average is computed using $\mathcal{F}_{\sigma}$ and
$B=\mu+\lambda$ is the bulk modulus.

In the isometric ensemble, we obtain (see Appendix A) a different result,

$\displaystyle\mathcal{F}_{\epsilon}$
$\displaystyle=\int\mathrm{d}^{2}r\left\\{\dfrac{\kappa}{2}(\nabla^{2}h)^{2}-\mathcal{E}h+\dfrac{Y}{2}\left(\dfrac{1}{2}\mathcal{P}_{ij}^{T}\partial_{i}h\partial_{j}h\right)^{2}\right.$
$\displaystyle\qquad\qquad\left.+\dfrac{B}{2}\left[\epsilon+\dfrac{1}{2A}\int\mathrm{d}^{2}r|\bm{\nabla}h|^{2}\right]^{2}\right\\}\;,$
(3.7)

where the $\epsilon$ subscript now refers to the fixed strain conditions and
now $\mathcal{F}_{\epsilon}$ is analogous to the Helmholtz free energy. As
before we include contributions from bending, the external field and the
nonlinear stretching terms. While the Young’s modulus penalizes bending
induced shear, in the presence of clamped boundaries, dilational stretching
induced by flexural deflections can no longer be accommodated by displacing
the boundary. Hence, global homogeneous dilations are a zero mode in the
isotensional ensemble, but _not_ in the isometric ensemble, and are penalized
by the bulk modulus in the latter. Upon expanding the final bracket and
neglecting an irrelevant constant, we obtain a quadratic term $\sim
B\epsilon|\bm{\nabla}h|^{2}$ that has been obtained previously Roldán  _et
al._ (2011) which mirrors the external tension term in Eq. 3.5. Importantly,
we also have an additional _nonlinear_ term of the form

$\dfrac{B}{8A}\int\mathrm{d}^{2}r\int\mathrm{d}^{2}r^{\prime}|\bm{\nabla}h|^{2}|\bm{\nabla}^{\prime}h|^{2}\;,$
(3.8)

that is independent of the strain imposed, but nonetheless arises only in
isometric ensemble. This highly nonlocal term couples the local tilts
$\sim\bm{\nabla}h$ of the membrane at arbitrarily distant points and has been
missed in previous studies Guitter  _et al._ (1988, 1989); Roldán  _et al._
(2011); Košmrlj and Nelson (2016). Anisotropic versions of this term do appear
in the description of micromechanical resonators as nonlinear beams Lifshitz
and Cross (2008) and have been included in a recent mean field analysis of a
uniaxially compressed ribbon Hanakata  _et al._ (2021). The consequences of
this term in the presence of thermal fluctuations are a major focus of this
paper.

We note some further unusual features of the new nonlocal term in Eq. 3.8.
Although it involves a double spatial integral, the whole term is extensive
(due to the factor of $1/A$), but it is importantly _not_ additive. As a
result, the membrane cannot be divided into a cumulative sum of macroscopic
parts which are roughly independent of each other in the thermodynamic limit.
Similar long-ranged interactions appear in models of compressible Sak (1974);
Bergman and Halperin (1976); de Moura  _et al._ (1976) or constrained Rudnick
_et al._ (1974) ferromagnets and can affect critical behaviour in some cases,
though without reference to ensemble inequivalence. In self-gravitating
systems Padmanabhan (1990) and mean-field models of magnets Barré  _et al._
(2001), long-ranged interactions are known to spoil the equivalence of
canonical and microcanonical ensembles, though typically in the context of
first-order phase transitions. Although the buckling instability under
compression can proceed as a continuous bifurcation, the highly nonlocal
interaction in Eq. 3.8 strongly suggests that the isotensional and isometric
ensembles may not be equivalent even in the thermodynamic limit
($A\to\infty$).

Before we proceed, we note the isometric ensemble variant of Eq. 3.6. The
average stress generated in the sheet due to the imposed strain is simply
given by

$\langle\sigma(\epsilon)\rangle=-\dfrac{k_{B}T}{A}\dfrac{\partial\ln\mathcal{Z}_{\epsilon}}{\partial\epsilon}=B\left(\epsilon+\dfrac{1}{2A}\int\mathrm{d}^{2}r\left\langle|\bm{\nabla}h|^{2}\right\rangle\right)\;,$
(3.9)

with the partition function
$\mathcal{Z}_{\epsilon}=\int\mathcal{D}h\;e^{-\mathcal{F}_{\epsilon}/k_{B}T}$
and the thermal average now performed using $\mathcal{F}_{\epsilon}$. Note
that the partition functions in the two ensembles are related,

$\mathcal{Z}_{\sigma}=\mathrm{const}.\;\int_{-\infty}^{\infty}\mathrm{d}\epsilon\;\mathcal{Z}_{\epsilon}\;e^{A\sigma_{0}\epsilon/k_{B}T}$
(3.10)

In order to incorporate both ensembles within the same calulation, we now work
with the following free energy

$\displaystyle\mathcal{F}$
$\displaystyle=\int\mathrm{d}^{2}r\left[\vphantom{\left(\dfrac{1}{2}\right)^{2}}\dfrac{\kappa}{2}(\nabla^{2}h)^{2}+\dfrac{\sigma}{2}|\bm{\nabla}h|^{2}-\mathcal{E}h\right.$
$\displaystyle\quad\left.+\dfrac{Y}{2}\left(\dfrac{1}{2}\mathcal{P}_{ij}^{T}\partial_{i}h\partial_{j}h\right)^{2}+\dfrac{v}{8A}|\bm{\nabla}h|^{2}\int\mathrm{d}^{2}r^{\prime}|\bm{\nabla}^{\prime}h|^{2}\right]\;.$
(3.11)

Upto unimportant additive constants, $\mathcal{F}=\mathcal{F}_{\epsilon}$ (Eq.
3.7) upon identifying $v=B$ and $\sigma=B\epsilon$, where $B=\mu+\lambda$ is
the bulk modulus. Alternately, if we set $v=0$ to switch off the nonlocal
nonlinear term and set $\sigma=\sigma_{0}$, then we find
$\mathcal{F}=\mathcal{F}_{\sigma}$ (Eq. 3.5). It is important to note that
setting $v=0$ is only a mathematical trick to obtain $\mathcal{F}_{\sigma}$
from Eq. 3.11. It does _not_ imply that the actual bulk modulus vanishes. As
we shall see later, the nonlocal nature of this new term guarantees that it
cannot be generated if initially absent, i.e., if we set $v=0$ in Eq. 3.11,
then it remains zero under the coarse-graining embodied in a renormalization
group transformation. Additionally we will see that the actual physical
elastic moduli ($\mu,\lambda$) renormalize identically in both ensembles, but
the buckling transition is described by two distinct fixed points (with
distinct critical exponents for some quantities!) depending on the ensemble.
So, in all that follows, we will use $v$ as a coupling constant that
distinguishes the two ensembles with $v=0$ being allowed in the isotensional
ensemble and $v>0$ only being allowed in the isometric ensemble. Only in the
latter case will $v$ be identical to the actual bulk modulus ($B$) of the
sheet.

##  4 Definition of scaling exponents near buckling

Before we analyze buckling criticality, we define our notation. Close to the
buckling transition, we expect power law scaling in a number of quantities,
the exponents for which we define below. The unprimed exponents below will
refer to the isotensional (constant stress) ensemble and the primed exponents
to the isometric (constant strain) ensemble. We will denote the buckling
threshold as $\sigma_{c}$ and $\epsilon_{c}$ in the isotensional and isometric
ensembles respectively.

###  4.1 Mechanical response

Upon approaching the buckling transition, the sheet develops a variety of
anomalous mechanical responses. In the absence of an external field
($\mathcal{E}=0$), up/down symmetry is spontaneously broken and the sheet
develops a finite $\langle h\rangle\neq 0$ when buckled. The average height
rises continuously at the transition, acting as an order parameter,

$\langle
h\rangle\propto\begin{cases}|\sigma_{0}-\sigma_{c}|^{\beta}\quad(\textrm{isotensional})\;,\\\
|\epsilon-\epsilon_{c}|^{\beta^{\prime}}\quad(\textrm{isometric})\;.\end{cases}$
(4.1)

The zero field susceptibility exhibits divergent scaling near buckling,

$\chi=\left.\dfrac{\partial\langle
h\rangle}{\partial\mathcal{E}}\right|_{\mathcal{E}=0}\propto\begin{cases}|\sigma_{0}-\sigma_{c}|^{-\gamma}\quad(\textrm{isotensional})\;,\\\
|\epsilon-\epsilon_{c}|^{-\gamma^{\prime}}\quad(\textrm{isometric})\;.\end{cases}$
(4.2)

We expect the exponents $\gamma,\gamma^{\prime}$ to be the same on either side
of the transition Goldenfeld (2018), though the amplitudes of the scaling
function can and will be different. The divergence of the susceptibility
signals the breakdown of linear response, which is also seen in the nonlinear
field dependence right at the buckling transition

$\langle
h\rangle\propto\begin{cases}\mathcal{E}^{1/\delta}\quad(\sigma_{0}=\sigma_{c}\;,\
\textrm{isotensional})\;,\\\
\mathcal{E}^{1/\delta^{\prime}}\quad(\epsilon=\epsilon_{c}\;,\,\textrm{isometric})\;.\end{cases}$
(4.3)

Finally, in conjunction with these out-of-plane responses, we also have a
concomitant nonlinear scaling in the in-plane mechanics, quantified by
anomalous stress-strain curves,

$\displaystyle\langle\epsilon\rangle$
$\displaystyle\propto\rm{const.}+(\sigma_{0}-\sigma_{c})^{1/\theta}\quad(\textrm{isotensional})\;,$
(4.4a) $\displaystyle\langle\sigma\rangle$
$\displaystyle\propto\rm{const.}+(\epsilon-\epsilon_{c})^{\theta^{\prime}}\quad(\textrm{isometric})\;,$
(4.4b)

at zero external field ($\mathcal{E}=0$). Note that the above only includes
the dominant singularity and neglects other contributions. The new exponents
$\theta,\theta^{\prime}\neq 1$ signal a violation of Hooke’s law.

###  4.2 Fluctuations and spatial scales

Apart from global quantites discussed above, local variables also develop
extended correlations when near the buckling transition. In the absence of an
external field ($\mathcal{E}=0$), the fluctuating height of the sheet has
spatial correlations with nontrival scaling properties. A nonzero external
stress or strain generically causes the height fluctuations to decay
exponentially over a finite correlation length $\xi$, albeit with different
large distance asymptotes depending on whether the sheet is buckled or flat.
This behaviour is also reflected in the normal-normal correlation function.
The unit normal to a surface specified by ${\bf X}({\bf r})=(x,y,h({\bf r}))$
in the Monge representation, is $\hat{{\bf
n}}=(-\partial_{x}h,-\partial_{y}h,1)/\sqrt{1+|\bm{\nabla}h|^{2}}$, which
allows us to simply relate the normal-normal and height-height correlation
functions as

$\langle\hat{{\bf n}}({\bf r})\cdot\hat{{\bf n}}({\bf 0})\rangle\simeq
1-\dfrac{1}{2}\left\langle\left|\bm{\nabla}h({\bf r})-\bm{\nabla}h({\bf
0})\right|^{2}\right\rangle\;,$ (4.5)

at lowest order in the height gradients. On either side of the buckling
transition, we have $\langle\hat{{\bf n}}({\bf r})\cdot\hat{{\bf n}}({\bf
0})\rangle\sim e^{-r/\xi}$, neglecting asymptotic constants and nonexponential
prefactors. The correlation length diverges at buckling as

$\xi\propto\begin{cases}|\sigma_{0}-\sigma_{c}|^{-\nu}\quad(\textrm{isotensional})\;,\\\
|\epsilon-\epsilon_{c}|^{-\nu^{\prime}}\quad(\textrm{isometric})\;.\end{cases}$
(4.6)

Right at the buckling transition, the normal correlations decay as a power law
and the sheet has critical fluctuations on all scales that cause the
correlation functions to behave anomalously. We define the translationally
invariant height and phonon correlators well away from the boundaries as

$\displaystyle G_{h}({\bf r})$ $\displaystyle=\langle h({\bf r})h({\bf
0})\rangle\;,$ (4.7) $\displaystyle\left[{\bf G}_{u}({\bf r})\right]_{ij}$
$\displaystyle=\langle u_{i}({\bf r})u_{j}({\bf 0})\rangle\;.$ (4.8)

Upon tuning to the transition, the Fourier transformed correlators [$G({\bf
q})=\int\mathrm{d}{\bf r}\;e^{-i{\bf q}\cdot{\bf r}}G({\bf r})$] exhibit a
power law scaling as $q=|{\bf q}|\to 0$. These averages define the well known
anomalous exponents $\eta$ and $\eta_{u}$ Nelson and Peliti (1987); Le Doussal
and Radzihovsky (1992); *le2018anomalous; David and Guitter (1988); Aronovitz
and Lubensky (1988); *aronovitz1989fluctuations; Bowick  _et al._ (1996)
through

$G_{h}({\bf
q})\propto\begin{cases}q^{-(4-\eta)}\quad(\sigma_{0}=\sigma_{c}\;,\,\textrm{isotensional})\;,\\\
q^{-(4-\eta^{\prime})}\quad(\epsilon=\epsilon_{c}\;,\,\textrm{isometric})\;,\end{cases}$
(4.9)

and for the in-plane phonons (irrespective of the tensor indices),

${\bf G}_{u}({\bf
q})\propto\begin{cases}q^{-(2+\eta_{u})}\quad(\sigma_{0}=\sigma_{c}\;,\,\textrm{isotensional})\;,\\\
q^{-(2+\eta^{\prime}_{u})}\quad(\epsilon=\epsilon_{c}\;,\,\textrm{isometric})\;.\end{cases}$
(4.10)

These results describe a divergent renormalization of the wave-vector
dependent bending rigidity $\kappa({\bf q})\sim q^{-\eta}$ and softening of
the elastic moduli $\mu({\bf q}),\lambda({\bf q})\sim q^{\eta_{u}}$ (with
analogous expressions with exponents $\eta^{\prime}$ and $\eta^{\prime}_{u}$
in the isometric ensemble).

##  5 Mean field theory

Here we neglect thermal fluctuations and analyze the buckling transition in
the mean-field limit, as appropriate at $T=0$. By minimizing the free energy
$\mathcal{F}$ (Eq. 3.11) over the surface profile $h({\bf r})$, we obtain the
Euler-Lagrange equation,

$\displaystyle\kappa\nabla^{4}h-\sigma\nabla^{2}h-Y\left[\mathcal{P}_{k\ell}^{T}\left(\dfrac{1}{2}\mathcal{P}_{ij}^{T}\partial_{i}h\partial_{j}h\right)\right]\partial_{k}\partial_{\ell}h$
$\displaystyle\quad-\dfrac{v}{2A}\nabla^{2}h\int\mathrm{d}^{2}r^{\prime}|\bm{\nabla}^{\prime}h|^{2}=\mathcal{E}\;,$
(5.1)

where $v=0$ and $v=B>0$ again allows us to distinguish ensembles. As both the
elastic sheet and the external load are isotropic, axisymmetry is assumed in
the following. We choose the eigenfunction of the linearized operator in Eq.
5.1 as an ansatz for the buckled height profile in a circular geometry,

$h_{0}({\bf r})=H_{0}J_{0}(q_{n}r)\;.$ (5.2)

Here, $r=|{\bf r}|$ is the distance from the center of the disc and $J_{0}(x)$
is the Bessel function of the first kind. The mode of buckling is controlled
by $q_{n}$ which is fixed by boundary conditions on the height. For
simplicity, we shall assume

$h_{0}(R)=0\implies J_{0}(q_{n}R)=0\;,$ (5.3)

where $R$ is the radius of the sheet. Other boundary conditions can also be
easily used with only minor quantitative changes in the results. A simple
Galerkin approximation Galerkin (1915) involves projecting Eq. 5.1 onto the
single mode ansatz, which then gives (the details of the calculation are
provided in Appendix B)

$q_{n}^{2}(\kappa
q_{n}^{2}+\sigma)H_{0}+q_{n}^{4}\left[\dfrac{v}{2}f(q_{n}R)+c_{0}Y\right]H_{0}^{3}=\mathcal{E}q_{n}R\;,$
(5.4)

where $c_{0}\approx 0.10567$ is a constant and $f(x)$ is a dimensionless
function given in Appendix B, with the asymptotics $f(x)\sim 2/(\pi x)$ for
$x\to\infty$. Similar to the simpler problem of the buckling of a ribbon at
$T=0$ Hanakata  _et al._ (2021), Eq. 5.4 resembles a Landau theory with
$H_{0}$ as the order parameter. It is the mechanical equivalent of a mean
field “equation of state”. Note, however, that the underlying Landau theory
has coefficients that depend on system size. In the absence of an external
field ($\mathcal{E}=0$), buckling occurs for a sufficiently negative $\sigma$
and spontaneously breaks up-down inversion symmetry. At the buckling
threshold, the lowest ($n=0$) mode (shown in Fig. 3) goes unstable first. The
buckling amplitude in either ensemble is given by

$|H_{0}|=\begin{cases}\left(\dfrac{\sigma_{c}-\sigma_{0}}{c_{0}Yq_{0}^{2}}\right)^{1/2}\quad(\sigma_{0}<\sigma_{c}\;,\,\textrm{isotensional})\;,\\\
\left(\dfrac{\epsilon_{c}-\epsilon}{c^{\prime}_{0}q_{0}^{2}}\right)^{1/2}\quad(\epsilon<\epsilon_{c}\;,\,\textrm{isometric})\;,\end{cases}$
(5.5)

Figure 3: Sketch of the first buckling mode with boundary conditions such that
$h_{0}(r=R)=0$, for a circular plate of radius $R$. The amplitude of the mode
at the center of the circular frame is $H_{0}$ and its wavevector is
$q_{0}\sim 1/R$. We assume hinged boundary conditions at $r=R$ for simplicity.
Qualitatively similar buckling modes appear for alternative boundary
conditions, such as for membranes that approach $r=R$ tangentially.

where $c^{\prime}_{0}=f(m_{0})+2c_{0}(Y/B)$ is weakly dependent on the
Poisson’s ratio through $Y/B$. As expected, we obtain a standard square root
scaling typical of pitchfork bifurcations ($\beta=\beta^{\prime}=1/2$). For
both $\sigma_{0}>\sigma_{c}$ or $\epsilon>\epsilon_{c}$, we have a stable flat
state with $H_{0}=0$. The critical stress for buckling is $\sigma_{c}=-\kappa
q_{0}^{2}$ and the critical strain is $\epsilon_{c}=-\kappa q_{0}^{2}/B$ in
the respective ensembles. The wavevector that first goes unstable with
decreasing $\sigma_{0}<0$ or $\epsilon<0$ is $q_{0}=m_{0}/R$, where
$m_{0}\approx 2.405$ is the smallest positive root of $J_{0}(m_{0})=0$ as
required by the boundary condition (Eq. 5.3). At threshold,
$f(m_{0})=J_{1}(m_{0})^{2}\approx 0.27$, is finite.

A couple of points are worth remarking on here. The buckling thresholds in
both ensembles involve compression ($\sigma_{c},\epsilon_{c}<0$) and are
$\propto 1/R^{2}$, vanishing as the area of the sheet becomes larger. Thus, in
the thermodynamic limit ($R\to\infty$), classical buckling is a thresholdless
long-wavelength ($q_{0}\sim 1/R\to 0$) instability. At the same time, the
buckling amplitude remains macroscopic ($|H_{0}|\propto R$). Note also that,
for a circular geometry, the buckled state acquires nonzero Gaussian curvature
due to the isotropic nature of the loading. This is the energetically
preferred state: a uniaxially buckled sheet that remains developable has a
higher energy for these circular loading conditions. As a result, even in the
isotensional ensemble, we pay stretching energy $\sim Y$ upon buckling (the
penalty associated with Gaussian curvature in Eq. 3.11), while in the
isometric ensemble, both the bulk ($B$) and Young’s ($Y$) moduli contribute.

Right at the transition, in the presence of an external field
($\mathcal{E}\neq 0$), we have an nonlinear response of height given by

$H_{0}=\begin{cases}\left(\dfrac{m_{0}\mathcal{E}}{c_{0}Yq_{0}^{4}}\right)^{1/3}\quad(\sigma_{0}=\sigma_{c}\;,\,\textrm{isotensional})\;,\\\
\left(\dfrac{2m_{0}\mathcal{E}}{c^{\prime}_{0}Bq_{0}^{4}}\right)^{1/3}\quad(\epsilon=\epsilon_{c}\;,\,\textrm{isometric})\;.\end{cases}$
(5.6)

As is typical of mean field models, we obtain $\delta=\delta^{\prime}=3$. Note
that, because $q_{0}\sim 1/R$, this response is strongly size dependent,
diverging as $R^{4/3}$ in a large sheet. Finally, we also have the zero field
susceptibility

$\chi=\left.\dfrac{\partial
H_{0}}{\partial\mathcal{E}}\right|_{\mathcal{E}=0}=\begin{cases}c_{\pm}\dfrac{m_{0}}{q_{0}^{2}|\sigma_{0}-\sigma_{c}|}\quad(\textrm{isotensional})\;,\\\
c_{\pm}\dfrac{m_{0}}{q_{0}^{2}B|\epsilon-\epsilon_{c}|}\quad(\textrm{isometric})\;,\end{cases}$
(5.7)

which diverges right at the transition with $\gamma=\gamma^{\prime}=1$. The
magnitude of this divergence is different on either side of the transition,
with $c_{+}=1$ above the transition and $c_{-}=1/2$ below the transition,
irrespective of the ensemble. Once again we find a strong size dependence,
with $\chi\propto R^{2}$ diverging as the area. Finally, a simple calculation
using Eqs. 3.6 and 3.9 also determines the stress-strain relation to be
$\sigma_{0}\propto\epsilon$, which sets the exponents defined by Eq. 4.4 to be
$\theta=\theta^{\prime}=1$ in both ensembles. Unlike the previous expressions,
the stress-strain relation is system size independent. Although the above
analysis was restricted to 2D membranes deforming in 3D space, all of the
mean-field results are qualitatively similar in arbitrary dimensions.

It is evident that the two ensembles have equivalent scaling behaviour in the
mean-field limit. In addition, some of the scaling functions have a nontrivial
size dependence, a feature peculiar to the buckling transition. This size
dependence is unusual from the point of view of conventional critical
phenomena Goldenfeld (2018), and as we shall see in Sec. 9, it has important
consequences for the critical exponents and their scaling relations. Below we
go beyond the mean field limit by including thermal fluctuations and show that
there are important changes.

##  6 Gaussian analysis

Here we shall primarily consider the simpler case of a flat unbuckled membrane
with $H_{0}=0$ with a vanishing symmetry-breaking external field
($\mathcal{E}=0$). We can rewrite Eq. 3.11 as
$\mathcal{F}=\mathcal{F}_{0}+\mathcal{F}_{\rm int}$ where

$\displaystyle\mathcal{F}_{0}$
$\displaystyle=\dfrac{1}{2}\int\mathrm{d}^{2}r\left[\kappa(\nabla^{2}h)^{2}+\sigma|\bm{\nabla}h|^{2}\right]\;,$
(6.1) $\displaystyle\mathcal{F}_{\rm int}$
$\displaystyle=\dfrac{Y}{2}\int\mathrm{d}^{2}r\left(\dfrac{1}{2}\mathcal{P}_{ij}^{T}\partial_{i}h\partial_{j}h\right)^{2}$
$\displaystyle\quad+\dfrac{v}{8A}\int\mathrm{d}^{2}r\int\mathrm{d}^{2}r^{\prime}|\bm{\nabla}h|^{2}|\bm{\nabla}^{\prime}h|^{2}\;.$
(6.2)

For small fluctuations, one might hope to neglect the nonlinear terms in
$\mathcal{F}_{\rm int}$. Upon Fourier transforming ($h_{{\bf
q}}\equiv\int\mathrm{d}{\bf r}\;e^{-i{\bf q}\cdot{\bf r}}h({\bf r})$), we
obtain the bare height-height correlation function

$G^{0}_{h}({\bf q})=\dfrac{1}{A}\langle|h_{{\bf
q}}|^{2}\rangle_{0}=\dfrac{k_{B}T}{\kappa q^{4}+\sigma q^{2}}\;,$ (6.3)

where $q=|{\bf q}|$ and the zero subscript denotes that the average is
performed in the noninteracting limit. We can similarly neglect all the
nonlinear interactions in $\mathcal{H}$ (Eq. 2.1) to find the bare in-plane
phonon correlation function,

$\displaystyle\left[{\bf G}^{0}_{u}({\bf q})\right]_{ij}$
$\displaystyle=\dfrac{1}{A}\langle u_{i}({\bf q})u_{j}(-{\bf q})\rangle_{0}$
$\displaystyle=\dfrac{k_{B}T}{\mu q^{2}}\mathcal{P}_{ij}^{T}({\bf
q})+\dfrac{k_{B}T}{(2\mu+\lambda)q^{2}}\mathcal{P}_{ij}^{L}({\bf q})\;,$ (6.4)

involving the longitudinal ($\mathcal{P}_{ij}^{L}({\bf q})=q_{i}q_{j}/q^{2}$)
and the transverse ($\mathcal{P}_{ij}^{T}({\bf
q})=\delta_{ij}-q_{i}q_{j}/q^{2}$) projection operators. As is evident, at the
Gaussian level, we have $\eta=\eta^{\prime}=0$ and
$\eta_{u}=\eta^{\prime}_{u}=0$ in both ensembles. From Eq. 6.3, we can easily
show that the correlation length in Gaussian limit is

$\xi=\sqrt{\dfrac{\kappa}{|\sigma|}}\propto|\sigma|^{-1/2}\;,$ (6.5)

which corresponds to $\nu=\nu^{\prime}=1/2$ in both ensembles. Note that, here
we work in the large sheet limit, which allows us to freely Fourier transform
and set $\sigma_{c},\epsilon_{c}\sim 0$.

We can now determine the importance of the nonlinear terms in Eq. 6.2 by
making a scale transformation. To do this, we rescale ${\bf r}\to b{\bf r}$,
$h\to b^{\zeta}h$ to get the following scaling dimensions for the bending
rigidity, tension and nonlinear couplings,

$y_{\kappa}=2\zeta-2\;,\quad y_{\sigma}=2\zeta\;,\quad
y_{Y}=y_{v}=4\zeta-2\;.$ (6.6)

where we have used the fact that the area $A\to b^{2}A$ under scaling. Here,
$D=2$ and $d=3$, but these scalings depend more generally on
dimensionality;their generalization for general $D$-dimensional manifolds
embedded in $d$-dimensions is given in Appendix C. In the Gaussian limit, we
have $\zeta=1$, as $h$ is simply the height with naïve dimensions of length.
This leaves the bending term scale-invariant ($y_{\kappa}=0$), but the
external load ($\sigma$), Young’s modulus ($Y$) and nonlocal coupling ($v$)
are all equally relevant perturbations for 2D membranes embedded in 3D:
$y_{\sigma}=y_{Y}=y_{v}=2>0$.

Hence, even at low temperatures when fluctuations may be small, we expect that
the nonlinear interactions eventually dominate in a large enough sheet. As
usual, a Ginzburg-like criterion determines the thermal length scale beyond
which such nonlinear fluctuations dominate Nelson and Peliti (1987); Kantor
and Nelson (1987); Košmrlj and Nelson (2016); Bowick  _et al._ (2017)

$\ell_{\rm{th}}=\sqrt{\dfrac{16\pi^{3}\kappa^{2}}{3k_{B}TY}}\;.$ (6.7)

Remarkably, at room temeprature, a monolayer of graphene or MoS2 has
$\ell_{\rm{th}}\sim 1-10~{}\rm{\AA}$, and thermal fluctuations matter already
on the atomic scale. Softer materials, such as naturally occuring 2D organic
polymers Schmidt  _et al._ (1993); Hermanson  _et al._ (2007); Klotz  _et al._
(2020) have a typical $\ell_{\rm{th}}\sim 0.1-1~{}\mu$m range which is much
larger due to their smaller Young’s moduli. As a result, the consequences of
thermal fluctuations are most dramatic in atomic crystals in contrast to the
other soft membranes. We can perturbatively account for such fluctuation
effects within a renormalization group framework that we implement below.

##  7 Perturbative Renormalization Group

We now implement a conventional Wilsonian renormalization group Goldenfeld
(2018) by iteratively integrating out a thin shell in momentum space of short
wavelength fluctuations. The cutoff in Fourier space is $\Lambda\sim 1/a$,
where $a$ is the microscopic lattice spacing. As an aside, we note that,
although the nonlocal term involving $v$ in Eq. 3.11 is quite unusual, it can
be treated straightforwardly within a standard Wilsonian treatment, as has
been done for related problems in, for example, compressible magnets Sak
(1974); Bergman and Halperin (1976); de Moura  _et al._ (1976). We perform a
systematic $\varepsilon=4-D$ expansion about the upper critical dimension
following previous works on unconstrained sheets Aronovitz and Lubensky
(1988); Guitter  _et al._ (1988). Although the full diagrammatic calculation
is presented in Appendix C, we describe the main results below. In Appendix D,
we separately provide a simple, but uncontrolled, one-loop calculation with
fixed internal ($D=2$) and external ($d=3$) dimensions that is qualitatively
correct, but numerically inaccurate.

###  7.1 Recursion relations

We carry out a perturbative low temperature expansion evaluation of thermal
fluctuations to one loop order. By integrating out fluctuations within a shell
of wavevectors $\Lambda/b\leq q\leq\Lambda$, where $b=e^{s}$ is a scale
factor, and $\Lambda^{-1}$ a short distance cut-off of order the lattice
spacing or membrane thickness, we compute corrections to the various
parameters in the model. As explained in Appendix C, even with the addition of
the new nonlinear term, the form of our elastic description in $\mathcal{F}$
(Eq. 3.11) remains unchanged at long-wavelengths under coarse-graining; only
coupling constants such as $\kappa,\sigma,Y,v$ and $\mathcal{E}$ get
renormalized. The fluctuation corrections can be cast as differential
recursion relations given below

$\displaystyle\dfrac{\mathrm{d}\kappa}{\mathrm{d}s}$
$\displaystyle=\kappa(2\zeta-\varepsilon)+\dfrac{5k_{B}T\Lambda^{2}}{192\pi^{2}(\kappa\Lambda^{2}+\sigma)}(Y+4\mu)\;,$
(7.1) $\displaystyle\dfrac{\mathrm{d}\sigma}{\mathrm{d}s}$
$\displaystyle=\sigma(2\zeta+2-\varepsilon)+\dfrac{k_{B}Tv\Lambda^{4}}{16\pi^{2}(\kappa\Lambda^{2}+\sigma)}\;,$
(7.2) $\displaystyle\dfrac{\mathrm{d}Y}{\mathrm{d}s}$
$\displaystyle=Y(4\zeta-\varepsilon)-\dfrac{5k_{B}TY^{2}\Lambda^{4}}{384\pi^{2}(\kappa\Lambda^{2}+\sigma)^{2}}\;,$
(7.3) $\displaystyle\dfrac{\mathrm{d}\mu}{\mathrm{d}s}$
$\displaystyle=\mu(4\zeta-\varepsilon)-\dfrac{k_{B}T\mu^{2}\Lambda^{4}}{96\pi^{2}(\kappa\Lambda^{2}+\sigma)^{2}}\;,$
(7.4) $\displaystyle\dfrac{\mathrm{d}v}{\mathrm{d}s}$
$\displaystyle=v(4\zeta-\varepsilon)-\dfrac{k_{B}Tv^{2}\Lambda^{4}}{16\pi^{2}(\kappa\Lambda^{2}+\sigma)^{2}}\;,$
(7.5) $\displaystyle\dfrac{\mathrm{d}\mathcal{E}}{\mathrm{d}s}$
$\displaystyle=\mathcal{E}(4-\varepsilon+\zeta)\;.$ (7.6)

The fluctuation corrections are evaluated here to leading order in
$\varepsilon=4-D$ and the codimension of the manifold is set to its physically
relevant value of $d_{c}=d-D=1$. Furthermore, as these equations are derived
in general dimension, we use the $D$-dimensional generalization of the Young’s
modulus ($Y=2\mu(2\mu+D\lambda)/(2\mu+\lambda)$) and the bulk modulus
($B=(2\mu/D)+\lambda$), which reduce to the standard 2D expressions for $D=2$,
as expected.

The renormalization equations for $\kappa$, $\mu$ and $Y$ in Eqs. 7.1, 7.4 and
7.3 are identical to those obtained previously Aronovitz and Lubensky (1988);
Guitter  _et al._ (1989), while the important coupled equations for $\sigma$
and $v$ are new results. The difference between the isometric and isotensional
ensembles captured by the presence of the nonlinear coupling $v$ is already
reflected in the modified renormalization group flows. We have also retained
the external field $\mathcal{E}$ here 111As we have set $d_{c}=1$ here, both
the height $h$ and the symmetry breaking field $\mathcal{E}$ are scalars
here.; this quantity renormalizes trivially without any graphical corrections
as it couples only to the average height ($\int\mathrm{d}{\bf r}\;h=h_{{\bf
q}={\bf 0}}$), which is left untouched by the elastic and geometric
nonlinearities.

The shear and Young’s moduli renormalize independently, as expected, but they
both contribute to the renormalization of the bending rigidity near $D=4$. One
can easily use the recursion relations for $\mu$ and $Y$ to obtain equivalent
ones for the $D$-dimensional versions of the bulk modulus
($B=(2\mu/D)+\lambda$) and the Poisson’s ratio
($\nu_{p}=\lambda/[2\mu+(D-1)\lambda]$), namely

$\displaystyle\dfrac{\mathrm{d}B}{\mathrm{d}s}$
$\displaystyle=B(4\zeta-\varepsilon)-\dfrac{k_{B}TB^{2}\Lambda^{4}}{16\pi^{2}(\kappa\Lambda^{2}+\sigma)^{2}}\;,$
(7.7) $\displaystyle\dfrac{\mathrm{d}\nu_{p}}{\mathrm{d}s}$
$\displaystyle=-\dfrac{k_{B}T\mu\Lambda^{4}}{192\pi^{2}(\kappa\Lambda^{2}+\sigma)^{2}}(1+\nu_{p})(1+3\nu_{p})\;.$
(7.8)

A couple of points are worth noting here. First, as expected, the bulk and
shear moduli also renormalize independently. Second, upon comparing Eq. 7.7
and Eq. 7.5, we immediately see that both $v$ and $B$ renormalize in identical
ways, guaranteeing that in the isometric ensemble, since $v=B$ at the
microscopic scale, they remain equal on larger scales as well. In contrast,
the isotensional ensemble is characterized by $v=0$ (which remains invariant
under renormalization), even though $B>0$.

The third important point concerns the Poisson’s ratio $\nu_{p}$ 222Note that,
we define the Poisson’s ratio simply using the elastic moduli here, as is
usually done Le Doussal and Radzihovsky (2018); Nelson  _et al._ (2004).
Alternate, yet related, definitions are possible that can yield different
results, see Ref. Burmistrov  _et al._ (2018b) for instance.. As is easily
seen from Eq. 7.8, we have a stable fixed point where
$\mathrm{d}\nu_{p}/\mathrm{d}s=0$ at $\nu_{p}=-1/3$ ($\nu_{p}=-1$ is
unphysical as it corresponds to $\lambda=-2\mu/D$ leading to a marginally
stable solid with vanishing bulk and Young’s moduli), which is exactly the
one-loop estimate for the universal Poisson’s ratio of a free standing elastic
membrane, in accord with previous self consistent calculations Le Doussal and
Radzihovsky (1992); *le2018anomalous and Monte-Carlo simulations Bowick  _et
al._ (1996); Falcioni  _et al._ (1997); *bowick2001universal; Cuerno  _et al._
(2016). This universal auxetic response is a characteristic property of the
flat phase of unconstrained thermalized membranes Bowick  _et al._ (2001);
Bowick and Travesset (2001). The reason we obtain this result from a simple
one-loop calculation near $D=4$ is that the structure of the one-loop
calculation is the same as that of the self-consistent calculation done by Le
Doussal and Radzihovsky (1992), and the higher order box diagrams are
convergent in that case. Furthermore, the universal Poisson’s ratio obtained
is independent of both internal and embedding dimensions of the membrane Le
Doussal and Radzihovsky (1992), hence we recover the same value even in an
$\varepsilon=4-D$ expansion.

To analyze these recursion relations, we introduce the following dimensionless
variables (the Poisson’s ratio is of course already dimensionless),

$\displaystyle
K=\dfrac{\kappa\Lambda^{2}}{\kappa\Lambda^{2}+\sigma}\;,\quad\bar{Y}=\dfrac{k_{B}T\Lambda^{4}}{(\kappa\Lambda^{2}+\sigma)^{2}}Y\;,$
$\displaystyle\bar{\mu}=\dfrac{k_{B}T\Lambda^{4}}{(\kappa\Lambda^{2}+\sigma)^{2}}\mu\;,\quad\bar{v}=\dfrac{k_{B}T\Lambda^{4}}{(\kappa\Lambda^{2}+\sigma)^{2}}v\;,$
(7.9)

which are appropriate near $D=4$ dimensions. For general $D$, we must replace
the factor of $\Lambda^{4}$ by $\Lambda^{D}$ to keep $\bar{Y}$, $\bar{\mu}$
and $\bar{v}$ dimensionless. As the external field $\mathcal{E}$ does not
influence the fixed points, we will set $\mathcal{E}=0$ for now. The effects
of $\mathcal{E}\neq 0$ will be addressed within a general scaling theory we
develop in Sec. 9.

The physical interpretation of $K$ is that it measures the relative importance
of bending to the external load. We note that $0\leq K<\infty$ and then
demarcate three distinct regimes based on the value for $K$ as follows:

  * •

$0<K<1$ : Tension dominated ($\sigma>0$),

  * •

$K\sim 1$ : Bending dominated ($\sigma\sim 0$),

  * •

$K>1$ : Compression dominated ($\sigma<0$).

The buckled state we found using mean-field theory thus occurs for $K>1$ in
the presence of compression. The recursion relations for these dimensionless
coupling constants then read

$\displaystyle\dfrac{\mathrm{d}K}{\mathrm{d}s}$
$\displaystyle=2(K-1)\left[K-\dfrac{5}{384\pi^{2}}(\bar{Y}+4\bar{\mu})\right]-\dfrac{\bar{v}K}{16\pi^{2}}\;,$
(7.10) $\displaystyle\dfrac{\mathrm{d}\bar{Y}}{\mathrm{d}s}$
$\displaystyle=\left[\varepsilon+4(K-1)-\dfrac{\bar{v}}{8\pi^{2}}-\dfrac{25\bar{Y}}{384\pi^{2}}-\dfrac{5\bar{\mu}}{24\pi^{2}}\right]\bar{Y}\;,$
(7.11) $\displaystyle\dfrac{\mathrm{d}\bar{\mu}}{\mathrm{d}s}$
$\displaystyle=\left[\varepsilon+4(K-1)-\dfrac{\bar{v}}{8\pi^{2}}-\dfrac{5\bar{Y}}{96\pi^{2}}-\dfrac{7\bar{\mu}}{32\pi^{2}}\right]\bar{\mu}\;.$
(7.12) $\displaystyle\dfrac{\mathrm{d}\bar{v}}{\mathrm{d}s}$
$\displaystyle=\left[\varepsilon+4(K-1)-\dfrac{3\bar{v}}{16\pi^{2}}-\dfrac{5\bar{Y}}{96\pi^{2}}-\dfrac{5\bar{\mu}}{24\pi^{2}}\right]\bar{v}\;.$
(7.13) $\displaystyle\dfrac{\mathrm{d}\nu_{p}}{\mathrm{d}s}$
$\displaystyle=-\dfrac{k_{B}T\mu\Lambda^{4}}{192\pi^{2}(\kappa\Lambda^{2}+\sigma)^{2}}(1+\nu_{p})(1+3\nu_{p})\;.$
(7.14)

As expected, the scale factor $\zeta$ for the flexural phonon field drops out
of the equations when cast in terms of dimensionless variables. Note that,
while we quote the renormalization group flow equations for
$\bar{Y},\bar{\mu}$ and $\nu_{p}$ (Eqs. 7.11, 7.12 and 7.14), only two of the
three are independent, as $\nu_{p}=(Y-2\mu)/[2\mu+(D-2)Y]$. For an
unconstrained membrane, $K$ remains fixed at unity, as $\sigma=\bar{v}=0$.
While in the isotensional ensemble, we can have any $\sigma\neq 0$ ($K\neq 1$)
with $\bar{v}=0$, in the isometric ensemble, we generally have _both_
$\sigma\neq 0$ and $\bar{v}\neq 0$. In the latter case, the system
spontaneously develops a thermally generated tension due to the geometric
confinement enforced by the clamped boundary conditions, as discussed below.
But first, we analyze the fixed points of the recursion relations.

###  7.2 Fixed points

Figure 4: A schematic of the full renormalization group flow diagram in the
three-dimensional parameter space of $\\{K,\bar{Y},\bar{v}\\}$. We fix the
Poisson’s ratio to its universal value $\nu_{p}=-1/3$ here, so that both the
shear and bulk moduli are determined by the Young’s modulus through
$\bar{\mu}=(D+1)\bar{Y}/4$ and $\bar{B}=(D+1)\bar{Y}/D(D+2)$. The three
isotensional fixed points (Gσ, Gκ and vKth) are shown as red points with
$0\leq K\leq 1$ and $\bar{v}=0$, while the new constrained fixed point CvKth
with $K>1$ and $\bar{v}\neq 0$ is shown in blue. An unphysical fixed point
with $\bar{v}\neq 0$ and $\bar{Y}=\bar{B}=\bar{\mu}=0$ is also present as a
green dot at the bottom, but this is irrelevant for our purposes. The red
plane at $\bar{v}=0$ on the left shows the accessible space of coupling
parameters within the often used isotensional ensemble. Under fixed stress
(isotensional) conditions, the thermal buckling transition occurs at
$\sigma=0$ (i.e., $K=1$) and is controlled by the conventional vKth fixed
point. However, for fixed strain (isometric) conditions, when $\bar{v}\neq 0$,
we flow instead to a new codimension-1 fixed point (CvKth) that now controls
the thermal buckling transition. The unstable renormalization group flow going
towards large $K>1$ corresponds to strong compression and postbuckling
behaviour. At low temperature, Gσ is a globally attracting and stable fixed
point which controls the properties of a tense flat membrane for both
ensembles.

We now enumerate the four physically relevant fixed points 333There are other
interacting fixed points present, including two that correspond to a fluid
membrane (with a vanishing shear modulus $\mu=0$) and a physically
unrealizable one ($B=0$ but $v\neq 0$), that are not relevant for us.
permitted in both ensembles, to $\mathcal{\mathcal{O}}(\varepsilon)$:

  1. (i)

Gσ: $K_{*}=0$, $\bar{Y}_{*}=0$, $\bar{\mu}_{*}=0$, $\bar{v}_{*}=0$.

  2. (ii)

Gκ: $K_{*}=1$, $\bar{Y}_{*}=0$, $\bar{\mu}_{*}=0$, $\bar{v}_{*}=0$.

  3. (iii)

vKth: $K_{*}=1$, $\bar{Y}_{*}=77\pi^{2}\varepsilon/25$,
$\bar{\mu}_{*}=96\pi^{2}\varepsilon/25$, $\bar{v}=0$ ($\nu_{p}=-1/3$).

  4. (iv)

CvKth: $K_{*}=1+\varepsilon/50$, $\bar{Y}_{*}=77\pi^{2}\varepsilon/25$,
$\bar{\mu}_{*}=96\pi^{2}\varepsilon/25$, $\bar{v}_{*}=16\pi^{2}\varepsilon/25$
($\nu_{p}=-1/3$).

We have set $d_{c}=1$ here as is physically relevant; the expressions for the
fixed points with arbitrary $d_{c}$ are presented in Appendix C. Of these
fixed points, only Gκ, Gσ and vKth are admissable in the isotensional
ensemble. The thermal Föppl-von Kármán fixed point vKth has been the focus of
virtually all studies to date. The _constrained_ thermal fixed point CvKth is
new and unique to the isometric ensemble. Both Gκ and Gσ are Gaussian
(noninteracting) fixed points that are bending and tension dominated
respectively. The conventional flat phase is described by vKth, and occurs for
a vanishing renormalized tension (hence $K_{*}=1$) that is appropriate for an
unconstrained fluctuating membrane. This fixed point has been extensively
studied previously Nelson and Peliti (1987); Le Doussal and Radzihovsky
(1992); Aronovitz and Lubensky (1988); Guitter  _et al._ (1988, 1989); Roldán
_et al._ (2011); Košmrlj and Nelson (2016) and it controls the buckling
transition in the absence of boundary constraints, i.e., the _isotensional_
ensemble as evidenced by $\bar{v}=0$ at the fixed point.

In contrast, a new constrained fixed point CvKth emerges in the isometric
ensemble with $\bar{v}\neq 0$, reflecting the geometric constraint imposed by
the clamped boundaries in the isometric ensemble. The new interacting fixed
point involves bare compression (as $K_{*}>1$), unlike the others, reflecting
the presence of a fluctuation induced spontaneous tension that appears only
when the boundary is constrained. We will discuss this feature in more detail
in Sec. 8. As we will show below, CvKth controls the buckling transition in
the _isometric_ ensemble.

We note that both vKth and CvKth are characterized by the universal Poisson’s
ratio $\nu_{p}=-1/3$. A schematic of the full renormalization group flow
diagram with the above fixed points is sketched in Fig. 4. The stability of
each fixed point can be analyzed by linearizing about it. The two Gaussian
fixed points, Gσ and Gκ, differ simply by the presence or absence of $\sigma$
and play a role in both ensembles. The tension-dominated fixed point Gσ is a
globally stable attractor that controls the low temperature phase of a tense
flat membrane. On the other hand, the bending dominated Gaussian fixed point
Gκ is unstable in all directions, as found previously for tense membranes
Guitter  _et al._ (1989); Roldán  _et al._ (2011); Košmrlj and Nelson (2016).
Both these fixed points along with their flow directions in different
invariant planar sections of the full parameter space, are shown in Fig. 5.

Within the invariant subspace of $\bar{v}=0$, associated with the isotensional
ensemble (see Fig. 4), the conventional flat phase fixed point vKth has only
one relevant direction corresponding to the external stress. As a result, it
controls the finite temperature buckling transition in the isotensional
ensemble, where the stress is tuned to zero at threshold. But once we allow
for a $\bar{v}\neq 0$, i.e., work in the isometric ensemble instead, we find
that vKth is now a codimension two fixed point, being _unstable_ to this new
nonlinear coupling. As sketched in Fig. 4, the system can now flow to the new
constrained fixed point CvKth instead, which has codimension one.
Perturbations in $\bar{Y},\bar{B}$ and $\bar{v}$ are all irrelevant at CvKth
and the only unstable or relevant direction is primarily along $K$. In other
words, within the isometric ensemble, the strain tuned buckling transition is
controlled by CvKth and _not_ vKth.

The identification of two distinct ensemble-dependent fixed points controlling
buckling is a significant achievement. The fact that the choice of fixed point
is picked by the mechanical ensemble, here decided by fixed strain or stress
boundary conditions, is quite intriguing. Although the isotensional and
isometric ensembles are dual to each other, they remain inequivalent even in
the thermodynamic limit, due to flexural phonons on all length scales at the
critical point. As mentioned in the introduction, this remarkable feature is
akin to Fisher renormalization of conventional critical exponents, which we
demonstrate explicitly in Sec. 9. Below we compute the flat to buckled phase
boundary and analyze the linearized flow in the vicinity of the two
interacting fixed points to extract critical exponents for the buckling
transition.

##  8 Buckling transition in $4-\varepsilon$ dimensions

Figure 5: A schematic of the renormalization group flows projected onto the
invariant attracting planes appropriate to the two ensembles. The Poisson’s
ratio is fixed to its universal value $\nu_{p}=-1/3$ in both cases. (a) In the
isotensional ensemble, $\bar{v}=0$ identically, and we have three fixed
points, Gσ, Gκ and vKth. The red lines are separatrices that connect the
various fixed points and demarcate their basins of attraction. The vertical
separatrix flowing into vKth at $K=1$ decribes the buckling transition in this
ensemble. The black streamlines are integral curves of the flow and the blue
curves highlight two representative trajectories bracketing the buckling
transition. (b) In the isometric ensemble, the relevant attractor is now a
plane characterized by $\bar{Y}=z_{*}\bar{v}$ with
$z_{*}=24/5+\mathcal{O}(\varepsilon)$. Once again, we have three important
fixed points, Gσ, Gκ and CvKth. The red curves are separatrices delimiting
stability basins for each fixed point and the separatrix flowing into CvKth
controls the buckling threshold in this ensemble. Unlike in the isotensional
case, this line is curved and bends towards $K>1$, signalling the generation
of spontaneous tension (Fig. 6). The black curves are streamlines of the local
renormalization group flow and the blue curves are representative flow
trajectories on either side of the buckling transition. In both ensembles, Gσ
is globally stable for flat and tense mebranes ($K<1$), while flows towards
larger values of $K(>1)$ lead to strong compression and buckling.

Let us now analyze the recursion relations given in Eqs. 7.10-7.14 in more
detail. For a given $\bar{Y}$, as $\bar{\mu}$ and $\nu_{p}$ are related, we
only have to consider one of them. From Eq. 7.14, we easily see that the fixed
point at $\nu_{p}=-1/3$ is stable and exponentially attracting for any finite
$\bar{Y}$. So we shall neglect perturbations in the Poisson’s ratio and fix
$\nu_{p}=-1/3$. This condition in turn fixes the shear modulus to be
$\bar{\mu}=(D+1)\bar{Y}/4$ and the bulk modulus to be
$\bar{B}=(D+1)\bar{Y}/D(D+2)$, allowing us to then concentrate on the flow in
the three-dimensional subspace of just $\\{K,\bar{Y},\bar{v}\\}$ parametrizing
the stable attractor.

As an aside, note that $\nu_{p}$ rapidly approaches its fixed point value of
$-1/3$ only when $\bar{Y}>0$, which is true in the vicinity of both the vKth
and CvKth fixed points. In contrast, for a tense membrane governed by Gσ,
$\bar{Y}(s)\propto e^{-s(4-\varepsilon)}\to 0$ approaches zero exponentially
fast on large scales. In this case, $\nu_{p}$ does not reach its universal
fixed point value and instead, the rapidly vanishing $\bar{Y}$ essentially
freezes $\nu_{p}$ at a value that depends on microscopic properties of the
material. Thus, while the large-distance Poisson’s ratio is universal at the
buckling transition, away from it, in a tense mebrane, it becomes nonuniversal
and depends on microscopic details, consistent with results for fluctuating
membranes under strong tension Roldán  _et al._ (2011); Košmrlj and Nelson
(2016); Los  _et al._ (2009); Burmistrov  _et al._ (2018a).

We shall now address buckling criticality in the two ensembles separately
below.

###  8.1 Isotensional ensemble

In the isotensional ensemble, we set both $\sigma=\sigma_{0}$ and $\bar{v}=0$.
The latter picks out a renormalization group invariant plane (see Figs. 4 and
5a) specific to this ensemble. The buckling transition then occurs at
$\sigma_{0}=\sigma_{c}=0$ (in the thermodynamic limit of infinite system
size), though $\sigma_{c}$ is nonzero for a finite sheet (see Sec. 5). Note
that, right at the unconstrained fixed point vKth, we do have $\sigma_{0}=0$
(i.e., $K=1$) even at finite temperature, a result that holds to all orders in
perturbation theory. Hence the critical stress at the buckling transition is
still given by its $T=0$ value,

$\sigma_{c}(T)=-\kappa q_{0}^{2}\;,$ (8.1)

and it does _not_ receive corrections from thermal fluctuations in the
isotensional ensemble. As before $q_{0}\sim 1/R$ is the smallest wavevector
determined by the boundary conditions and the size of the sheet (Sec. 5).

To compute anomalous scaling exponents at the transition, we use a standard
renormalization group matching procedure Rudnick and Nelson (1976) to relate
correlation functions evaluated near the transition to those further away from
the critical point. Under scaling, ${\bf r}=b{\bf r}^{\prime}$ and $h({\bf
r})=b^{\zeta}h({\bf r}^{\prime})$, conversely, ${\bf q}={\bf q}^{\prime}/b$
and $h_{{\bf q}}=h^{\prime}_{{\bf q}^{\prime}}b^{D+\zeta}$ (in
$D$-dimensions). Upon setting $b=e^{s}$, we have the following scaling
relation for the height-height correlation function ($G_{h}({\bf
q})\equiv\langle|h_{{\bf q}}|^{2}\rangle/V_{D}$, where $V_{D}$ is the
$D$-dimensional volume of the manifold)

$G_{h}({\bf
q})=\exp\left\\{\int_{0}^{s}\mathrm{d}s^{\prime}[D+2\zeta(s^{\prime})]\right\\}G_{h}({\bf
q}e^{s};s)\;,$ (8.2)

where $G_{h}({\bf k};s)$ is computed using all the parameters evaluated at
scale $s$. We now choose $s=s^{*}$ such that $|{\bf
q}|e^{s^{*}}=\ell_{\rm{th}}^{-1}$, set by the thermal length (Eq. 6.7). This
condition allows us to write,

$\displaystyle G_{h}({\bf q})$
$\displaystyle=\exp\left[Ds^{*}+2\int_{0}^{s^{*}}\mathrm{d}s\;\zeta(s)\right]G_{h}(\ell_{\rm{th}}^{-1};s^{*})$
$\displaystyle=k_{B}T\ell_{\rm{th}}^{4}\dfrac{K(s^{*})}{\kappa(s^{*})}\exp\left[Ds^{*}+2\int_{0}^{s^{*}}\mathrm{d}s\;\zeta(s)\right]\;.$
(8.3)

Here we have used the fact that on small scales
($\ell\lesssim\ell_{\rm{th}}$), fluctuation corrections are negligible and a
Gaussian or mean field treatment is valid. To evaluate the renormalized
bending rigidity at scale $s^{*}$, we need the following flow equation as well

$\dfrac{\mathrm{d}\kappa}{\mathrm{d}s}=\kappa\left[2\zeta-\varepsilon+\dfrac{5(\bar{Y}+4\bar{\mu})}{192\pi^{2}K}\right]\;.$
(8.4)

It is convenient to chose the height rescaling factor $\zeta(s)$ to keep
$\kappa(s)$ fixed, which leads to

$\zeta(s)=\dfrac{\varepsilon}{2}-\dfrac{5}{384\pi^{2}K(s)}\left[\bar{Y}(s)+4\bar{\mu}(s)\right]\;.$
(8.5)

Right at the buckling transition, the coupling constants flow to the fixed
point vKth. The height correlator defines a renormalized bending rigidity via

$\kappa_{R}({\bf q})^{-1}=q^{4}\dfrac{G_{h}({\bf q})}{k_{B}T}\;.$ (8.6)

Upon using Eq. 8.4 and Eq. 8.3 right at buckling, this then gives the well-
known diverging bending rigidity,

$\kappa_{R}({\bf
q})=\kappa(q\ell_{\rm{th}})^{-\eta}\;,\quad\eta=\dfrac{12}{25}\varepsilon+\mathcal{O}(\varepsilon^{2})\;.$
(8.7)

The anomalous exponent $\eta$ that we obtain matches earlier calculations
performed for unconstrained membranes Aronovitz and Lubensky (1988); Guitter
_et al._ (1988). While we don’t expect the one-loop approximation of $\eta$ to
be numerically accurate in the physically relevant case of $D=2$ dimensions,
we nonetheless obtain a reasonable value of $\eta=24/25\approx 0.96$ upon
setting $\varepsilon=2$ in Eq. 8.7. More sophisticated calculations involving
self consistent Le Doussal and Radzihovsky (1992); *le2018anomalous and
nonperturbative techniques Kownacki and Mouhanna (2009) give $\eta\simeq
0.82-0.85$ which compares well with the exponent value measured in numerical
simulations Bowick  _et al._ (1996); Los  _et al._ (2009); Roldán  _et al._
(2011); Bowick  _et al._ (2017) and recent experiments Blees  _et al._ (2015).

The elastic moduli also experience a scale dependent renormalization, though
they now get _softer_ on larger scales. The renormalized Young’s modulus
scales as

$Y_{R}({\bf
q})=Y(q\ell_{\rm{th}})^{\eta_{u}}\;,\quad\eta_{u}=\dfrac{\varepsilon}{25}+\mathcal{O}(\varepsilon^{2})\;,$
(8.8)

The other elastic moduli ($\mu,\lambda$ and $B$) all scale in the same way,
with the same $\eta_{u}$ exponent. Both $\eta$ and $\eta_{u}$ are related via
a Ward identity $\eta_{u}+2\eta=\varepsilon$, which is a consequence of
rotational invariance Guitter  _et al._ (1988, 1989); Aronovitz and Lubensky
(1988); Aronovitz  _et al._ (1989); Le Doussal and Radzihovsky (1992).

If the external tension is small but nonzero, then we perturb slightly away
from vKth. We write $K=1+\delta K$ and linearize in $\delta K$ to get

$\dfrac{\mathrm{d}\delta
K}{\mathrm{d}s}=\left(2-\dfrac{12}{25}\varepsilon\right)\delta
K\equiv(2-\eta)\delta K\;,$ (8.9)

which is exactly true by virtue of the definition of $\eta$. For a small
external stress ($|\sigma_{0}|\ll\kappa\Lambda^{2}$), $K\approx 1$, with
$\delta K(s)=\delta K(0)e^{(2-\eta)s}$ growing with $s$, as expected of a
relevant perturbation. Eventually, we reach a scale $s^{*}$ at which the
external stress has grown large and is now comparable to the bending rigidity,
i.e., $|\sigma(s^{*})|\approx\kappa(s^{*})\Lambda^{2}$), which defines the
correlation length $\xi$ via $s^{*}=\ln(\xi/a)$ ($a\sim 1/\Lambda$ is a
lattice cutoff) beyond which the stress dominates bending. After incorporating
a nonzero $\sigma_{c}$ appropriate to a finite system, we have $\delta
K(0)\propto|\sigma_{0}-\sigma_{c}|$ which gives

$\xi\propto|\sigma_{0}-\sigma_{c}|^{-\nu}\;,\quad\nu=\dfrac{1}{2-\eta}=\dfrac{1}{2}+\dfrac{3\varepsilon}{25}+\mathcal{O}(\varepsilon^{2})\;,$
(8.10)

for the isotensional ensemble. On short length scales ($\ell\ll\xi$), the
system is controlled by the bending-dominated vicinty of the vKth fixed point,
while on larger scales ($\ell\gg\xi$) when $\sigma_{0}>0$, the system is
dominated by external tension and the Gσ fixed point. Hence, for
$\sigma_{0}>0$, we have the following length scale dependence (for $D=2$)

$\dfrac{\kappa_{R}(\ell)}{\kappa}\propto\begin{cases}1\;,\quad(\ell<\ell_{\rm{th}})\\\
(\ell/\ell_{\rm{th}})^{\eta}\;,\quad(\ell_{\rm{th}}<\ell<\xi)\\\
(\xi/\ell_{\rm{th}})^{\eta}\ln(\ell/\xi)\;,\quad(\ell>\xi)\end{cases}\;,$
(8.11)

with similar expressions for the elastic moduli Košmrlj and Nelson (2016). A
similar dependence without the $\ln(\ell/\xi)$ term also exists for $2<D<4$
dimensions.

###  8.2 Isometric ensemble

Now that we have recapitulated the isotensional results, let us move onto the
more interesting isometric ensemble. We now set $\sigma=B\epsilon$ and
identify $\bar{v}=\bar{B}$. As mentioned earlier, the $\bar{v}=0$ plane
defining the isotensional ensemble is unstable to finite $\bar{v}$
perturbations, leading us to consider the full 3D space of parameters
$\\{K,\bar{Y},\bar{v}\\}$. For $\bar{v}>0$, i.e., in the isometric ensemble
now, we can further reduce dimensionality by writing $\bar{Y}=z\bar{v}$, which
gives

$\dfrac{\mathrm{d}z}{\mathrm{d}s}=\dfrac{\bar{Y}}{16\pi^{2}}\left(1-\dfrac{5z}{24}\right)\;,$
(8.12)

to leading order in $\mathcal{O}(\varepsilon)$. With $\bar{Y}>0$ as before,
this equation has an exponentially stable fixed point given by
$z=24/5+\mathcal{O}(\varepsilon)$ that attracts all the renormalization group
flows for $\bar{v}>0$. This stable fixed point in $z$ hence defines an
attracting invariant plane in the $\\{K,\bar{Y},\bar{v}\\}$ space (see Fig.
5b), only accessible within the isometric ensemble. Note that the new
constrained fixed point CvKth lies on this plane as well, a welcome feature
that guarantees that when $\bar{v}\neq 0$ microscopically and equal to
$\bar{B}$ on short scales (as it should be in the isometric ensemble), their
equality is retained on larger scales. More importantly, the relation
$\bar{v}=\bar{B}>0$ serves as an invariant attractor under coarse-graining in
the isometric ensemble. Within this plane, we have two coupled flow equations,
namely

$\displaystyle\dfrac{\mathrm{d}K}{\mathrm{d}s}$
$\displaystyle=2K(K-1)+\dfrac{\bar{v}}{16\pi^{2}}(12-13K)\;,$ (8.13)
$\displaystyle\dfrac{\mathrm{d}\bar{v}}{\mathrm{d}s}$
$\displaystyle=\bar{v}\left[\varepsilon+4(K-1)-\dfrac{27\bar{v}}{16\pi^{2}}\right]\;,$
(8.14)

as shown in Fig. 5b.

Upon dividing the above two equations, we obtain
$\mathrm{d}K/\mathrm{d}\bar{v}$ which we numerically integrate to obtain the
basin of attraction of Gσ and CvKth. The stable and unstable manifolds are
obtained as integral curves of the stable and unstable eigendirections at
CvKth, which correspond to separatrices shown in red in Fig. 5b. The
separatrix connects the unstable fixed point Gκ at $K=1$ to the constrained
thermal buckling transition fixed point at CvKth, and it delineates the
stability region for a clamped membrane. All parameter values that fall to the
left of this separatrix flow into the stable Gσ fixed point, leading to a
sheet that is flat and tense on large scales. In the opposite case, parameter
values starting to the right of the red separatrix flow away to larger values
of $K$, signalling strong compression ($\sigma=B\epsilon<0$) and buckling of
the membrane. Representative flow trajectories illustrating this are shown in
blue in Fig. 5b. The separatrices are computed as the solution of a boundary
value problem and hence don’t admit a simple analytical solution. Nonetheless,
we can obtain some asymptotic results informed by the algebraic structure of
the recursion relations, in conjunction with the numerical solution.

Figure 6: The numerically computed buckling threshold $|\epsilon_{c}(T)|$
extrapolated to $\varepsilon=2$ ($D=2$). The critical buckling strain gets
shifted to more negative (compressive) values at finite temperature due to the
generation of a thermally induced spontaneous tension. We have subtracted out
the zero temperature buckling threshold $\epsilon_{c}(0)=-\kappa q_{0}^{2}/B$.
Near the buckling transition, we have used the universal Poisson’s ratio
($\nu_{p}=-1/3$) to relate the bulk and Young’s modulus in $D=2$. At low
temperature (large $\ell_{\rm th}$), $\epsilon_{c}(T)\sim T\ln T$, while at
higher temperature (smaller $\ell_{\rm th}$), we have
$\epsilon_{c}(T)\sim-\sqrt{T}$ (Eq. 8.15). While this plot is obtained by
setting $\varepsilon=2$ in the resursion relations obtained by expanding
around $D=4$ dimensions, a qualitatively similar curve is obtained from a
fixed dimension calculation with $D=2$ and $d=3$ (Appendix D).

Upon using the definitions of $K$ and $\bar{B}$ and the relation
$\sigma=B\epsilon$, we obtain the critical curve for the buckling transition
strain ($\epsilon_{c}(T)$) in terms of the elastic constants and temperature,
as plotted in Fig. 6 for $D=2$ ($\varepsilon=2$). For low and high
temperatures, we find simple asymptotic expansions for the buckling threshold
$\epsilon_{c}(T)<0$ (in $D=2$),

$\displaystyle|\epsilon_{c}(T)|$ $\displaystyle\simeq\dfrac{\kappa
q_{0}^{2}}{B}+\dfrac{k_{B}T}{8\pi\kappa}\left[2\ln\left(\dfrac{a}{\ell_{\rm{th}}}\right)+c_{1}\right]\;,\quad(\ell_{\rm{th}}\gg
a)\;,$ (8.15a) $\displaystyle|\epsilon_{c}(T)|$
$\displaystyle\simeq\dfrac{\kappa
q_{0}^{2}}{B}+\dfrac{k_{B}T}{8\pi\kappa}\left(\dfrac{\ell_{\rm{th}}}{a}\right)\left[c_{2}-c_{3}\dfrac{\ell_{\rm{th}}}{a}\right]\;,\quad(\ell_{\rm{th}}\ll
a)\;,$ (8.15b)

where $a\sim\Lambda^{-1}$ is the lattice cutoff and $c_{1,2,3}$ are numerical
constants that must be computed by numerical integration of the recursion
relations. While here, we extrapolated our perturbative solution to
$\varepsilon=2$, we have confirmed that the same asymptotic expressions for
the buckling strain, with only $c_{1,2,3}$ modified, are also obtained within
a fixed dimension calculation with $D=2$ and $d=3$ (Appendix D). As $T\to 0$,
$\epsilon_{c}(T)\to-\kappa q_{0}^{2}/B$, which is the zero temperature
buckling instability threshold with $q_{0}\sim 1/R$ being the smallest
available mode in a system of size $R$ (Sec. 5). We have utilized the fact
that, near the transition $\nu_{p}=-1/3$, which relates $Y$ and $B$ via
$B=3Y/8$ (in $D=2$), allowing us to write the above in terms of the thermal
length $\ell_{\rm{th}}$. As $\ell_{\rm{th}}\sim T^{-1/2}$ (Eq. 6.7), for $D=2$
we have $|\epsilon_{c}(T)|\sim T\ln(1/T)$ at low temperature and
$|\epsilon_{c}(T)|\sim\sqrt{T}$ at high temperature, as shown in Fig. 6. For
general $D$, as $T\to 0$, we find that $|\epsilon_{c}(T)|\sim T[1+{\rm
const.}\;(T^{2/\varepsilon-1}-1)/(2-\varepsilon)]$ and the linear $T$
dependence dominates at small $T$ for $0<\varepsilon<2$ ($2<D<4$), while an
additional logarithmic term $\sim T\ln(T)$ appears when $\varepsilon=2$
($D=2$). Note that the high temperature asymptotics depends only weakly on
dimension.

We pause here to comment on the above results. Unlike in the isotensional
ensemble, where the buckling threshold did _not_ receive any correction from
thermal fluctuations (Eq. 8.1), in the isometric ensemble, the buckling
threshold gets pushed to higher values of compression (as $\epsilon_{c}(T)<0$
and $|\epsilon_{c}(T)|$ increases with $T$) at higher temperature. In other
words, the sheet spontaneously develops a tension due to thermal fluctuations
in the isometric ensemble. A freely fluctuating sheet wants to shrink in-plane
for entropic reasons, but the clamped boundaries resist this shrinkage,
thereby putting the sheet under tension. As a result, the externally imposed
strain now has to compensate and overcome this thermally induced tension in
order to cause buckling. This effect is absent in the isotensional ensemble
because the boundaries are free to fluctuate, allowing the sheet to freely
shrink with increasing temperature, albeit against a constant external stress.

We now compute critical scaling exponents at the buckling transition. Here, by
tuning right to the buckling threshold, we approach the constrained fixed
point CvKth. Upon evaluating the height-height correlator, we obtain the
renormalized bending rigidity to be

$\kappa_{R}({\bf
q})=\kappa(q\ell_{\rm{th}})^{-\eta^{\prime}}\;,\quad\eta^{\prime}=\dfrac{12}{25}\varepsilon+\mathcal{O}(\varepsilon^{2})\;.$
(8.16)

Remarkably, we obtain the same anomalous scaling exponent here as in the
isotensional ensemble (Eq. 8.7). As we will show later in Sec. 4, we expect
$\eta=\eta^{\prime}$ from general scaling arguments, which we also confirm for
arbitrary $d_{c}$ within a lowest order systematic $\varepsilon=4-D$ expansion
in Appendix C and Table 1.

A similar analysis of the phonon correlator or equivalently the nonlinear
stretching term also provides the renormalized Young’s modulus,

$Y_{R}({\bf
q})=Y(q\ell_{\rm{th}})^{\eta^{\prime}_{u}}\;,\quad\eta^{\prime}_{u}=\dfrac{\varepsilon}{25}+\mathcal{O}(\varepsilon^{2})\;.$
(8.17)

As before, $\eta^{\prime}_{u}$ satisfies the Ward identity
$\eta^{\prime}_{u}+2\eta^{\prime}=\varepsilon$. A consequence of the equality
$\eta=\eta^{\prime}$ is that $\eta_{u}=\eta^{\prime}_{u}$ as well, which is
verified here to leading order in an expansion in $\varepsilon=4-D$.

Distinct critical exponents appear, however, when we perturb away from the
buckling threshold, with one relevant direction that flows away from the fixed
point CvKth. Upon linearizing about this fixed point, we obtain and
diagonalize the resulting Jacobian matrix to produce the following eigenvalues
valid to $\mathcal{O}(\varepsilon)$,

$y_{0}=2-\dfrac{13\varepsilon}{25}\;,\
y_{1}=y_{2}=-\dfrac{\varepsilon}{25}\;,\ y_{3}=-\varepsilon\;.$ (8.18)

We have three irrelevant directions with negative eigenvalues ($y_{1,2,3}<0$)
and one relevant direction with a positive eigenvalue ($y_{0}>0$). If we write
$K(s)=K_{*}+\delta K(s)$, where $K_{*}=1+\varepsilon/50$ is the fixed point
value and $\delta K(s)$ is a small deviation, then we find that $\delta
K(s)\simeq\delta K(0)e^{y_{0}s}$ grows with scale as the renormalization group
flow proceeds away from CvKth along the outgoing separatrix. Note that $\delta
K(0)\propto(\epsilon_{c}(T)-\epsilon)$ is controlled by the distance to the
buckling transition. This relation can be easily obtained by expanding
$\sigma/(\kappa\Lambda^{2})=K^{-1}-1$ to linear order in $\delta K$ and
setting $\sigma=B\epsilon$ as appropriate in the isometric ensemble. Upon
using standard renormalization group arguments, we obtain the divergent
correlation length to be

$\xi\propto|\epsilon-\epsilon_{c}(T)|^{-\nu^{\prime}}\;,\quad\nu^{\prime}=\dfrac{1}{y_{0}}=\dfrac{1}{2}+\dfrac{13}{100}\varepsilon+\mathcal{O}(\varepsilon^{2})\;.$
(8.19)

Strikingly, we obtain a distinct value of $\nu^{\prime}$ here as compared to
the value of $\nu$ obtained in the isotensional ensemble. In Sec. 4, we in
fact demonstrate through general scaling arguments that
$\nu<2/D<\nu^{\prime}$, which is satisfied to leading order in $\varepsilon$
in our systematic expansion (see also Table 1). While, we don’t expect this
leading order perturbative calculation to remain numerically accurate for
$\varepsilon=2$, the exponent inequality continues to hold in $D=2$, where it
reduces to $\nu<1<\nu^{\prime}$. Note that, the requirement $\nu^{\prime}>1$
for $D=2$ in the isometric ensemble implies an exceptionally strong divergence
for a correlation length in critical phenomena.

The difference in exponents ($\nu\neq\nu^{\prime}$) directly demonstrates that
the universality class for buckling within the isometric ensemble is distinct
from that in the isotensional ensemble, as advertised in the introduction.

Our analysis of the renormalization group flow at the two fixed points
associated with buckling in the two dual ensembles is now complete. While the
calculation presented here was performed within a systematic $\varepsilon=4-D$
expansion at fixed manifold codimension $d_{c}$, we provide the general
results for arbitrary $d_{c}$ in Appendix C. We also present a simpler, yet
uncontrolled, and hence inaccurate one-loop approximation for the fixed points
and exponents directly in fixed internal and embedding dimension ($D=2,d=3$)
in Appendix D. Although, we directly compute only the scaling exponents
associated with the fluctuation spectra ($\eta,\eta^{\prime}$ or equivalently
$\eta_{u},\eta^{\prime}_{u}$) and the correlation length ($\nu,\nu^{\prime}$),
the other exponents defined in Sec. 4 can be obtained through various exponent
identities derived below. All the exponents for both ensembles are listed to
leading order in an $\varepsilon=4-D$ expansion for arbitrary codimension
$d_{c}=d-D$ in Table 1. We also use the most accurate estimates for the $\eta$
exponent in the physically relevant dimensions of $D=2$ and $d=3$, obtained
through the self-consistent screening approximation Le Doussal and Radzihovsky
(1992) along with the scaling relations derived in Sec. 9 to directly quote
our best estimates for the various buckling exponents in both ensembles, in
Table 1.

Below, we present a general scaling theory valid near the buckling transition
and derive relations between various exponents, which acquire nonstandard
forms due to the unusual size dependence exhibited by the buckling transition.
This framework will also allow us to explicitly demonstrate that the
distinction between the isotensional and isometric buckling universality
classes constitutes a mechanical variant of Fisher renormalization Fisher
(1968).

##  9 Scaling Relations & Fisher Renormalization

In this section, we continue to work in the general setting of a
$D$-dimensional elastic manifold embedded in $d$-dimensional Euclidean space.
As before, the codimension $d_{c}=d-D>0$ counts the number of directions in
which the elastic material can deform extrinsically, i.e., the flexural modes.
Close to the thermalized buckling transition, we have universal scaling laws,
just as in conventional critical phenomena, compactly captured by the scaling
form of the free energy itself. Standard renormalization group arguments show
that the free energy density defined by
$F=-(k_{B}T/V_{D})\ln[\int\mathcal{D}h\;e^{-\mathcal{F}/k_{B}T}]$ ($V_{D}$ is
the $D$-dimensional volume, which is just $V_{2}=A$ the area for $D=2$) has a
singular part $F_{s}$ which obeys the following scaling relation close to the
transition Goldenfeld (2018)

$F_{s}=b^{-D}\tilde{\Psi}_{F}\left(\Delta\sigma
b^{1/\nu},\mathcal{E}b^{y_{\mathcal{E}}}\right)\;,\\\ $ (9.1)

where $b$ is a scale factor and $\tilde{\Psi}_{F}$ is a scaling function that
implicitly depends on the system size via $R/b$, the bending rigidity via
$\kappa b^{-\eta}$ and the elastic moduli via $\\{Y,B\\}b^{\eta_{u}}$. We
suppress this dependence to ease notation, but these quantities are important
as they give rise to nonstandard scaling relations later. Eq. 9.1 allows us to
map the physics near the buckling transition onto the mean field theory
derived in Sec. 5. A finite external field $\mathcal{E}$ is a strongly
relevant perturbation, and has been retained with its scaling exponent
$y_{\mathcal{E}}>0$. For convenience, we will work within the isotensional
ensemble where the distance from the buckling threshold is given by
$\Delta\sigma=\sigma_{0}-\sigma_{c}(T)$ 444Note that $\sigma_{c}(T)$ is in
general size dependent in a finite sheet, due to the zero temperature buckling
threshold $\propto 1/R^{2}$, see Sec. 5. In addition, fluctuation effects will
lead to a shift in the buckling threshold $\propto 1/R^{1/\nu}$ as per
standard finite size scaling arguments Amit and Martin-Mayor (2005), which we
expect to dominate over the zero temperature threshold in a large sheet. A
similar argument holds in the isometric ensemble as well.. Equivalent results
for the isometric ensemble will be quoted directly as they follow immediately
by replacing $\Delta\sigma$ with $B\Delta\epsilon=B(\epsilon-\epsilon_{c}(T))$
and exchanging the unprimed exponents for the primed ones. This connection
holds for all the exponent identities, except for the stress-strain exponents
$\theta$ and $\theta^{\prime}$, which require a minor modification due to
their definition (Eq. 4.4), as will be clear later on.

By choosing $b=|\Delta\sigma|^{-\nu}\propto\xi$, we scale out the
$\Delta\sigma$ dependence to obtain

$F_{s}=|\Delta\sigma|^{\nu
D}\Psi_{F}\left(\dfrac{\mathcal{E}}{|\Delta\sigma|^{\phi}}\right)\;,$ (9.2)

where $\Psi_{F}(x)=\tilde{\Psi}_{F}(1,x)$. The crossover or gap exponent is
given by $\phi=\nu y_{\mathcal{E}}$ (correspondingly
$\phi^{\prime}=\nu^{\prime}y^{\prime}_{\mathcal{E}}$ in the isometric
ensemble). The full scaling form,

$F_{s}=|\Delta\sigma|^{\nu
D}\Psi_{F}\left(\dfrac{\mathcal{E}}{|\Delta\sigma|^{\phi}},\kappa|\Delta\sigma|^{\nu\eta},\dfrac{\\{Y,B\\}}{|\Delta\sigma|^{\nu\eta_{u}}},R|\Delta\sigma|^{\nu}\right)\;,$
(9.3)

is a function of five distinct variables, four of which we have suppressed in
Eq. 9.2.

The height field has a scaling dimension $\zeta$. Right at buckling, we expect
$\langle h({\bf r})^{2}\rangle\sim\int_{1/R}\mathrm{d}^{D}q/(\kappa_{R}({\bf
q})q^{4})\sim R^{2\zeta}$, which gives Aronovitz and Lubensky (1988); Guitter
_et al._ (1988, 1989); Aronovitz  _et al._ (1989); Le Doussal and Radzihovsky
(1992)

$\zeta=\dfrac{4-D-\eta}{2}\;,\quad\zeta^{\prime}=\dfrac{4-D-\eta^{\prime}}{2}\;.$
(9.4)

Similarly, the requirement that the rotationally invariant nonlinear strain
tensor $u_{ij}$ renormalize correctly leads to a Ward identity exponent
relation Guitter  _et al._ (1988, 1989); Aronovitz and Lubensky (1988);
Aronovitz  _et al._ (1989); Le Doussal and Radzihovsky (1992), now stated in
general $D$

$2\eta+\eta_{u}=4-D\;,\quad 2\eta^{\prime}+\eta^{\prime}_{u}=4-D\;.$ (9.5)

These are well known identities, which we will use below in deriving
additional exponent relations.

Similar to the mean field treatment of buckled ribbons in Ref. Hanakata  _et
al._ (2021), the average height $\langle h\rangle$ serves as an order
parameter for the buckling transition here. By definition, we obtain (for
general $\mathcal{E}$)

$\langle h\rangle=-\dfrac{\partial
F_{s}}{\partial\mathcal{E}}=|\Delta\sigma|^{\nu
D-\phi}\Psi_{h}\left(\dfrac{\mathcal{E}}{|\Delta\sigma|^{\phi}}\right)$ (9.6)

where $\Psi_{h}(x)=\Psi^{\prime}_{F}(x)$. Now, we identify $\langle
h\rangle\sim b^{\zeta}\propto|\Delta\sigma|^{-\nu\zeta}$, which gives the
exponent relations

$\phi=\dfrac{\nu}{2}(4+D-\eta)\;,\quad\phi^{\prime}=\dfrac{\nu^{\prime}}{2}(4+D-\eta^{\prime})\;.$
(9.7)

As a consistency check, one can easily confirm that this relation for the gap
exponents is equivalent to demanding a trivial renormalization of the external
field $\mathcal{E}$ (see Eq. 7.6 in Sec. 7 for the $D=2$ version).

In the zero field limit ($\mathcal{E}=0$), for strong compression, we
spontaneously break symmetry by buckling, which leads to a finite $\langle
h\rangle$. As $\Delta\sigma\to 0$ in a finite system, the rescaled bulk and
Young’s moduli ($B$ and $Y$) diverge as $|\Delta\sigma|^{-\eta_{u}\nu}$, while
the rescaled system size ($R|\Delta\sigma|^{\nu}$) and bending rigidity
($\kappa|\Delta\sigma|^{\nu\eta}$) become vanishingly small. But the latter
can’t be set to zero quite yet. We know from our mean field analysis that upon
buckling, $\langle h\rangle\propto R/Y^{1/2}$ (Eq. A.4 in Sec. 5). As a
result, the elastic moduli act as relevant scaling variables at the transition
while the system size behaves as a _dangerously irrelevant variable_ Amit and
Peliti (1982) whose scaling nontrivially affects the critical exponents. For
$\mathcal{E}=0$, this observation leads to

$\lim_{\Delta\sigma\to
0}\Psi_{h}(0)\propto\dfrac{R}{Y^{1/2}}|\Delta\sigma|^{\nu(1+\eta_{u}/2)}\;.\\\
$ (9.8)

Note that this interesting size dependence for the scaling function is not a
statement of conventional finite-size scaling Amit and Martin-Mayor (2005); it
is instead a unique feature arising from the long wavelength nature of the
buckling transition. By combining Eqs. 9.4, 9.5 with this asymptotic
behaviour, we can compute the order parameter exponents $\beta,\beta^{\prime}$
for the buckled height (Eq. 4.1) in our two ensembles to be

$\beta=\nu\left(1-\dfrac{\eta}{2}\right)\;,\quad\beta^{\prime}=\nu^{\prime}\left(1-\dfrac{\eta^{\prime}}{2}\right)\;.$
(9.9)

This exponent identity is new and distinct from the usual hyperscaling
relation that relates $\beta$ and $\eta$ in conventional critical phenomena
(the latter in our current notation would read as $\beta=-\nu\zeta<0$, which
is obviously wrong).

We can similarly compute the susceptibility exponents $\gamma,\gamma^{\prime}$
defined in Sec. 4. From our mean field analysis, we know that $\chi\propto
R^{2}$ (Eq. 5.7). Upon taking a derivative of Eq. 9.6 and evaluating it at
$\mathcal{E}=0$, we obtain

$\lim_{\Delta\sigma\to 0}\Psi^{\prime}_{h}(0)\propto
R^{2}|\Delta\sigma|^{2\nu}\;.\\\ $ (9.10)

The susceptibility scaling exponents then satisfy

$\gamma=\nu(2-\eta)\;,\quad\gamma^{\prime}=\nu^{\prime}(2-\eta^{\prime})\;.$
(9.11)

Although we recover the standard Fisher’s identity, its appearance is in fact
nontrivial, as is easily seen by noting that Fisher’s identity reflects the
equilibrium fluctuation-response relation Chaikin  _et al._ (1995),

$k_{B}T\chi=\int\mathrm{d}^{D}r\langle h({\bf r})h({\bf 0})\rangle\;,$ (9.12)

which by a naïve application of scaling would give $\chi\sim\xi^{4-\eta}$
resulting in $\gamma=\nu(4-\eta)$ instead of Eq. 9.11. The resolution, as
before, lies in the nontrivial size dependence of the correlation function
integral (because $h({\bf r})\sim R$), which leads to
$k_{B}T\chi\sim\xi^{4-\eta}(R/\xi)^{2}\propto\xi^{2-\eta}$ and then correctly
producing Eq. 9.11. In addition, combining Eqs. 9.9 and 9.11 results in
another unusual exponent identity

$\gamma=2\beta\;,\quad\gamma^{\prime}=2\beta^{\prime}\;.$ (9.13)

Note that the inequality $\nu\neq\nu^{\prime}$ (discussed below) will
necessarily imply additional exponent differences between the two ensembles,
such as $\phi\neq\phi^{\prime}$, $\beta\neq\beta^{\prime}$ and
$\gamma\neq\gamma^{\prime}$ (see Table 1 for a summary).

Next we move on to the nonlinear response in the presence of an external
field. For finite $\mathcal{E}$ we now take the limit $\Delta\sigma\to 0$,
which requires that we focus on the $x\to\infty$ asymptotics of
$\Psi_{h}(x)\sim x^{1/\delta}$. Once again, we must be careful with regard to
scaling with the elastic moduli and the system size. The mean field analysis
(Eq. 5.6) dictates that the prefactor to the nonlinear response itself scales
as $(R^{4}/Y)^{1/3}$. Upon taking this dependence into account, we have

$\lim_{\Delta\sigma\to
0}\Psi_{h}\left(\dfrac{\mathcal{E}}{|\Delta\sigma|^{\phi}}\right)\propto\left(\dfrac{\mathcal{E}}{|\Delta\sigma|^{\phi}}\right)^{1/\delta}\left(\dfrac{R^{4}|\Delta\sigma|^{4\nu}}{Y|\Delta\sigma|^{-\nu\eta_{u}}}\right)^{1/3}\;,\\\
$ (9.14)

By demanding that the $\Delta\sigma$ dependence cancel as $\Delta\sigma\to 0$,
we determine $\delta$ and $\delta^{\prime}$. Remarkably, upon using the other
identites from Eqs. 9.4, 9.5 and 9.7, we obtain

$\delta=\delta^{\prime}=3\;,$ (9.15)

which is an exact result, independent of dimension! In Appendix E, we show how
the above results can all be combined into a simpler scaling form for $\langle
h\rangle$ with a _single_ size dependent scaling variable and a modified gap
exponent.

Exponent | Relation | Isotensional ensemble | Isometric ensemble
---|---|---|---
$D=2,\,d=3^{\dagger}$ | $D=4-\varepsilon,\,d=D+d_{c}$ | $D=2,\,d=3^{\dagger}$ | $D=4-\varepsilon,\,d=D+d_{c}$
$\eta,\,\eta^{\prime}$ | $\eta^{\prime}=\eta$ | $0.821$ | $\dfrac{12\varepsilon}{24+d_{c}}$ | $0.821$ | $\dfrac{12\varepsilon}{24+d_{c}}$
$\eta_{u},\eta^{\prime}_{u}$ | $\begin{aligned} 2\eta+\eta_{u}&=4-D\\\ 2\eta^{\prime}+\eta^{\prime}_{u}&=4-D\end{aligned}$ | $0.358$ | $\dfrac{d_{c}\varepsilon}{24+d_{c}}$ | $0.358$ | $\dfrac{d_{c}\varepsilon}{24+d_{c}}$
$\nu,\nu^{\prime}$ | $\begin{aligned} \nu&=1/(2-\eta)\\\ \nu^{\prime}&=\nu/(D\nu-1)\end{aligned}$ | $0.848$ | $\dfrac{1}{2}+\dfrac{3\varepsilon}{24+d_{c}}$ | $1.218$ | $\dfrac{1}{2}+\dfrac{(12+d_{c})\varepsilon}{(96+4d_{c})}$
$\beta,\beta^{\prime}$ | $\begin{aligned} \beta&=\nu(1-\eta/2)\\\ \beta^{\prime}&=\nu^{\prime}(1-\eta^{\prime}/2)\end{aligned}$ | $\dfrac{1}{2}$ | $\dfrac{1}{2}$ | $0.718$ | $\dfrac{1}{2}+\dfrac{d_{c}\varepsilon}{96+4d_{c}}$
$\gamma,\gamma^{\prime}$ | $\begin{aligned} \gamma&=\nu(2-\eta)=2\beta\\\ \gamma^{\prime}&=\nu^{\prime}(2-\eta^{\prime})=2\beta^{\prime}\end{aligned}$ | $1$ | $1$ | $1.436$ | $1+\dfrac{d_{c}\varepsilon}{48+2d_{c}}$
$\delta,\delta^{\prime}$ | $\begin{aligned} \delta\beta&=\gamma+\beta\\\ \delta^{\prime}\beta^{\prime}&=\gamma^{\prime}+\beta^{\prime}\end{aligned}$ | $3$ | $3$ | $3$ | $3$
$\theta,\theta^{\prime}$ | $\begin{aligned} &\theta=1/(\nu D-1)\\\ \theta^{\prime}&=\nu^{\prime}D-1,\,\theta^{\prime}=\theta\end{aligned}$ | $1.436$ | $1+\dfrac{d_{c}\varepsilon}{48+2d_{c}}$ | $1.436$ | $1+\dfrac{d_{c}\varepsilon}{48+2d_{c}}$
$\phi,\phi^{\prime}$ | $\begin{aligned} \phi&=\nu(4+D-\eta)/2\\\ \phi^{\prime}&=\nu^{\prime}(4+D-\eta^{\prime})/2\end{aligned}$ | $2.196$ | $2+\dfrac{(12-d_{c})\varepsilon}{96+4d_{c}}$ | $3.155$ | $2+\dfrac{3(4+d_{c})\varepsilon}{4(24+d_{c})}$

† Exponents computed using $\eta=4/(1+\sqrt{15})$ obtained from the self-
consistent screening approximation Le Doussal and Radzihovsky (1992).

Table 1: The two buckling universality classes. We list all the scaling
exponents along with the relevant scaling identites for the isotensional
(unprimed exponents) and isometric (primed exponents) ensembles. The exponents
are obtained within an $\varepsilon$-expansion ($D=4-\varepsilon,d=D+d_{c}$),
accurate to $\mathcal{O}(\varepsilon)$ and also in physical dimensions
($D=2,d=3$) using the best self-consistent estimates for $\eta$ Le Doussal and
Radzihovsky (1992). Only $\beta=1/2$, $\gamma=1$ (isotensional) and
$\delta=\delta^{\prime}=3$ (both ensembles) are exact to all orders and
independent of dimensionality.

Finally, we study the anomalous stress-strain curve exponents
$\theta,\theta^{\prime}$ defined by Eqs. 4.4. These quantities are simpler as
they approach finite limits when the system size $R\to\infty$. Here we
distinguish the isotensional and the isometric ensemble as $\theta$ and
$\theta^{\prime}$ are defined differently in the two. We now note that the
scaling variable is $B\Delta\epsilon$ in the isometric ensemble, set
$\mathcal{E}=0$ in Eq. 9.2 and employ the definitions in Eqs. 3.6, 3.9 to get

$\displaystyle\langle\epsilon\rangle$ $\displaystyle=-\dfrac{\partial
F_{s,\sigma}}{\partial\sigma_{0}}\propto|\Delta\sigma|^{\nu
D-1}\sim|\Delta\sigma|^{1/\theta}\quad(\textrm{isotensional})\;,$ (9.16a)
$\displaystyle\dfrac{\langle\sigma\rangle}{B}$ $\displaystyle=\dfrac{\partial
F_{s,\epsilon}}{\partial(B\epsilon)}\propto|B\Delta\epsilon|^{\nu^{\prime}D-1}\sim|\Delta\epsilon|^{\theta^{\prime}}\quad(\textrm{isometric})\;.$
(9.16b)

In this mechanical context $\theta,\theta^{\prime}$ take on the role usually
played by energy scaling in conventional critical phenomena. We now obtain
exponent relations analogous to Josephson’s hyperscaling relation Goldenfeld
(2018)

$\theta=\dfrac{1}{\nu D-1}\;,\quad\theta^{\prime}=\nu^{\prime}D-1\;.$ (9.17)

We shall see later that $\theta,\theta^{\prime}>1$ (see Table 1), which leads
to a crucial distinction between the two ensembles. In the isotensional
ensemble, for $1/\theta<1$, the anomalous sublinear response dominates any
linear Hookean response as $\Delta\sigma\to 0$ Košmrlj and Nelson (2016). In
contrast, in the isometric ensemble, for $\theta^{\prime}>1$, the dominant
strain response is in fact the nonsingular linear term as $\Delta\epsilon\to
0$. This dichotomy reflects a crucial physical consequence of the different
boundary conditions: in the isotensional ensemble, the sheet is infinitely
compliant to homogeneous dilations or contractions in the plane, which is a
zero mode of the system, but in the isometric ensemble, the clamped boundary
conditions prohibit this zero mode and the sheet has a finite compliance to
homogeneous isotropic distortions.

We now discuss one more scaling relation that is _only_ true in the
isotensional ensemble and _not_ in the isometric ensemble. In the isotensional
case, as we saw before in Sec. 7, the in-plane tension $\sigma=\sigma_{0}$
does not receive any graphical corrections as $\bar{v}=0$ identically. This is
true in any dimension and to all orders in perturbation theory, as a
consequence of the fact that $\sigma_{0}$ is the sole term that breaks
rotational invariance, while the bending and nonlinear stretching terms
preserve rotational symmetry. This nonrenormalization condition then implies

$\nu=\dfrac{1}{2-\eta}\;,$ (9.18)

as we already saw by explicit calculation in Eq. 8.10. Because we must have
$\eta,\eta_{u}>0$, Eq. 9.5 implies that $\eta\leq(4-D)/2$ (for $D\leq 4$).
This inequality then shows that $\nu\leq 2/D$ always. An additional important
consequence of Eq. 9.18 is that both $\beta$ and $\gamma$ take on their mean-
field values,

$\beta=\dfrac{1}{2}\;,\quad\gamma=1\;,$ (9.19)

in any dimension. Note that such relations do _not_ hold in the isometric
ensemble as can be seen in an $\varepsilon$-expansion as displayed in Table 1.
Eq. 9.18 determining $\nu$ when plugged into the hyperscaling relation (Eq.
9.17) also gives $\theta=(2-\eta)/(D-2+\eta)$, which is consistent with the
scaling expected from
$\langle\epsilon\rangle\sim\int(\mathrm{d}^{D}r/V_{D})\langle|\bm{\nabla}h|^{2}\rangle\sim\xi^{2\zeta-2}$.
This result agrees with previously derived scaling relations in arbitrary
dimensions Guitter  _et al._ (1989); Aronovitz and Lubensky (1988) and leads
to $\theta=(2-\eta)/\eta$ when $D=2$ Košmrlj and Nelson (2016).

This completes the derivation of the various scaling exponents in the two
ensembles. The values of the exponents computed within an
$\varepsilon$-expansion ($D=4-\varepsilon,d=D+d_{c}$)and by using estimates
from a self-consistent calculation Le Doussal and Radzihovsky (1992) in $D=2$
and $d=3$ dimensions are displayed in Table 1. With these scaling identites in
hand, we can finally address the last key result of the paper, which is to
show that the two ensembles are related to each other via a mechanical analog
of Fisher renormalization Fisher (1968). As the isotensional and isometric
ensembles are thermodynamic duals of each other, the corresponding free energy
densities are related to each other via a Legendre transformation (in the
thermodynamic limit $V_{D}\to\infty$)

$F_{\epsilon}=\min_{\sigma_{0}}\left(F_{\sigma}+\sigma_{0}\epsilon\right)=F_{\sigma}(\sigma_{*})+\sigma_{*}\epsilon\;,$
(9.20)

where $\sigma_{*}$ solves $\epsilon=-\partial
F_{\sigma}/\partial\sigma_{0}|_{\sigma_{0}=\sigma_{*}}$, and the free energy
densities are given by
$F_{\sigma,\epsilon}=-(k_{B}T/V_{D})\ln\mathcal{Z}_{\sigma,\epsilon}$. Now,
near the buckling transition (at zero symmetry-breaking field
$\mathcal{E}=0$), we can use the scaling theory developed above to obtain

$\Delta\epsilon\sim|\sigma_{*}-\sigma_{c}|^{1/\theta}\implies\sigma_{*}-\sigma_{c}\sim|\Delta\epsilon|^{\theta}\;,$
(9.21)

where we have assumed $\theta>1$ and retained only the leading order term as
$\Delta\epsilon\to 0$. Upon combining Eq. 9.21 with the scaling of the
singular part of the free energy densities
$F_{\epsilon}\sim|\Delta\epsilon|^{1+\theta^{\prime}}$ and
$F_{\sigma}\sim|\Delta\sigma|^{1+1/\theta}$ 555The free energy densities also
have nonsingular terms that need to be properly accounted for. In the limit
$\Delta\sigma\to 0$ and $\Delta\epsilon\to 0$, we have the following
expansions for the respective free energy densities:
$\mathcal{F}_{\sigma}=a_{0}+a_{1}\Delta\sigma+a_{2}|\Delta\sigma|^{1+1/\theta}+\mathcal{O}(\Delta\sigma^{2})$
and
$\mathcal{F}_{\epsilon}=b_{0}+b_{1}\Delta\epsilon+b_{2}|\Delta\epsilon|^{1+\theta^{\prime}}+\mathcal{O}(\Delta\epsilon^{2})$,
where $a_{0,1,2}$ and $b_{0,1,2}$ are constants., we require that both sides
of Eq. 9.20 scale in the same way as $\Delta\sigma\to 0$. This constraint
gives the equality

$\theta=\theta^{\prime}\;.$ (9.22)

Although $\theta=\theta^{\prime}$, from Eq. 9.17 we immediately see that the
correlation length exponents must differ in the two ensembles,
$\nu\neq\nu^{\prime}$. The simple form of the nontrivial relation in Eq. 9.22
reflects the definition of $\theta,\theta^{\prime}$ in Eq. 4.4. Eqs. 9.20,
9.21 now lead to

$\displaystyle\beta^{\prime}=\beta\theta\;,\quad\gamma^{\prime}=\gamma\theta\;,$
(9.23a)
$\displaystyle\phi^{\prime}=\phi\theta\;,\quad\nu^{\prime}=\nu\theta\;.$
(9.23b)

The last of these relations can be solved using Eq. 9.17 to explicitly give
the important connection $\nu^{\prime}=\nu/(\nu D-1)$. With the help of Eq.
9.18, this relation simplifies to

$\nu^{\prime}=\dfrac{1}{D-2+\eta}\;,$ (9.24)

which was obtained previously Guitter  _et al._ (1989), without however
recognizing a possible distinction in ensembles. In $D=2$, we get
$\nu^{\prime}=1/\eta$, which upon noting that $\eta<1$ leads to
$\nu^{\prime}>1$, which has been observed in old Monte-Carlo simulations of
thermalized buckling in clamped sheets Guitter  _et al._ (1990). In general
$D$, we can show that the isotensional and isometric ensembles have differing
correlation length exponents such that

$\nu<\dfrac{2}{D}<\nu^{\prime}\;.$ (9.25)

which is consistent with our renormalization group results in
$D=4-\varepsilon$ dimensions. Finally, demanding that Eq. 9.23 be consistent
with Eqs. 9.9, 9.11 leads to equality of the eta exponents in the two
ensembles

$\eta^{\prime}=\eta\;,\quad\eta^{\prime}_{u}=\eta_{u}\;,$ (9.26)

the latter being a consequence of the Ward identity (Eq. 9.5).

We note that the difference in some of the exponents between the ensembles is
not simply because the control variables ($\Delta\epsilon$ and $\Delta\sigma$)
have different dimensions. In fact, both the scaling variables we use have
dimensions of stress and are equivalent to each other: we use $\Delta\sigma$
in the isotensional ensemble and $B\Delta\epsilon$ (not $\Delta\epsilon$) in
the isometric ensemble. Hence the difference in exponents between the
ensembles is due to a genuine change in the fixed point and its associated
universality class, as confirmed by our renormalization group calculations.

These results are reminiscent of the Fisher renormalization of critical
exponents due to hidden variables Fisher (1968), and are also related to the
problems of a constrained Rudnick  _et al._ (1974) or a compressible magnets
Sak (1974), where the presence of a constraint (much like Eq. 9.20) leads to
modified exponents. In conventional critical phenomena, such as in 3D magnets
or superfluid He, Fisher renormalization doesn’t affect the numerical values
of exponents by much as it usually involves dividing the conventional
exponents by $1-\alpha$, where $\alpha$ is the specific heat exponent, which
is often a rather small correction Fisher (1968). Here, however, the exponent
$\theta$ replaces $1-\alpha$ in the mechanical context, allowing for a much
stronger distinction in critical behaviour between the two ensembles.

In Table 1, we see that to leading order in an $\varepsilon$-expansion, all
the equalities in Eqs. 9.22 and 9.23 are satisfied. Within the
$\varepsilon$-expansion, we find that the anomalous exponents are equal to
leading order in the two ensembles, i.e., $\eta=\eta^{\prime}$ and
$\eta_{u}=\eta^{\prime}_{u}$, as expected from our scaling considerations (Eq.
9.26). All the exponents in both ensembles are recapitulated in Table 1 along
with their exponent identities.

##  10 Discussion

Although the study of thermalized membranes is more than three decades old
Nelson  _et al._ (2004), it has been revitalized in recent years by enhanced
interest in 2D materials such as graphene and MoS2. Motivated by the ability
to study extreme mechanics in such ultrathin materials Blees  _et al._ (2015),
we have investigated the impact of thermal fluctuations on a classic (circa
1757!) Euler buckling instability of thin plates. By viewing the finite
temperature buckling transition through the lens of critical phenomena, we
have uncovered new exponent relations and remarkable phenomena that tie
together geometry, mechanics and fluctuations in a thin elastic sheet.

Near a thermodynamic continuous phase transition, fluctuations emerge on all
scales and physics becomes universal. As we have shown, a similar situation
arises on the verge of a mechanical instability, such as buckling, though with
some surprises. The long-wavelength nature of the buckling transition leads to
unusual critical scaling behaviour reflected in the system size dependence of
the mechanical response. Additionally, buckling can be actuated under either
isotensional (constant stress) or isometric (constant strain) loading, which
as we have found, actually constitute separate universality classes. This
remarkable feature highlights the importance of oft-neglected boundary
conditions that when clamped can induce a novel thermally generated
spontaneous tension and modify important scaling exponents. Our work
demonstrates that the inequivalence of mechanical ensembles distinguished by
their boundary conditions exemplifies the phenomenon of Fisher renormalization
Fisher (1968) in a mechanical context.

We emphasize the salient role of geometry in isotropic thermalized buckling.
Much of the phenomena discussed here arise due to the inevitable geometric
coupling between in-plane stretching and out of plane bending, ubiquitous in
thin plates but absent in lower dimensional counterparts, such as slender
filaments. As a consequence, there is no analogue of our results in the finite
temperature buckling of beams and polymers Odijk (1998); *hansen1999buckling;
*baczynski2007stretching; *emanuel2007buckling; *bedi2015finite;
*stuij2019stochastic. Single-molecule measurements with polymers have noted an
inequivalence of similar mechanical ensembles Lubensky and Nelson (2002);
Keller  _et al._ (2003); *sinha2005inequivalence, though in this case due to
finite size effects. In contrast for thin sheets, the ensemble inequivalence
at buckling survives the thermodynamic limit, as it instead originates from
the tensionless flat phase being a critical phase, with fluctuations on all
scales.

Our work is directly relevant to recent experiments that have probed the
mechanics of 2D materials in various geometries Blees  _et al._ (2015);
Nicholl  _et al._ (2015); *nicholl2017hidden. Boundary manipulation is a
popular way to induce strong deformations and morphologies Bao  _et al._
(2009); Vandeparre  _et al._ (2011). When coupled to the material’s electronic
properties, this also allows for strain engineering synthetic gauge fields
Pereira  _et al._ (2010); *guinea2010energy. Our results suggest that the
interplay of thermal fluctuations and boundary constraints can be important in
many of these contexts, allowing for enhanced control of the emergent
mechanics in nanoscale systems. Anisotropic buckling is particularly relevant
in such solid-state devices, either in terms of uniaxially compressed ribbons
Košmrlj and Nelson (2016) or sheets crushed in the presence of a background
aligning field Le Doussal and Radzihovsky (2021), both of which pose
challenging directions for future research. It is also known that the
influence of boundaries can often persist in macroscopically large diffuse
regions in slender elastic bodies Schroll  _et al._ (2011), an effect that is
amplified by the geometry of plates Lobkovsky and Witten (1997);
*cerda2004elements; *barois2014curved and shells Mahadevan  _et al._ (2007);
*santangelo2013nambu. It would be interesting to explore the consequence of
thermal fluctuations in these cases, where boundary effects are again
particularly important.

While much of our analysis focused on the buckling transition, far beyond the
threshold, the sheet adopts a curved profile whose mechanical description is
akin to that of thin shells. The ensuing curvature can be manipulated to
control localized deformations Vaziri and Mahadevan (2008);
*evans2017geometric and fluctuation driven nonlinear response Paulose  _et
al._ (2012); *kovsmrlj2017statistical. Such curved geometries offer tunable
mechanisms to modulate the mechanical and vibrational properties of
electromechanical resonators Lifshitz and Cross (2008), another attractive
direction for future research. Postbuckled states also often exhibit
bistability and sudden snap-through transitions that exhibit critical slowing
down even in the absence of fluctuations Gomez  _et al._ (2017). It would be
of interest to extend our analysis to incorporate such hysteritic and dynamic
effects at finite temperature, but this remains a formidable challenge. We
hope that this work will encourage future explorations at the rich
intersection of geometry, statistical mechanics and elastic instabilities.

_Note added_ : We recently became aware of related work by Leo Radzihovsky and
Pierre Le Doussal Le Doussal and Radzihovsky (2021) on the buckling transition
in a thermalized membrane subjected to an external aligning field, which
unlike our work, breaks rotational symmetry explicitly in the bulk.

###### Acknowledgements.

SS acknowledges the support of the Harvard Society of Fellows. We thank Mark
Bowick, Paul Hanakata, Andrej Košmrlj, John Toner and Leo Radzihovsky for
illuminating discussions. This work was also supported by the NSF through the
Harvard Materials Research Science and Engineering Center, via Grant No.
DMR-2011754 as well as by Grant No. DMR-1608501.

##  Appendix A Integrating out the in-plane phonons

As ${\bf u}$ appears quadratically in the Hamiltonian (Eq. 2.1), we can
exactly integrate it out. To do this, we separate the average strain
($u^{0}_{ij}$) from the nonzero wavelength deformations (i.e., ${\bf
q}\neq{\bf 0}$, denoted by the prime on the ${\bf q}$ integral), to write

$\displaystyle u_{ij}({\bf r})=u_{ij}^{0}+\int_{{\bf q}}^{\prime}\;e^{i{\bf
q}\cdot{\bf r}}\left[\dfrac{i}{2}(q_{i}u_{j}({\bf q})+q_{j}u_{i}({\bf
q}))+\mathcal{A}_{ij}({\bf q})\right]\;,$ (A.1)
$\displaystyle\textrm{with}\quad
u^{0}_{ij}=\dfrac{1}{A}\int\mathrm{d}^{2}r\;u_{ij}\;,\quad\mathcal{A}_{ij}({\bf
r})=\dfrac{1}{2}\partial_{i}h\partial_{j}h\;,$ (A.2)

where $\int_{{\bf q}}=\int\mathrm{d}^{2}q/(2\pi)^{2}$. Note that, while there
are only two independent in-plane phonon degrees of freedom ($u_{i}({\bf q})$)
for nonzero wavevector, the homogeneous part of the strain tensor
($u^{0}_{ij}$) has _three_ independent components, corresponding to the three
distinct modes of macroscopically deforming a 2D solid. It is well known that
only the transverse component of $\mathcal{A}_{ij}$ is important, as the rest
can be absorbed into a global translation zero mode (constant displacement)
Nelson  _et al._ (2004). The total Hamiltonian now takes the form
$\mathcal{H}=\mathcal{H}^{\prime}+\mathcal{H}^{0}$

$\displaystyle\mathcal{H}^{\prime}$
$\displaystyle=\int\dfrac{\mathrm{d}^{2}q}{(2\pi)^{2}}\left[\dfrac{\kappa}{2}q^{4}\left|h_{\bf
q}\right|^{2}+\dfrac{1}{2}u_{i}({\bf q})\left(\mu
q^{2}\mathcal{P}^{T}_{ij}+(2\mu+\lambda)q^{2}\mathcal{P}^{L}_{ij}\right)u_{j}(-{\bf
q})\right]-\mathcal{E}h_{{\bf q}={\bf 0}}$
$\displaystyle+\int^{\prime}\dfrac{\mathrm{d}^{2}q}{(2\pi)^{2}}\left[i\mu\left(q_{i}u_{j}({\bf
q})+q_{j}u_{i}(-{\bf q})\right)\mathcal{A}_{ij}(-{\bf q})+i\lambda
q_{i}u_{i}({\bf q})\mathcal{A}_{kk}(-{\bf
q})+\dfrac{1}{2}\left(2\mu|\mathcal{A}_{ij}({\bf
q})|^{2}+\lambda|\mathcal{A}_{kk}({\bf q})|^{2}\right)\right]\;,$ (A.3)
$\displaystyle\mathcal{H}^{0}$
$\displaystyle=\dfrac{A}{2}\left[2\mu(u^{0}_{ij})^{2}+\lambda(u^{0}_{kk})^{2}\right]-A\sigma^{\rm
ext}_{ij}\left(u^{0}_{ij}-\mathcal{A}^{0}_{ij}\right)\;,$ (A.4)

where $\mathcal{H}^{\prime}$ includes all contributions from the ${\bf
q}\neq{\bf 0}$ in-plane phonon modes and $\mathcal{H}^{0}$ includes all the
terms corresponding to the ${\bf q}={\bf 0}$ phonon modes. Here we have used
the transverse ($\mathcal{P}^{T}_{ij}({\bf q})=\delta_{ij}-q_{i}q_{j}/q^{2}$)
and longitudinal ($\mathcal{P}^{L}({\bf q})=q_{i}q_{j}/q^{2}$) projection
operators and written $\mathcal{A}^{0}_{ij}=(1/A)\int\mathrm{d}{\bf
r}\;\mathcal{A}_{ij}({\bf r})$. The ${\bf q}={\bf 0}$ and ${\bf q}\neq{\bf 0}$
in-plane phonon modes clearly decouple from each other, so when we integrate
them out, the total free energy is simply
$\mathcal{F}=\mathcal{F}^{\prime}+\mathcal{F}^{0}$, where
$\mathcal{F}^{\prime}$ arises from integrating out the ${\bf q}\neq{\bf 0}$
phonons and $\mathcal{F}^{0}$ arises from the ${\bf q}={\bf 0}$ phonon modes.
The former is a standard calculation Nelson  _et al._ (2004), which gives

$\mathcal{F}^{\prime}=\int\mathrm{d}^{2}r\left\\{\dfrac{\kappa}{2}(\nabla^{2}h)^{2}+\dfrac{Y}{2}\left(\dfrac{1}{2}\mathcal{P}^{T}_{ij}\partial_{i}h\partial_{j}h\right)^{2}-\mathcal{E}h\right\\}\;.$
(A.5)

Note that this part of the calculation is common to both ensembles. For the
zero mode calculation, we consider the two ensembles separately. In the
isotensional ensemble, $\sigma^{\rm ext}_{ij}=\sigma_{0}\delta_{ij}$, which
gives

$\mathcal{H}^{0}_{\sigma}=\dfrac{A}{2}\left[2\mu(\tilde{u}^{0}_{ij})^{2}+(\mu+\lambda)(u^{0}_{kk})^{2}\right]-A\sigma_{0}\left(u^{0}_{kk}-\mathcal{A}^{0}_{kk}\right)\;,$
(A.6)

where we have decomposed $u^{0}_{ij}$ into its deviatoric part
($\tilde{u}^{0}_{ij}=u_{ij}^{0}-\delta_{ij}u_{kk}^{0}/2$) and its trace
$u^{0}_{kk}$. Both the shear ($\tilde{u}_{ij}^{0}$) and the dilation
($u^{0}_{kk}$) components of the homogeneous strain can be integrated over
freely now to obtain the zero mode contribution to the free energy,

$\mathcal{F}^{0}_{\sigma}=A\sigma_{0}\mathcal{A}^{0}_{kk}=\dfrac{\sigma_{0}}{2}\int\mathrm{d}^{2}r\;|\bm{\nabla}h|^{2}\;,$
(A.7)

in the isotensional ensemble Košmrlj and Nelson (2016). In the isometric
ensemble, we set $\sigma^{\rm ext}_{ij}=0$ and instead have
$(1/A)\int\mathrm{d}{\bf r}\bm{\nabla}\cdot{\bf u}=\epsilon$. The homogeneous
part of the strain tensor is then given by

$u^{0}_{ij}=\tilde{u}^{0}_{ij}+\dfrac{\delta_{ij}}{2}\left(\epsilon+\mathcal{A}^{0}_{kk}\right)\;,$
(A.8)

where we have once again separated out the devaitoric shear component
($\tilde{u}^{0}_{ij}$). We immediately see that, while in the isotensional
ensemble, all three components of $u^{0}_{ij}$ were freely integrated over, in
the isometric ensemble, only _two_ out of the three degrees of freedom can be
freely integrated over. The clamped boundary conditions prevent homogeneous
dilations or contractions, but the two homogeneous shear deformations in
$\tilde{u}^{0}_{ij}$ continue to be zero modes. Upon integrating out
$\tilde{u}^{0}_{ij}$, we obtain

$\displaystyle\mathcal{F}^{0}_{\epsilon}$
$\displaystyle=\dfrac{A}{2}(\mu+\lambda)\left(\epsilon+\mathcal{A}^{0}_{kk}\right)^{2}$
$\displaystyle=\dfrac{B}{2A}\int\mathrm{d}^{2}r\left[\epsilon+\dfrac{1}{2A}\int\mathrm{d}^{2}r^{\prime}|\bm{\nabla}^{\prime}h|^{2}\right]^{2}\;,$
(A.9)

with the bulk modulus $B=\mu+\lambda$, in the isometric ensemble. By adding
together $\mathcal{F}^{\prime}$ from Eq. A.5 with
$\mathcal{F}^{0}_{\sigma,\epsilon}$ (Eqs. A.7 and A.9), we get the total free
energy in the two ensembles (Eqs. 3.5 and 3.7 in the main text).

##  Appendix B Mean field equation of state

For the mean field calculation we use a single mode Galerkin approximation
using $h_{0}({\bf r})=H_{0}J_{0}(q_{n}r)$ for a general wavevector $q_{n}$.
The linear terms are easily diagonalized by $h_{0}({\bf r})$ which is an
eigenfunction of the Laplacian,

$\nabla^{2}h_{0}({\bf r})=-q_{n}^{2}h_{0}({\bf r})\;.$ (B.1)

The nonlinear terms are computed as follows. We first have the integral

$\int\dfrac{\mathrm{d}^{2}r}{A}|\bm{\nabla}h_{0}|^{2}=H_{0}^{2}q_{0}^{2}f(q_{0}R)\;,$
(B.2)

where the dimensionless function is given by

$\displaystyle f(x)$