# Structure of fluctuating thin sheets under random forcing

Chanania Steinbock  Eytan Katzav  Racah Institute of Physics, The Hebrew
University, Jerusalem 9190401, Israel  Arezki Boudaoud  LadHyX, CNRS, Ecole
Polytechnique, Institut Polytechnique de Paris, Palaiseau Cedex 91128 France

###### Abstract

We propose a mathematical model to describe the athermal fluctuations of thin
sheets driven by the type of random driving that might be experienced prior to
weak crumpling. The model is obtained by merging the Föppl-von Kármán
equations from elasticity theory with techniques from out-of-equilibrium
statistical physics to obtain a nonlinear strongly coupled $\phi^{4}$-Langevin
field equation with spatially varying kernel. With the aid of the self-
consistent expansion (SCE), this equation is analytically solved for the
structure factor of a fluctuating sheet. In contrast to previous research
which has suggested that the structure factor follows an anomalous power-law,
we find that the structure factor in fact obeys a logarithmically corrected
rational function. Numerical simulations of our model confirm the accuracy of
our analytical solution.

##  I Introduction

Thin sheets are able to exhibit a remarkable variety of structures with
extreme variation in complexity. For example, the same sheet of paper can be
used to construct a simple greeting card, an intricate work of origami or a
“mere” crumpled ball. Though obtaining the crumpled ball certainly requires
less deliberation, the underlying structure is no less complex than that of
the origami artwork. Indeed, though the structure of a randomly crumpled sheet
has been the subject of active research, both experimentally [1, 2, 3, 4, 5,
6, 7, 8, 9] and through modelling [1, 10, 11, 12], a mathematical theory
derived from first principles remains lacking.

Early work on crumpling focused on membranes or so called “tethered sheets” in
which a two-dimensional triangular lattice of spheres with fixed distances
between neighbouring spheres is randomly fluctuated in space [13, 14, 15, 16,
17, 18]. The great advantage of such models is that they are readily
generalised to include the effects of self-avoidance, an important property in
the context of strong crumpling [10]. An interesting insight that was used in
these early works is a useful analogy with two-dimensional models of
Heisenberg ferromagnetism, where one could think of the normals to the surface
as Heisenberg spins. In this framework, a flat sheet corresponds to a
ferromagnetic state while a crumpled sheet to a paramagnetic state.

The subsequent development of a mathematical theory which included elastic
effects led to a description of the various types of singular structures which
can arise within thin sheets, notably, structures such as ridges and d-cones
[19, 20, 21, 22, 23, 24, 25]. From the theory of elasticity, it is known that
thin sheets undergo two types of deformation – energetically cheap bending and
energetically expensive stretching. The theory of ridges and d-cones describes
from first principles how stretching becomes highly localised in thin sheets
to minimise the energy expenditure and these become the creases of a folded
sheet. This theory, however, has only been successfully applied to
characterise situations with relatively few ridges.

In this paper, we propose a different approach, combining methods from
elasticity theory with out-of-equilibrium statistical physics. From elasticity
theory, we consider the Föppl-von Kármán equations – the theory of linear
elasticity applied to thin sheets while also accounting for the geometric
nonlinear effects of deformation [26]. The result is a pair of coupled
nonlinear partial differential equations which describe the deformation of
thin sheets and though these equations are difficult to solve in general, it
has indeed been shown that ridges and d-cones are particular solutions [20,
21]. By studying a dynamic noisy variant of these equations, we are able to
model a thin sheet subjected to a continually varying athermal random force.
In this manner, we are able to begin characterising the structure of a
randomly fluctuating thin sheet, a structure which could conceivably resemble
to some degree the structure of a crumpled sheet. The underlying logic here
being that it is the pure elastic deformations that set the fundamental modes
of the sheet and these are then plastified into creases by some irreversible
process when the curvature is large. From this perspective, a weakly crumpled
sheet is a nonlinear echo of its elastic fluctuations. Similar arguments are
indeed often invoked in the study of irreversible processes in mechanical
systems. A famous example being the fragmentation of thin rods, described in
[27, 28], where the points that a buckled thin rod snaps are the points of
maximal curvature of the purely elastically deformed rod.

Indeed, the question of characterising the structure of fluctuating thin
sheets is a fundamental one in its own right with ramifications spanning a
diversity of fields [29]. Applications range from the biophysics of cell
membranes [30], to the properties of graphene [31, 32, 33, 34, 35, 36], to the
character of wave turbulence [37, 38, 39, 40, 41, 42, 43, 44, 45]. In recent
years, it has become apparent that the precise structure of fluctuating thin
sheets is important for understanding their mechanical and elastic properties
[46, 47, 48, 49, 50, 51] as well as their acoustic [52, 53] and optical
emissions [54]. Accordingly, we consider the lack of an out-of-equilibrium
model for fluctuating surfaces as a deficiency which we begin to remedy here.

The paper is organised as follows. In Section II, we outline and prepare the
variant of the Föppl-von Kármán equations which we intend to study while in
Section III, we use a method known as the self-consistent expansion (SCE) to
obtain analytical solutions for the structure factor of a fluctuating sheet.
In Section IV, our solution is compared with numerical simulations of the
dynamic Föppl-von Kármán equations and the conclusions are discussed in
Section V. Various technical details of the solution derived in Section III
are postponed to an Appendix.

##  II The Overdamped Dynamic Föppl-von Kármán Equations

The Föppl-von Kármán equations

$\displaystyle P_{ex}$
$\displaystyle=\frac{Eh^{3}}{12\left(1-\nu^{2}\right)}\nabla^{4}\xi$
$\displaystyle-h\left(\frac{\partial^{2}\xi}{\partial
x^{2}}\frac{\partial^{2}\chi}{\partial y^{2}}+\frac{\partial^{2}\xi}{\partial
y^{2}}\frac{\partial^{2}\chi}{\partial x^{2}}-2\frac{\partial^{2}\xi}{\partial
x\partial y}\frac{\partial^{2}\chi}{\partial x\partial y}\right)$ (1)
$\displaystyle 0$
$\displaystyle=\nabla^{4}\chi+E\left[\frac{\partial^{2}\xi}{\partial
x^{2}}\frac{\partial^{2}\xi}{\partial
y^{2}}-\left(\frac{\partial^{2}\xi}{\partial x\partial y}\right)^{2}\right]$
(2)

have long been used to study the equilibrium out-of-plane displacement
$\xi\left(x,y\right)$ of a thin sheet subject to an external pressure
$P_{ex}$. Here, $h$ denotes the thickness of the sheet, $E$ its Young’s
modulus and $\nu$ its Poisson ratio. The scalar field $\chi\left(x,y\right)$
is the Airy stress potential of the deformation. The Monge parameterisation
used, in which the vertical displacement $\xi(x,y)$ is written as a function
of $x$ and $y$, implies that we consider deformations that are mostly flat and
thus we will only focus in the continuation on weak fluctuations. Clearly,
self-avoidance of any real sheet further strengthens this assumption.

To study a driven system, we simply apply Newton’s second law to each element
of the sheet with density $\rho$

$h\rho\frac{\partial^{2}\xi}{\partial
t^{2}}=-P_{ex}+P_{damping}+P_{driving}\,.$ (3)

Though similar equations have been studied in the context of wave turbulence
[37, 38, 39, 40, 41, 42, 43, 44, 45], our approach deviates from this prior
research in that we consider specific damping and driving pressures. In
particular, for the damping, we consider ordinary fluid friction
$P_{damping}=-\alpha\frac{\partial\xi}{\partial t}$ where $\alpha$ denotes the
coefficient of friction while for the driving, we consider zero-mean conserved
Gaussian noise $P_{driving}=\eta\left(\vec{r},t\right)$ with noise amplitude
$D$, that is,

$\displaystyle\left\langle\eta\left(\vec{r},t\right)\right\rangle$
$\displaystyle=0$ (4)
$\displaystyle\left\langle\eta\left(\vec{r},t\right)\eta\left(\vec{r}\,^{\prime},t^{\prime}\right)\right\rangle$
$\displaystyle=-D\delta\left(t-t^{\prime}\right)\nabla^{2}\delta\left(\vec{r}-\vec{r}\,^{\prime}\right)\,.$
(5)

While other driving forces such as thermal fluctuations could certainly be
considered, our aim here is to develop a minimal model for the type of
athermal fluctuations one might expect immediately prior to crumpling. In the
spirit of critical phenomena models that incorporate the symmetries and
conserved quantities of their systems, we impose conserved noise on the sheet
as this ensures that the sheet’s center of mass will not wander in space.
Indeed, any more specific form of noise which is consistent with the
conservation of center of mass could also be considered. Applying the
principle of parsimony however, any more specific form would need additional
justification, for example, by being experimentally motivated.

Finally, in further contrast to the wave turbulence context, we focus on the
overdamped limit in which the inertia term
$h\rho\frac{\partial^{2}\xi}{\partial t^{2}}$ can be neglected relative to the
damping term. Ultimately, this gives the overdamped dynamic Föppl-von Kármán
equation

$\alpha\frac{\partial\xi}{\partial
t}+\frac{Eh^{3}}{12\left(1-\nu^{2}\right)}\nabla^{4}\xi\\\
-h\left(\frac{\partial^{2}\xi}{\partial x^{2}}\frac{\partial^{2}\chi}{\partial
y^{2}}+\frac{\partial^{2}\xi}{\partial y^{2}}\frac{\partial^{2}\chi}{\partial
x^{2}}-2\frac{\partial^{2}\xi}{\partial x\partial
y}\frac{\partial^{2}\chi}{\partial x\partial
y}\right)=\eta\left(\vec{r},t\right)\,.$ (6)

Together with equation (2), this equation describes the stochastic time
varying deformation of a sheet fluctuating under the influence of conserved
Gaussian noise. Alternatively, for a sheet of dimensions $L\times L$,
equations (2) and (6) can be written in Fourier space and combined into the
single equation

$\alpha\frac{\partial\tilde{\xi}_{\vec{n}}}{\partial
t}+\frac{\left(2\pi\right)^{4}}{L^{4}}\frac{Eh^{3}}{12\left(1-\nu^{2}\right)}\left|\vec{n}\right|^{4}\tilde{\xi}_{\vec{n}}\\\
+\frac{\left(2\pi\right)^{4}}{L^{4}}\frac{hE}{2}\sum_{\vec{\ell}_{1}\neq\vec{n}}\sum_{\vec{\ell}_{2}}\sum_{\vec{\ell}_{3}}V_{\vec{n},\vec{\ell}_{1},\vec{\ell}_{2},\vec{\ell}_{3}}\tilde{\xi}_{\vec{\ell}_{1}}\tilde{\xi}_{\vec{\ell}_{2}}\tilde{\xi}_{\vec{\ell}_{3}}=\tilde{\eta}_{\vec{n}}\left(t\right)\,,$
(7)

where the kernel of the sum
$V_{\vec{n},\vec{\ell}_{1},\vec{\ell}_{2},\vec{\ell}_{3}}$, given by

$V_{\vec{n},\vec{\ell}_{1},\vec{\ell}_{2},\vec{\ell}_{3}}=\delta_{\vec{n},\vec{\ell}_{1}+\vec{\ell}_{2}+\vec{\ell}_{3}}\frac{\left|\vec{n}\times\vec{\ell}_{1}\right|^{2}\left|\vec{\ell}_{2}\times\vec{\ell}_{3}\right|^{2}}{\left|\vec{n}-\vec{\ell}_{1}\right|^{4}}\,,$
(8)

is simply the Fourier transform of the transverse projection operator of the
sheet deformation [29, 55], having denoted
$\vec{n}\times\vec{\ell}=n_{x}\ell_{y}-n_{y}\ell_{x}$. Here,
$\tilde{\xi}_{\vec{n}}=\tilde{\xi}_{\vec{n}}(t)$ and
$\tilde{\eta}_{\vec{n}}\left(t\right)$ denote the Fourier components of
$\xi\left(\vec{r},t\right)$ and $\eta\left(\vec{r},t\right)$ respectively,
that is
$\xi\left(\vec{r},t\right)=\sum_{\vec{n}}\tilde{\xi}_{\vec{n}}(t)e^{i\frac{2\pi}{L}\vec{n}\cdot\vec{r}}$
and similarly for $\eta\left(\vec{r},t\right)$. The sums in equation (7) are
carried out over all integer lattice points of $\mathbb{R}^{2}$ excluding the
point $\vec{\ell}_{1}=\vec{n}$. As the nonlinear term of equation (7) is
essentially a quartic interaction, this equation can also be thought of as a
type of Langevin $\phi^{4}$-field theory [56] decorated by a non-trivial
spatially varying kernel
$V_{\vec{n},\vec{\ell}_{1},\vec{\ell}_{2},\vec{\ell}_{3}}$. It is worth noting
that the zeroth mode, $\vec{n}=0$, of equation (7) vanishes and thus
$\tilde{\xi}_{\vec{0}}$ remains constant in time. Accordingly, not only is the
noise conservative (by design) but the deterministic forces are as well,
though this conservation is achieved in a different manner than model B within
the classical classification of Hohenberg and Halperin [57] where a Laplacian
is also applied to the deterministic forces.

It will be convenient to work with a nondimensionalised form of equation (7).
To this end, we can define the dimensionless length, time, vertical
displacement and noise in Fourier space

$\displaystyle\bar{r}$ $\displaystyle=\vec{r}/L$ (9) $\displaystyle\bar{t}$
$\displaystyle=\left[\left(2\pi\right)^{6}hDE/\left(\alpha^{3}L^{8}\right)\right]^{1/2}t$
(10) $\displaystyle\bar{\xi}_{\vec{n}}(\bar{t})$
$\displaystyle=\left[\left(2\pi\right)^{2}\alpha
hE/D\right]^{1/4}\tilde{\xi}_{\vec{n}}(t)$ (11)
$\displaystyle\bar{\eta}_{\vec{n}}\left(\bar{t}\right)$
$\displaystyle=\left[\alpha^{3}L^{16}/\left(\left(2\pi\right)^{10}hD^{3}E\right)\right]^{1/4}\tilde{\eta}_{\vec{n}}\left(t\right)$
(12)

which gives rise to the nondimensionalised equation

$\frac{\partial\bar{\xi}_{\vec{n}}}{\partial\bar{t}}+g\left|\vec{n}\right|^{4}\bar{\xi}_{\vec{n}}\\\
+\frac{1}{2}\sum_{\vec{\ell}_{1}\neq\vec{n}}\sum_{\vec{\ell}_{2}}\sum_{\vec{\ell}_{3}}V_{\vec{n},\vec{\ell}_{1},\vec{\ell}_{2},\vec{\ell}_{3}}\bar{\xi}_{\vec{\ell}_{1}}\bar{\xi}_{\vec{\ell}_{2}}\bar{\xi}_{\vec{\ell}_{3}}=\bar{\eta}_{\vec{n}}\left(\bar{t}\,\right)$
(13)

containing the single dimensionless parameter

$g=\frac{2\pi}{12\left(1-\nu^{2}\right)}\sqrt{\frac{\alpha h^{5}E}{D}}$ (14)

and nondimensionalised zero-mean conserved Gaussian noise with covariance

$\left\langle\bar{\eta}_{\vec{n}}\left(\bar{t}\,\right)\bar{\eta}_{\vec{n}^{\prime}}\left(\bar{t}\,^{\prime}\right)\right\rangle=\left|\vec{n}\right|^{2}\delta_{\vec{n},-\vec{n}^{\prime}}\delta\left(\bar{t}-\bar{t}\,^{\prime}\right)\,.$
(15)

Since the dimensionless parameter $g$ scales like $h^{5/2}$, $g$ will be small
for sufficiently thin sheets. Indeed, a typical sheet of aluminium foil or
steel would have a thickness $h\sim 0.5\mathrm{\,mm}$ and Young’s modulus
$E\sim 200\mathrm{\,GPa}$. The magnitude of the driving force needed to induce
fluctuations has been measured in [40] as being around $D\sim 5\times
10^{-3}\mathrm{\,N^{2}/Hz}$. The frequency dependent frictional damping rate
$\gamma(f)$ of a fluctuating steel sheet has been measured in [42] as being
constant for small frequencies and slowly growing for large frequencies at a
rate of $\gamma\sim f^{0.6}$. Since we intend to describe the large scale
features of a fluctuating sheet and $\gamma$ anyway grows slowly with large
frequencies, little is lost by treating $\gamma$ as constant. From [42], we
thus take $\gamma\sim 5\mathrm{\,Hz}$. Since the coefficient of friction is
simply $\alpha=\rho h\gamma$ where $\rho\sim 8\times 10^{3}\mathrm{\,kg\cdot
m^{-3}}$ denotes the density of the sheet and $\left(1-\nu^{2}\right)\sim 1$,
we obtain under typical circumstances that $g\sim 0.1$ which is indeed
significantly smaller than unity. For even thinner sheets or more violent
driving forces, $g$ will only shrink and thus it is clear that the physically
interesting case is indeed $g\ll 1$.

Equation (13) is a nonlinear Langevin equation for a type of quartic
interaction and thus describes the sheet’s structure both in-equilibrium and
out-of-equilibrium, that is, in principle, equation (13) can be used to obtain
both static moments (equal time) such as the sheet’s structure factor

$S_{\vec{n}}=\left\langle\bar{\xi}_{\vec{n}}(\bar{t}\,)\bar{\xi}_{-\vec{n}}(\bar{t}\,)\right\rangle\,$
(16)

and time-dependent moments such as the time-dependent correlation function

$S_{\vec{n}}(\bar{t}_{1},\bar{t}_{2})=\left\langle\bar{\xi}_{\vec{n}}(\bar{t}_{1})\bar{\xi}_{-\vec{n}}(\bar{t}_{2})\right\rangle\,.$
(17)

Though the methods used in the continuation can be applied to both static and
time-dependent quantities, in the following, we will focus only on determining
the static structure factor and will defer the slightly more involved analysis
of time-dependent quantities to the future.

It is interesting to mention here that based on the analogy to the Heisenberg
model, the static structure factor is also related to the correlations between
normals to the surface. As explained in [55] the tipping angles $\theta(x,y)$
of the normals relative to the $z$ axis obey approximately
$\langle\theta^{2}(x,y)\rangle\simeq\langle|\nabla\xi(x,y)|^{2}\rangle$.

Examining equation (13), one sees that the small parameter $g$ appears in
front of the linear part and thus we find that the fluctuations in the
vertical displacement are strongly coupled to the nonlinear term. On the other
hand, the linear term also grows rapidly with $\vec{n}$ since it is
proportional to $\left|\vec{n}\right|^{4}$. As such, any perturbative
expansion around the linear part of equation (13) which merely treats the non-
linear term as a small correction would be meaningless. At the same time, the
linear term itself can also not simply be neglected. Accordingly, a more
sophisticated approach is required.

##  III The Self-Consistent Expansion (SCE)

The self-consistent expansion (SCE) can be conceived of as a renormalised
perturbation theory [58] capable of providing series approximations even in
the presence of strong coupling. It has been used to solve systems similar to
our own such as the KPZ equation and its variations [59, 60, 61, 62, 63, 64,
65, 66, 67], fracture and wetting fronts [68, 69] and turbulence [70].
Additionally, the SCE is known to provide an extremely successful solution to
the zero-dimensional $\phi^{4}$-theory giving good results at low orders and
providing exact convergence at high orders [71, 72]. Accordingly, it is
reasonable to suppose that the SCE might shed light on our system as well.
Though both simple self-consistent arguments [55] and more sophisticated self-
consistent epsilon expansions in $d=4-\varepsilon$ dimensions [73, 74] have
previously been applied to closely related problems, the approach presented
here has the advantages of being natural, transparent and directly applicable
to the problem at hand.

###  III.1 The Fokker-Planck Equation

To apply the SCE to our system, it is convenient to write its corresponding
Fokker-Planck equation [75]

$\frac{\partial P}{\partial\bar{t}}=\mathcal{O}P$ (18)

where $P=P\left(\left\\{\bar{\xi}_{\vec{n}}\right\\},\bar{t}\,\right)$ denotes
the probability that the system will have Fourier components
$\left\\{\bar{\xi}_{\vec{n}}\right\\}$ at time $\bar{t}$ and $\mathcal{O}$
denotes the Fokker-Planck operator for our system. For equation (13), the
Fokker-Planck equation takes the explicit form

$\frac{\partial
P}{\partial\bar{t}}=\frac{1}{2}\sum_{\vec{n}}\left|\vec{n}\right|^{2}\frac{\partial^{2}P}{\partial\bar{\xi}_{\vec{n}}\partial\bar{\xi}_{-\vec{n}}}+g\sum_{\vec{n}}\left|\vec{n}\right|^{4}\frac{\partial}{\partial\bar{\xi}_{\vec{n}}}\left(\bar{\xi}_{\vec{n}}P\right)\\\
+\frac{1}{2}\sum_{\vec{n}}\frac{\partial}{\partial\bar{\xi}_{\vec{n}}}\left[P\sum_{\vec{\ell}_{1}\neq\vec{n}}\sum_{\vec{\ell}_{2}}\sum_{\vec{\ell}_{3}}V_{\vec{n},\vec{\ell}_{1},\vec{\ell}_{2},\vec{\ell}_{3}}\bar{\xi}_{\vec{\ell}_{1}}\bar{\xi}_{\vec{\ell}_{2}}\bar{\xi}_{\vec{\ell}_{3}}\right]\,,$
(19)

where we have denoted $\bar{\xi}_{\vec{n}}=\bar{\xi}_{\vec{n}}(\bar{t}\,)$ for
the sake of conciseness. Multiplying this equation by any functional
$\mathbb{F}\left(\left\\{\bar{\xi}_{\vec{n}}\right\\}\right)$ of the Fourier
components of the field $\xi$ and integrating over all $\bar{\xi}_{\vec{n}}$
gives after some integration by parts

$\frac{\partial\left\langle\mathbb{F}\right\rangle}{\partial\bar{t}}=\frac{1}{2}\sum_{\vec{n}}\left|\vec{n}\right|^{2}\left\langle\frac{\partial^{2}\mathbb{F}}{\partial\bar{\xi}_{\vec{n}}\partial\bar{\xi}_{-\vec{n}}}\right\rangle-g\sum_{\vec{n}}\left|\vec{n}\right|^{4}\left\langle\frac{\partial\mathbb{F}}{\partial\bar{\xi}_{\vec{n}}}\bar{\xi}_{\vec{n}}\right\rangle\\\
-\frac{1}{2}\sum_{\vec{n}}\sum_{\vec{\ell}_{1}\neq\vec{n}}\sum_{\vec{\ell}_{2}}\sum_{\vec{\ell}_{3}}V_{\vec{n},\vec{\ell}_{1},\vec{\ell}_{2},\vec{\ell}_{3}}\left\langle\frac{\partial\mathbb{F}}{\partial\bar{\xi}_{\vec{n}}}\bar{\xi}_{\vec{\ell}_{1}}\bar{\xi}_{\vec{\ell}_{2}}\bar{\xi}_{\vec{\ell}_{3}}\right\rangle$
(20)

where we have defined the expectation values

$\left\langle\mathbb{F}\right\rangle=\int\prod_{\vec{n}}d\bar{\xi}_{\vec{n}}\,\mathbb{F}\left(\left\\{\bar{\xi}_{\vec{n}}\right\\}\right)P\left(\left\\{\bar{\xi}_{\vec{n}}\right\\},\bar{t}\,\right)\,.$
(21)

In this paper we will only consider static quantities and thus the left-hand
side of equation (20) will vanish. When the function
$\mathbb{F}\left(\left\\{\bar{\xi}_{\vec{n}}\right\\}\right)$ is chosen such
that $\left\langle\mathbb{F}\right\rangle$ is a moment of
$P\left(\left\\{\bar{\xi}_{\vec{n}}\right\\},\bar{t}\,\right)$, equation (20)
becomes an equation relating different moments to each other. For example, for
$\mathbb{F}=\bar{\xi}_{\vec{n}}\bar{\xi}_{\vec{n}^{\prime}}$, equation (20)
reads

$0=\left|\vec{n}\right|^{2}\delta_{\vec{n},-\vec{n}^{\prime}}-g\left(\left|\vec{n}\right|^{4}+\left|\vec{n}^{\prime}\right|^{4}\right)\left\langle\bar{\xi}_{\vec{n}}\bar{\xi}_{\vec{n}^{\prime}}\right\rangle\\\
-\frac{1}{2}\sum_{\vec{\ell}_{2}}\sum_{\vec{\ell}_{3}}\left[\sum_{\vec{\ell}_{1}\neq\vec{n}}V_{\vec{n},\vec{\ell}_{1},\vec{\ell}_{2},\vec{\ell}_{3}}\left\langle\bar{\xi}_{\vec{n}^{\prime}}\bar{\xi}_{\vec{\ell}_{1}}\bar{\xi}_{\vec{\ell}_{2}}\bar{\xi}_{\vec{\ell}_{3}}\right\rangle\right.\\\
\left.+\sum_{\vec{\ell}_{1}\neq\vec{n}^{\prime}}V_{\vec{n}^{\prime},\vec{\ell}_{1},\vec{\ell}_{2},\vec{\ell}_{3}}\left\langle\bar{\xi}_{\vec{n}}\bar{\xi}_{\vec{\ell}_{1}}\bar{\xi}_{\vec{\ell}_{2}}\bar{\xi}_{\vec{\ell}_{3}}\right\rangle\right]\,,$
(22)

which expresses the second moment
$\left\langle\bar{\xi}_{\vec{n}}\bar{\xi}_{\vec{n}^{\prime}}\right\rangle$ in
terms of the fourth moment
$\left\langle\bar{\xi}_{\vec{n}}\bar{\xi}_{\vec{\ell}_{1}}\bar{\xi}_{\vec{\ell}_{2}}\bar{\xi}_{\vec{\ell}_{3}}\right\rangle$.
We could obtain an expression for the fourth moments by subbing
$\mathbb{F}=\bar{\xi}_{\vec{n}_{1}}\bar{\xi}_{\vec{n}_{2}}\bar{\xi}_{\vec{n}_{3}}\bar{\xi}_{\vec{n}_{4}}$
into equation (20) though the resulting expression would need knowledge of the
sixth moment. Reminiscent of the BBGKY hierarchy [76, 77, 78], this property
is generic and deciding how to implement closure to these equations is a
complicated question. Ideally, we would like to argue that the higher order
moments only contribute at a higher order such that they can be neglected at
zeroth order and used to perturbatively correct subsequent orders. As
discussed at the end of section II, the strong coupling of the nonlinear term
guarantees that such a perturbative expansion is destined to fail.

###  III.2 The Self-Consistent Expansion

Instead, the SCE methodology argues that though we cannot neglect the higher
order moments relative to the lower order ones, there might exist some other
linear theory which we can expand around instead. That is, we can define a
linear operator

$\mathcal{O}_{0}P=\frac{1}{2}\sum_{\vec{n}}\left|\vec{n}\right|^{2}\frac{\partial^{2}P}{\partial\bar{\xi}_{\vec{n}}\partial\bar{\xi}_{-\vec{n}}}+\sum_{\vec{n}}\Gamma_{\left|\vec{n}\right|}\frac{\partial}{\partial\bar{\xi}_{\vec{n}}}\left(\bar{\xi}_{\vec{n}}P\right)$
(23)

where $\Gamma_{\left|\vec{n}\right|}$ is a free parameter which will be
determined in the continuation and rewrite equation (18) as

$\frac{\partial
P}{\partial\bar{t}}=\mathcal{O}_{0}P+\left(\mathcal{O}-\mathcal{O}_{0}\right)P\,.$
(24)

At this stage, one can think of $\Gamma_{\left|\vec{n}\right|}$ as an unknown
renormalisation of the the bending rigidity, related to
$\kappa_{R}(\mathbf{q})$ in [55], and its precise significance will become
apparent as the calculation proceeds.

Explicitly in terms of moments, this equation becomes

$0=\frac{1}{2}\sum_{\vec{n}}\left|\vec{n}\right|^{2}\left\langle\frac{\partial^{2}\mathbb{F}}{\partial\bar{\xi}_{\vec{n}}\partial\bar{\xi}_{-\vec{n}}}\right\rangle-\sum_{\vec{n}}\Gamma_{\left|\vec{n}\right|}\left\langle\frac{\partial\mathbb{F}}{\partial\bar{\xi}_{\vec{n}}}\bar{\xi}_{\vec{n}}\right\rangle\\\
-\sum_{\vec{n}}\left(g\left|\vec{n}\right|^{4}-\Gamma_{\left|\vec{n}\right|}\right)\left\langle\frac{\partial\mathbb{F}}{\partial\bar{\xi}_{\vec{n}}}\bar{\xi}_{\vec{n}}\right\rangle\\\
-\frac{1}{2}\sum_{\vec{n}}\sum_{\vec{\ell}_{1}\neq\vec{n}}\sum_{\vec{\ell}_{2}}\sum_{\vec{\ell}_{3}}V_{\vec{n},\vec{\ell}_{1},\vec{\ell}_{2},\vec{\ell}_{3}}\left\langle\frac{\partial\mathbb{F}}{\partial\bar{\xi}_{\vec{n}}}\bar{\xi}_{\vec{\ell}_{1}}\bar{\xi}_{\vec{\ell}_{2}}\bar{\xi}_{\vec{\ell}_{3}}\right\rangle\,.$
(25)

$\Gamma_{\left|\vec{n}\right|}$ can be thought of as an effective coupling
constant such that corrections to this linear theory are kept small. The value
of $\Gamma_{\left|\vec{n}\right|}$ will be determined in the continuation
though due to the isotropic character of our system, we have however already
assumed that it can only depend on the size of $\vec{n}$ and not its
direction. Now if $\left\langle\mathbb{F}\right\rangle^{\left(m\right)}$
denotes an $m^{th}$ order expansion of $\left\langle\mathbb{F}\right\rangle$,
then by assumption, the latter terms which follow from the nonlinear operator
$\left(\mathcal{O}-\mathcal{O}_{0}\right)$ will contribute at a higher order
and thus we can devise an iterative scheme

$0=\frac{1}{2}\sum_{\vec{n}}\left|\vec{n}\right|^{2}\left\langle\frac{\partial^{2}\mathbb{F}}{\partial\bar{\xi}_{\vec{n}}\partial\bar{\xi}_{-\vec{n}}}\right\rangle^{\left(m\right)}-\sum_{\vec{n}}\Gamma_{\left|\vec{n}\right|}\left\langle\frac{\partial\mathbb{F}}{\partial\bar{\xi}_{\vec{n}}}\bar{\xi}_{\vec{n}}\right\rangle^{\left(m\right)}\\\
-\sum_{\vec{n}}\left(g\left|\vec{n}\right|^{4}-\Gamma_{\left|\vec{n}\right|}\right)\left\langle\frac{\partial\mathbb{F}}{\partial\bar{\xi}_{\vec{n}}}\bar{\xi}_{\vec{n}}\right\rangle^{\left(m-1\right)}\\\
-\frac{1}{2}\sum_{\vec{n}}\sum_{\vec{\ell}_{1}\neq\vec{n}}\sum_{\vec{\ell}_{2}}\sum_{\vec{\ell}_{3}}V_{\vec{n},\vec{\ell}_{1},\vec{\ell}_{2},\vec{\ell}_{3}}\left\langle\frac{\partial\mathbb{F}}{\partial\bar{\xi}_{\vec{n}}}\bar{\xi}_{\vec{\ell}_{1}}\bar{\xi}_{\vec{\ell}_{2}}\bar{\xi}_{\vec{\ell}_{3}}\right\rangle^{\left(m-1\right)}$
(26)

relating higher order moments to lower order ones. This equation can now be
used to obtain any moment up to any order. For example, to obtain the second
moment at zeroth order, it is sufficient to sub
$\mathbb{F}=\bar{\xi}_{\vec{n}}\bar{\xi}_{\vec{n}^{\prime}}$ and $m=0$ into
equation (26), in which case the $(m-1)$ terms drop out and one immediately
obtains

$\left\langle\bar{\xi}_{\vec{n}}\bar{\xi}_{\vec{n}^{\prime}}\right\rangle^{\left(0\right)}=\frac{\left|\vec{n}\right|^{2}}{2\Gamma_{\left|\vec{n}\right|}}\delta_{\vec{n},-\vec{n}^{\prime}}\,.$
(27)

Similarly, subbing in
$\mathbb{F}=\bar{\xi}_{\vec{n}_{1}}\bar{\xi}_{\vec{n}_{2}}\bar{\xi}_{\vec{n}_{3}}\bar{\xi}_{\vec{n}_{4}}$
and $m=0$ ultimately gives the four-point function at zeroth order as

$\left\langle\bar{\xi}_{\vec{n}_{1}}\bar{\xi}_{\vec{n}_{2}}\bar{\xi}_{\vec{n}_{3}}\bar{\xi}_{\vec{n}_{4}}\right\rangle^{\left(0\right)}=\left\langle\bar{\xi}_{\vec{n}_{1}}\bar{\xi}_{\vec{n}_{2}}\right\rangle^{\left(0\right)}\left\langle\bar{\xi}_{\vec{n}_{3}}\bar{\xi}_{\vec{n}_{4}}\right\rangle^{\left(0\right)}\\\
+\left\langle\bar{\xi}_{\vec{n}_{1}}\bar{\xi}_{\vec{n}_{3}}\right\rangle^{\left(0\right)}\left\langle\bar{\xi}_{\vec{n}_{2}}\bar{\xi}_{\vec{n}_{4}}\right\rangle^{\left(0\right)}+\left\langle\bar{\xi}_{\vec{n}_{1}}\bar{\xi}_{\vec{n}_{4}}\right\rangle^{\left(0\right)}\left\langle\bar{\xi}_{\vec{n}_{2}}\bar{\xi}_{\vec{n}_{3}}\right\rangle^{\left(0\right)}\,,$
(28)

which is just Isserlis’ theorem [79, 80], also known as Wick’s theorem [78],
for the four-point function of a multivariate Gaussian distribution.

To study the effect of the nonlinearity, we proceed to higher orders. Subbing
in $\mathbb{F}=\bar{\xi}_{\vec{n}}\bar{\xi}_{\vec{n}^{\prime}}$ and $m=1$
gives the following expression for the two-point function at first order

$\left(\Gamma_{\left|\vec{n}\right|}+\Gamma_{\left|\vec{n}^{\prime}\right|}\right)\left\langle\bar{\xi}_{\vec{n}}\bar{\xi}_{\vec{n}^{\prime}}\right\rangle^{\left(1\right)}=\left|\vec{n}\right|^{2}\delta_{\vec{n},-\vec{n}^{\prime}}\\\
-2\left(g\left|\vec{n}\right|^{4}-\Gamma_{\left|\vec{n}\right|}\right)\left\langle\bar{\xi}_{\vec{n}}\bar{\xi}_{\vec{n}^{\prime}}\right\rangle^{\left(0\right)}\\\
-\frac{1}{2}\sum_{\begin{subarray}{c}\vec{\ell}_{1}\neq\vec{n}\\\
\vec{\ell}_{2},\vec{\ell}_{3}\end{subarray}}V_{\vec{n},\vec{\ell}_{1},\vec{\ell}_{2},\vec{\ell}_{3}}\left\langle\bar{\xi}_{\vec{n}^{\prime}}\bar{\xi}_{\vec{\ell}_{1}}\bar{\xi}_{\vec{\ell}_{2}}\bar{\xi}_{\vec{\ell}_{3}}\right\rangle^{\left(0\right)}\\\
-\frac{1}{2}\sum_{\begin{subarray}{c}\vec{\ell}_{1}\neq\vec{n}^{\prime}\\\
\vec{\ell}_{2},\vec{\ell}_{3}\end{subarray}}V_{\vec{n}^{\prime},\vec{\ell}_{1},\vec{\ell}_{2},\vec{\ell}_{3}}\left\langle\bar{\xi}_{\vec{n}}\bar{\xi}_{\vec{\ell}_{1}}\bar{\xi}_{\vec{\ell}_{2}}\bar{\xi}_{\vec{\ell}_{3}}\right\rangle^{\left(0\right)}$
(29)

where we have already found expressions for the zeroth order two-point and
four-point functions. Subbing in these expressions together with the kernel
$V_{\vec{n},\vec{\ell}_{1},\vec{\ell}_{2},\vec{\ell}_{3}}$ allows us to
simplify the sums to obtain

$\left\langle\bar{\xi}_{\vec{n}}\bar{\xi}_{\vec{n}^{\prime}}\right\rangle^{\left(1\right)}=\left\langle\bar{\xi}_{\vec{n}}\bar{\xi}_{\vec{n}^{\prime}}\right\rangle^{\left(0\right)}+\frac{\left\langle\bar{\xi}_{\vec{n}}\bar{\xi}_{\vec{n}^{\prime}}\right\rangle^{\left(0\right)}}{\Gamma_{\left|\vec{n}\right|}}\times\\\
\times\left(\Gamma_{\left|\vec{n}\right|}-g\left|\vec{n}\right|^{4}-\frac{1}{2}\sum_{\vec{\ell}\neq\vec{n}}\frac{\left|\vec{\ell}\,\right|^{2}\left|\vec{n}\times\vec{\ell}\,\right|^{4}}{\Gamma_{\left|\vec{\ell}\,\right|}\left|\vec{n}-\vec{\ell}\,\right|^{4}}\right)\,.$
(30)

This equation expresses the two-point function at first order as a correction
to the zeroth order two-point function. If we wish our expansion for the two-
point function to be meaningful, it must be the case that corrections to the
zeroth order approximation be small. Accordingly, self-consistency demands
that we select $\Gamma_{\left|\vec{n}\right|}$ so as to ensure that this be
the case however since we have complete freedom to choose
$\Gamma_{\left|\vec{n}\right|}$, we can go even further and demand that the
first order correction vanish entirely. There are of course other ways to
select $\Gamma_{\left|\vec{n}\right|}$ however previous work [71] has shown
this method to be extremely successful as it introduces closure into our
hierarchy of moment equations in a manner that ensures they are self-
consistently correct up to first order. Though we will not need to do so here,
this approach can be applied at any order to obtain a corresponding self-
consistent $m^{th}$ order theory. Each such theory will give a different
expression for $\Gamma_{\left|\vec{n}\right|}$ thus
$\Gamma_{\left|\vec{n}\right|}$ is not strictly an effective coupling constant
or renormalised bending rigidity but rather the coefficient of a particular
linear theory whose $m^{th}$ order expansion closely resembles the $m^{th}$
order expansion of the nonlinear theory we are interested in. For the first
order theory we are developing here, we obtain

$\Gamma_{\left|\vec{n}\right|}=g\left|\vec{n}\right|^{4}+\frac{1}{2}\sum_{\vec{\ell}\neq\vec{n}}\frac{\left|\vec{\ell}\,\right|^{2}\left|\vec{n}\times\vec{\ell}\,\right|^{4}}{\Gamma_{\left|\vec{\ell}\,\right|}\left|\vec{n}-\vec{\ell}\,\right|^{4}}\,,$
(31)

which is a two-dimensional nonlinear discrete integral equation for
$\Gamma_{\left|\vec{n}\right|}$. While solving such equations is in general
extremely challenging, in this particular case, the equation is indeed
amenable to analytic methods. In the appendix, we provide a detailed
derivation of the approximate solution

$\Gamma_{\left|\vec{n}\right|}\simeq
n^{4}\sqrt{\frac{3\pi}{4}\left[A\left(g\right)-\ln\left(\frac{n}{n_{max}}\right)\right]}\times\\\
\times\left[1+B\left(g\right)\frac{n}{n_{max}}\right]\,,$ (32)

where $n_{max}$ is some upper-frequency cut-off which must be imposed on the
equation. $A\left(g\right)$ and $B\left(g\right)$ are constants which only
depend on $g$ and are determined by the implicit solutions to equations (53)
and (54) in the appendix. Since the physically interesting case is $g\ll 1$,
it is more convenient to provide the small $g$ expansions

$\displaystyle A$ $\displaystyle\simeq$ $\displaystyle
0.137+0.336g+0.243g^{2}+0.112g^{3}+O\left(g^{4}\right),$ (33) $\displaystyle
B$ $\displaystyle\simeq$
$\displaystyle-0.265+0.360g-0.395g^{2}+0.311g^{3}+O\left(g^{4}\right).\qquad$
(34)

Combining with equation (27), we obtain that the structure factor up to first
order is

$\left\langle\bar{\xi}_{\vec{n}}\bar{\xi}_{\vec{n}^{\prime}}\right\rangle^{\left(1\right)}=\\\
\frac{\delta_{\vec{n},-\vec{n}^{\prime}}}{n^{2}\sqrt{3\pi\left[A\left(g\right)-\ln\left(\frac{n}{n_{max}}\right)\right]}\left[1+B\left(g\right)\frac{n}{n_{max}}\right]}\,.$
(35)

In Fig. 1, this structure factor
$S_{n}^{(1)}=\left\langle\bar{\xi}_{\vec{n}}\bar{\xi}_{-\vec{n}}\right\rangle^{\left(1\right)}$
is plotted for various values of $g$ together with a guide line for $\sim
n^{-2}$. Since the log-log plots appear straight over large intervals, they
give the illusion of a power law with anomalous exponent close to but not
exactly $-2$. In actual fact, it is the logarithmic and polynomial corrections
which causes the deviation from the non-anomalous power law. While for small
$n$, all values of $g$ tend towards the same asymptote, for large $n$, the
polynomial correction, determined by the coefficient $B\left(g\right)$,
becomes increasingly important, even giving rise to a small increasing tail
for extremely small values of $g$.

Figure 1: The first order approximation for the structure factor $S_{\vec{n}}$
of the sheet as given by equation (35) for various values of $g$ (solid and
dashed lines). The dot-dashed line is a guideline for the power law $\sim
n^{-2}$. For small $n$, the logarithmic correction can be seen to introduce a
small deviation from the power law, while for large $n$, the polynomial
correction $\left[1+B\left(g\right)\frac{n}{n_{max}}\right]$ introduces a
small increasing tail which becomes more pronounced with shrinking $g$.

Before moving to the next section, we would like to discuss the
renormalisation of the coupling constant $g$. Simple power counting of the
linear terms shows that the nonlinear quartic interaction scales as $L^{2-d}$,
where $d$ is the dimensionality of the sheet. Since $d=2$ for the problem of
interest, the nonlinear term is marginal and its relevance or irrelevance at
long distances can only be decided by a renormalization group calculation.
This point is also discussed in [55] where the assumption that the coupling
constant is not renormalised is assumed. It turns out that the self-consistent
expansion does not indicate any renormalisation of the coupling constant $g$
itself – at least not at the leading order used in this paper. It is however
important to stress that this option cannot be ruled out.

##  IV Comparison with Simulations

Equation (13), which describes the fluctuations of the sheet height in Fourier
space, is an ordinary Langevin field equation and thus in principle can be
simulated over say a square lattice by forward iteration with a small
discretisation of the time step. In practice, the conserved noise can easily
be computed in the Fourier space but computing the quartic interaction term
directly can be time consuming. To circumvent this problem, a pseudo-spectral
method, as described in [44], can be implemented. This entails expressing the
quartic interaction as the Fourier transform of its real space counterpart
which can be calculated far more quickly though this also has the effect of
imposing periodic boundary conditions on our fluctuating sheet. The fact that
a maximum frequency $n_{max}$ must be chosen to perform simulations poses no
problem as our solution, given by equation (35), indeed assumes the existence
of such a cut-off. Physically, this corresponds to the fact that no sheet can
be probed with infinite resolution. The time step $\delta t$ however must be
taken to be sufficiently small to avoid numerical instability. The maximum
size $\delta t$ can be taken while ensuring stability is determined by the
maximum frequency used, with larger values of $n_{max}$ requiring smaller time
steps. Accordingly $n_{max}$ cannot be too large if we wish to simulate
equation (13) over sufficiently long periods of time to obtain meaningful
results for the structure factor $S_{\vec{n}}$. The results shown in Fig. 2
were obtained for the frequencies
$\vec{n}\in\left[-20,20\right]\times\left[-20,20\right]$, ie. a $41\times 41$
square lattice such that the maximum frequency is $n_{max}=20\sqrt{2}\approx
28$. To obtain numerically stable results, a time step of $\delta t=10^{-6}$
was used and for each value of $g\in\left\\{1,0.1,0.01,0\right\\}$, $10$
simulations were run over $10^{6}$ time steps.

Figure 2: Numerical results for the structure factor
$S_{\vec{n}}=\left\langle\bar{\xi}_{\vec{n}}\bar{\xi}_{-\vec{n}}\right\rangle$
derived from simulating equation (13) with various values of $g$: (a) $g=1$,
(b) $g=0.1$, (c) $g=0.01$, (d) $g=0$. The simulations are not completely
isotropic thus for each simulation, the structure factor has been computed at
angles of $0^{\circ}$ (blue plus signs), $15^{\circ}$ (red circles),
$30^{\circ}$ (green triangles) and $45^{\circ}$ (purple squares) relative to
the “x-axis”. The solid black lines are fits of equation (35) to each data set
while the dashed lines are guidelines proportional to $n^{-2}$.

To validate the functional form of the structure factor $S_{\vec{n}}$
predicted by our theory, equation (35) with $A\left(g\right)$ and
$B\left(g\right)$ as free parameters was fitted (solid black lines) to each
data set shown in Fig. 2 (coloured shapes). It can be seen that no matter the
value of $g$, equation (35) indeed captures the shape of the structure factor.
The dashed line in each figure is a guideline proportional to $n^{-2}$ from
which it can easily be seen how the logarithmic correction in equation (35)
slightly modifies the slope of the structure factor.

The periodic boundary conditions over a square lattice have the effect of
breaking the isotropy of the system and thus the structure factor
$S_{\vec{n}}$ is found to also depend on the angle of $\vec{n}$ and not just
its magnitude. For this reason, in each figure, separate data has been plotted
for angles of $0^{\circ}$, $15^{\circ}$, $30^{\circ}$ and $45^{\circ}$.
Curiously, the degree of anisotropy decreases with increasing $g$ vanishing
completely as $g$ approaches 1. Even though our theory is derived under the
assumption of complete isotropy, the functional form of $S_{\vec{n}}$ remains
a superb fit for each data set in the presence of the anisotropy. Indeed, for
any given $g$, the anisotropy only seems to effect the tail parameter,
$B\left(g\right)$, while the parameter $A\left(g\right)$, related to the
logarithmic correction, is entirely unaffected by the angle of $\vec{n}$. This
robustness to slightly anisotropic circumstances can be viewed as further
validation of the functional form given by equation (35).

Figure 3: Comparisons between the theoretical values of $A(g)$ and $B(g)$
given by equations (33) and (34) (solid lines) and the values obtained from
fitting equation (35) to the simulation data (various shapes). The dashed
lines in the latter plot are angular dependent scalings $f(\theta)B(g)$ of
equation (34). Equation (33) and the scaled forms of equation (34) beautifully
predict the fitted value of $A(g)$ and $B(g)$ over a wide range of $g$ values.

Comparisons between the theoretical values of $A(g)$ and $B(g)$, given by
equations (33) and (34), and the values obtained from fitting equation (35) to
our simulation data are shown in Fig. 3. The theoretical values predicted by
equation (33) for $A(g)$ (solid line) can be seen to be extremely accurate
over a wide range of $g$ values, only breaking down as $g$ approaches 1. In
contrast, due to the anisotropic nature of $B(g)$, the theoretical values
predicted by equation (34) (solid line) do not capture any of the particular
fit values. One can however obtain superb fits by simply scaling the function
$B(g)$ given by equation (34) by a constant, that is, by finding for each
angle $\theta$ the quantity $f(\theta)$ which minimises the difference between
$f(\theta)B(g)$ and the numerical results for $B(g)$ 111We find
$f(0^{\circ})=-0.40$, $f(15^{\circ})=1.35$, $f(30^{\circ})=2.75$,
$f(45^{\circ})=3.40$. The dashed lines in Fig. 3 show these optimal scalings
and can be seen to be excellent fits for all angles over the entire range of
$g$ values. As mentioned above, the origin of this angular dependent scaling
is due to the periodic boundary conditions used in our simulation though it is
curious that this anisotropy only effects $B(g)$ and not $A(g)$. Though the
particular nature of $f(\theta)$ is not fully understood, these results inform
us that, at least in this case, a more precise expression for the structure
factor would be

$\left\langle\bar{\xi}_{\vec{n}}\bar{\xi}_{\vec{n}^{\prime}}\right\rangle^{\left(1\right)}=\\\
\frac{\delta_{\vec{n},-\vec{n}^{\prime}}}{n^{2}\sqrt{3\pi\left[A\left(g\right)-\ln\left(\frac{n}{n_{max}}\right)\right]}\left[1+B\left(g\right)\frac{n}{n_{max}/f(\theta)}\right]}$
(36)

where $A(g)$ and $B(g)$ are still given by equations (33) and (34)
respectively. This expression suggests the angular scaling is essentially just
a renormalisation of the maximum frequency supported by the sheet in each
direction as a result of the boundary conditions. As elaborated on in the
appendix, the function $B(g)$ is essentially determined by matching the
general solution of equation (31) with its large $\vec{n}$ boundary condition
and thus it is perhaps unsurprising that an anisotropic renomalisation of the
maximum supported frequency should imply a rescaling of
$B(g)\frac{n}{n_{max}}$.

##  V Discussion

We have combined the Föppl-von Kármán equations from elasticity theory with
out-of-equilibrium statistical physics to develop a mathematical model for
fluctuating thin sheets. Solving this model, we have shown that the structure
factor of such sheets will follow nontrivial scaling – a logarithmically
corrected power law for small frequencies which becomes a rational function
for large frequencies. The solution detailed in the appendix suggests that the
problem requires an upper frequency cutoff and also provides definite
relationships between the coefficients appearing in our solution and the
physical properties such as sheet thickness $h$ and Young’s modulus $E$. This
opens up the possibility for experimental study through the measurement of
such coefficients.

Indeed, while experimental investigation of fluctuating sheets has been
performed in the context of studying wave turbulence [38, 39, 40, 41, 42], we
are not aware of any attempt to measure the structure factor of the
fluctuating sheet itself. On the other hand, the structure factor of a
crumpled sheet has been measured several times [1, 3, 4, 5] and in each of
these instances, anomalous power-laws have been loosely fitted to the data.
Observed exponents cannot be simply explained by the structure factor of a
single ridge $\sim 1/q^{4}$ and as random crumpling is merely the
plastification of a particular set of random deformations of a thin sheet, it
is not unreasonable to hypothesise a connection between the structure of a
fluctuating sheet and that of a crumpled sheet. Accordingly, our results
suggest that the structure factor of crumpled sheets might be better described
by rational functions with logarithmic corrections, than by crude anomalous
power-laws. The wide range of applications of fluctuating sheets including
extremely thin sheets such as graphene suggests that many potential
experiments to test our results can be developed.

Furthermore, it is not only our results which have indicated the importance of
logarithmic corrections to the structure factor and renormalised bending
rigidity of thin sheets. In recent years, work carried out on thermally
fluctuating sheets in the presence of external tension [48, 50] have also
observed logarithmic effects, though these are not identical with ours. In
particular, our logarithmic correction occurs with a square root while theirs
occurs unadorned. This prevalence amongst related yet distinct models is
supportive of our argument that more sophisticated analyses which go beyond
simple power laws are needed for a proper understanding of such systems.

Perhaps the biggest limitation of our theory is the absence of self-avoidance.
As has been shown in [10], self-avoidance plays a significant role in the
context of strong crumpling. On the other hand, in the context of weak
crumpling, instances of self-avoidance will occur infrequently. Accordingly,
neglecting self-avoidance simply limits the scope of our theory rather than
its validity.

As described in the introduction, much research into crumpled paper has
focused on the formation of singular structures such as ridges and d-cones and
their interactions. Since these structures are known to solve the Föppl-von
Kármán equations [20, 21], it is interesting to ask whether the approach taken
in this paper can be connected directly to this research. For example, we
wonder whether our approach can be reframed as some sort of many body model
whose constituent parts are the various ridges and d-cones? If this were the
case, our approach could be conceived of as a type of “many crease”
generalisation of the theory of singular structures in thin sheets.

Further, we re-emphasise that the structure factor of fluctuating thin sheets
is an input into much contemporary research including the mechanical and
elastic properties of deformed sheets [46, 47, 48, 49, 50, 51] and their
acoustic [52, 53] and optical emissions [54]. While such research has so far
tended to assume ordinary power-laws, our results suggest that a wider family
of inputs needs to be considered including rational functions, possibly with
logarithmic corrections. Indeed, exploration of the effect of this wider range
of functions on acoustic emissions could suggest novel questions of inference
such as, to paraphrase Kac’s provocative question [82], can one hear the shape
of a fluctuating or crumpled sheet?

The approach we have taken in this paper naturally generalises to investigate
other types of thin sheet and forcing. For example, instead of beginning with
the flat sheet Föppl-von Kármán equations, we could begin with the Föppl-von
Kármán equations for shells [83, 49] or other types of curved surfaces [84]
and this would allow us to investigate the effect of the underlying undeformed
geometry on the fluctuating structure. Similarly, it would be of great
interest to incorporate more realistic modelling of the forcing such as
driving from the boundaries which would amount to complex spatio-temporal
noise. Accordingly, the approach we have introduced here is extremely generic,
providing a first step to analytically describing the behaviour of a wide
range of different types of fluctuating surface.

###### Acknowledgements.

The authors wish to thank one of the anonymous referees for insightful
comments and analysis that helped to clarify the origin of the $\sqrt{\ln{q}}$
correction to the structure factor. This work was supported by the Israel
Science Foundation Grant No. 1682/18.

*

##  Appendix A Solution to the Integral Equation

In this appendix, we outline a method for solving equation (31)

$\Gamma_{\left|\vec{n}\right|}=g\left|\vec{n}\right|^{4}+\frac{1}{2}\sum_{\vec{\ell}\neq\vec{n}}\frac{\left|\vec{\ell}\,\right|^{2}\left|\vec{n}\times\vec{\ell}\,\right|^{4}}{\Gamma_{\left|\vec{\ell}\,\right|}\left|\vec{n}-\vec{\ell}\,\right|^{4}}\,,$
(37)

a two-dimensional nonlinear discrete integral equation for the parameter
$\Gamma_{\left|\vec{n}\right|}$. To obtain an approximate solution, we can
consider its continuous limit. In particular, let
$\vec{q}=\vec{n}/L,\vec{k}=\vec{\ell}/L$ and
$\Gamma\left(q\right)=\Gamma_{\left|\vec{n}\right|}/\left(\Lambda
L\right)^{4}$ such that

$\Gamma\left(q\right)=g\frac{\left|\vec{q}\,\right|^{4}}{\Lambda^{4}}+\frac{1}{2}\frac{1}{\Lambda^{8}}\sum_{\vec{k}\neq\vec{q}}\frac{1}{L^{2}}\frac{\left|\vec{k}\,\right|^{2}\left|\vec{q}\times\vec{k}\,\right|^{4}}{\Gamma\left(k\right)\left|\vec{q}-\vec{k}\,\right|^{4}}\,.$
(38)

Here, $\Lambda$ denotes some frequency scale which ensures that
$\Gamma\left(q\right)$ remains dimensionless. In the large sheet limit of
$L\rightarrow\infty$, the sum becomes an integral and we obtain a two-
dimensional continuous integral equation for $\Gamma\left(q\right)$

$\Gamma\left(q\right)=g\frac{\left|\vec{q}\,\right|^{4}}{\Lambda^{4}}+\frac{1}{2}\frac{1}{\Lambda^{8}}\int
d^{2}k\frac{\left|\vec{k}\,\right|^{2}\left|\vec{q}\times\vec{k}\,\right|^{4}}{\Gamma\left(k\right)\left|\vec{q}-\vec{k}\,\right|^{4}}\,.$
(39)

Since $\Gamma\left(k\right)$ only depends on the size of $\vec{k}$ and not its
direction, the angular part of the integral can be carried out exactly giving

$\Gamma\left(q\right)=g\frac{q^{4}}{\Lambda^{4}}+\frac{3\pi}{8}\frac{1}{\Lambda^{8}}\left(\int_{0}^{q}dk\,\frac{k^{7}}{\Gamma\left(k\right)}+q^{4}\int_{q}^{\Lambda}dk\,\frac{k^{3}}{\Gamma\left(k\right)}\right)\,.$
(40)

Our task has thus been reduced to finding the solution to this one-dimensional
continuous integral equation. In the continuation, it will become apparent
that this integral equation only has a solution over finite intervals
$q\in\left[0,q_{max}\right]$ where $q_{max}<\infty$ is some finite upper
frequency cut-off. As such, we have chosen to identify the previously
arbitrary frequency scale $\Lambda$ with the upper cut-off and introduced it
into the integral bounds.

To solve this equation, differentiate it twice with respect to $q$ to obtain

$\frac{d\Gamma}{dq}=4g\frac{q^{3}}{\Lambda^{4}}+\frac{3\pi}{2}\frac{q^{3}}{\Lambda^{8}}\int_{q}^{\Lambda}dk\,\frac{k^{3}}{\Gamma\left(k\right)}$
(41)

and

$\frac{d^{2}\Gamma}{dq^{2}}=12g\frac{q^{2}}{\Lambda^{4}}+\frac{3\pi}{2}\frac{1}{\Lambda^{8}}\left(3q^{2}\int_{q}^{\Lambda}dk\,\frac{k^{3}}{\Gamma\left(k\right)}-\frac{q^{6}}{\Gamma\left(q\right)}\right)\,.$
(42)

These equations can now be combined to eliminate the integral giving rise to a
second order nonlinear ODE

$\frac{d^{2}\Gamma}{dq^{2}}=\frac{3}{q}\frac{d\Gamma}{dq}-\frac{3\pi}{2}\frac{1}{\Lambda^{8}}\frac{q^{6}}{\Gamma\left(q\right)}\,.$
(43)

Apart from being far easier to solve than the integral equation, this equation
has the added benefit that the previously explicit $g$ dependence has dropped
out. Instead, the $g$ dependence now enters in the form of boundary
conditions. In particular, the limits $q\rightarrow 0$ in equation (40) and
$q\rightarrow\Lambda$ in equation (41) give rise to the boundary conditions

$\displaystyle\Gamma\left(q=0\right)$ $\displaystyle=0\,,$ (44)
$\displaystyle\left.\frac{d\Gamma}{dq}\right|_{q=\Lambda}$
$\displaystyle=\dfrac{4g}{\Lambda}\,.$ (45)

The general solutions to equation (43) are easily obtained numerically with
arbitrary starting values. These solutions all exhibit movable singularities
[85] such that for sufficiently large $q$, $\Gamma\left(q\right)$ plummets to
0 and beyond this movable singularity, the solution is no longer defined. On
this basis, the aforementioned upper frequency cut-off $\Lambda$ was
identified.

To investigate the solution to equation (43), we begin by studying the small
$q$ behaviour. Expanding $\Gamma\left(q\right)$ at lowest order in a power
series

$\Gamma\left(q\right)=Aq^{\nu}+O\left(q^{\nu+1}\right)$ (46)

results in the equation

$\left(4-\nu\right)\nu
A^{2}q^{2\nu}+O\left(q^{2\nu+1}\right)=\frac{3\pi}{2}\frac{1}{\Lambda^{8}}q^{8}\,.$
(47)

Comparing the exponents of $q$ on the left and right-hand sides of this
equation, we find that only $\nu=4$ can solve the equation at lowest order,
however the factor of $\left(4-\nu\right)$ on the left-hand side ensures that
$\nu=4$ will also fail to actually balance the left and right-hand sides. This
observation forces the conclusion that the small $q$ expansion of
$\Gamma\left(q\right)$ is not a power series. To obtain the actual small $q$
behaviour, the following heuristic argument can be made. Setting
$\Gamma(q)=q^{4}u(q)$, equation (43) can be written as an ODE for $u(q)$

$q^{2}u\frac{d^{2}u}{dq^{2}}+\frac{5}{2}q\frac{d\left(u^{2}\right)}{dq}=-\frac{3\pi}{2}\frac{1}{\Lambda^{8}}\,.$
(48)

Neglecting the second derivative term, this equation is exactly solvable,
having solution

$u\left(q\right)=\sqrt{A-\frac{3\pi}{5}\frac{1}{\Lambda^{8}}\ln\left(q\right)}$
(49)

where $A$ is some constant. Accordingly, we might expect the small $q$
behaviour of $\Gamma(q)$ to be logarithmically corrected by a term
proportional to $\sqrt{\ln(q)}$ and indeed, without neglecting any terms, one
can find upon substitution that for small $q$, equation (43) is solved by

$\Gamma\left(q\right)=\left(\frac{q}{\Lambda}\right)^{4}\sqrt{\frac{3\pi}{4}\left[A-\ln\left(\frac{q}{\Lambda}\right)\right]}\times\\\
\times\left[1+O\left(\frac{\ln\left(\ln\left(\Lambda/q\right)\right)}{\ln\left(\Lambda/q\right)}\right)\right]$
(50)

where $A$ is a constant determined by the boundary condition at $q=\Lambda$.
For example, requiring that this equation match with equation (45) gives

$A=\frac{1}{8}+\frac{2}{3\pi}g^{2}+\frac{g}{\sqrt{6\pi}}\sqrt{1+\frac{8}{3\pi}g^{2}}\,.$
(51)

Since equation (50) is established around $q=0$ while equation (45) imposes a
boundary condition far from $q=0$, there is little reason to trust this
approximation as we move away from $q=0$. To obtain reliable results across
the entire domain $q\in\left[0,\Lambda\right]$, we can also expand the
solution to equation (43) around $q=\Lambda$ and stitch the expansions
together. To this end, we propose a type of two-point expansion that uses the
known analytic behaviour around the two end points [86]

$\Gamma\left(q\right)\simeq\left(\frac{q}{\Lambda}\right)^{4}\sqrt{\frac{3\pi}{4}\left[A-\ln\left(\frac{q}{\Lambda}\right)\right]}\left[1+B\frac{q}{\Lambda}\right]$
(52)

where $B$ is another constant. For $q$ close to 0, this modification has no
effect while for $q$ close to $\Lambda$, $B$ can be chosen to ensure that
$\Gamma\left(q\right)$ be correct up to zeroth order. In particular, subbing
this expression into the boundary condition given by equation (45) gives

$\frac{16}{\sqrt{3\pi}}g\sqrt{A}=\left(8A-1\right)+\left(10A-1\right)B\,,$
(53)

while subbing it into equation (43) and expanding around $q=\Lambda$ implies
that at lowest order, we must have

$8A=\left(1+8A+B+12AB-20A^{2}B\right)\left(1+B\right)\,.$ (54)

This pair of equations determines the constants $A$ and $B$. It is worth
mentioning that the equations have multiple real solutions but only one of
these solutions actually generates an approximation for $\Gamma\left(q\right)$
which approximately satisfies the original integral equation given by equation
(40). This useful solution is plotted in Fig. 4 as a function of $g$.
Unfortunately, this single solution ceases to exist for values of $g$ greater
than roughly 0.912 however since we are interested in small $g$, this poses
little difficulty to us in practice. Furthermore, as $g$ grows, it turns out
that using equation (50) with equation (51) becomes an increasingly good
approximation. Considering higher order corrections is also an option though
this results in little gain at the cost of significantly increased complexity.

Figure 4: $A$ and $B$ as functions of $g$. The solid lines are the exact
solution obtained from equations (53) and (54). The dashed lines are the power
series expansions around $g=0$ given by equations (56) and (57).

It is not difficult to numerically evaluate equations (53) and (54) for any
particular value of $g$. Since we are primarily interested in the case where
$g$ is small, however, it is more convenient to have a small $g$ expansion for
$A$ and $B$. This also has the added advantage of avoiding the spurious
solutions of equations (53) and (54). To this end, we expand

$A(g)=\sum_{n=0}^{\infty}A_{n}g^{n}\qquad B(g)=\sum_{n=0}^{\infty}B_{n}g^{n}$
(55)

where the $A_{n}$ and $B_{n}$ are constants. Subbing these expressions into
equations (53) and (54) and equating coefficients of $g$ order by order, we
obtain up to $O\left(g^{4}\right)$

$A(g)\simeq 0.137+0.336g+0.243g^{2}+\\\ +0.112g^{3}+O\left(g^{4}\right)\,,$
(56) $B(g)\simeq-0.265+0.360g-0.395g^{2}+\\\
+0.311g^{3}+O\left(g^{4}\right)\,.$ (57)

These approximations are plotted as the dashed lines in Fig. 4 and can be seen
to be faithful for values of $g\lesssim 0.4$. Fig. 5 shows a comparison
between equation (52) with $A(g)$ and $B(g)$ given by equations (56) and (57)
(dots) and direct numerical solutions of equation (43) (solid and dashed
lines). The numerical solutions were obtained with $\Lambda=1$ via a shooting
method [87] in which numerical integration of equation (43) starting at
$q=\Lambda$ was performed with varying initial values of $\Gamma(q=\Lambda)$
until solutions were found with the desired boundary condition given by
equation (44). As can be seen, for small $g$, equation (52) with $A(g)$ and
$B(g)$ given by equations (56) and (57) is a superb approximation for the
numerical solution though the quality decreases once $g$ reaches the size of
unity. It is understood that this is because the series expansions for $A(g)$
and $B(g)$ are only valid for small $g$ and indeed, as stated above, for non-
small values of $g$ a better approximation is simply equation (50) with $A(g)$
given by equation (51). Finally, we note that solutions to equation (43) with
differing frequency cut-offs $\Lambda$ are self-similar, that is, if
$\Gamma(q)$ is a solution for a given maximum frequency $\Lambda$, then the
scaled solution $\Gamma(aq)/a^{4}$ with $a>0$ will be a solution with maximum
frequency $a\Lambda$. Accordingly, the validity of the comparison shown in
Fig. 5 is not limited to the case where $\Lambda=1$ since all other cases can
be directly obtained by merely scaling the results.

Figure 5: Comparison between numerical solutions of equation (43) (solid and
dashed lines) and the theory predicted by equation (52) with $A(g)$ and $B(g)$
given by equations (56) and (57) (dots). For small $g$, the theory exhibits
excellent agreement with the numerical solutions.

Returning to the discrete system we began with, we obtain an approximate
solution for $\Gamma_{\left|\vec{n}\right|}$ is

$\Gamma_{\left|\vec{n}\right|}\simeq\left(\Lambda
L\right)^{4}\Gamma\left(q=\frac{n}{L}\right)\,,$ (58)

or explicitly

$\Gamma_{\left|\vec{n}\right|}\simeq
n^{4}\sqrt{\frac{3\pi}{4}\left[A-\ln\left(\frac{n}{n_{max}}\right)\right]}\times\\\
\times\left[1+B\frac{n}{n_{max}}\right]\,,$ (59)

where $n_{max}=\Lambda L$ denotes a discrete upper-frequency cut-off.

## References

  * Plouraboué and Roux [1996] F. Plouraboué and S. Roux, Experimental study of the roughness of crumpled surfaces, Physica A 227, 173 (1996).
  * Matan  _et al._ [2002] K. Matan, R. B. Williams, T. A. Witten, and S. R. Nagel, Crumpling a thin sheet, Phys. Rev. Lett. 88, 076101 (2002).
  * Blair and Kudrolli [2005] D. L. Blair and A. Kudrolli, Geometry of crumpled paper, Phys. Rev. Lett. 94, 166107 (2005).
  * Balankin  _et al._ [2006] A. S. Balankin, O. S. Huerta, R. Cortes Montes de Oca, D. S. Ochoa, J. Martínez Trinidad, and M. A. Mendoza, Intrinsically anomalous roughness of randomly crumpled thin sheets, Phys. Rev. E 74, 061602 (2006).
  * Andresen  _et al._ [2007] C. A. Andresen, A. Hansen, and J. Schmittbuhl, Ridge network in crumpled paper, Phys. Rev. E 76, 026108 (2007).
  * Balankin and Huerta [2008] A. S. Balankin and O. S. Huerta, Entropic rigidity of a crumpling network in a randomly folded thin sheet, Phys. Rev. E 77, 051124 (2008).
  * Deboeuf  _et al._ [2013] S. Deboeuf, E. Katzav, A. Boudaoud, D. Bonn, and M. Adda-Bedia, Comparative study of crumpling and folding of thin sheets, Phys. Rev. Lett. 110, 104301 (2013).
  * Balankin  _et al._ [2013] A. S. Balankin, A. Horta Rangel, G. García Pérez, F. Gayosso Martinez, H. Sanchez Chavez, and C. L. Martínez-González, Fractal features of a crumpling network in randomly folded thin matter and mechanics of sheet crushing, Phys. Rev. E 87, 052806 (2013).
  * Gottesman  _et al._ [2018] O. Gottesman, J. Andrejevic, C. H. Rycroft, and S. M. Rubinstein, A state variable for crumpled thin sheets, Communications Physics 1, 1 (2018).
  * Vliegenthart and Gompper [2006] G. Vliegenthart and G. Gompper, Forced crumpling of self-avoiding elastic sheets, Nature materials 5, 216 (2006).
  * Sultan and Boudaoud [2006] E. Sultan and A. Boudaoud, Statistics of crumpled paper, Phys. Rev. Lett. 96, 136103 (2006).
  * Andrejevic  _et al._ [2021] J. Andrejevic, L. M. Lee, S. M. Rubinstein, and C. H. Rycroft, A model for the fragmentation kinetics of crumpled thin sheets, Nature communications 12, 1470 (2021).
  * Kantor  _et al._ [1986] Y. Kantor, M. Kardar, and D. R. Nelson, Statistical mechanics of tethered surfaces, Phys. Rev. Lett. 57, 791 (1986).
  * Kantor  _et al._ [1987] Y. Kantor, M. Kardar, and D. R. Nelson, Tethered surfaces: Statics and dynamics, Phys. Rev. A 35, 3056 (1987).
  * Kantor and Nelson [1987] Y. Kantor and D. R. Nelson, Crumpling transition in polymerized membranes, Phys. Rev. Lett. 58, 2774 (1987).
  * Kardar and Nelson [1987] M. Kardar and D. R. Nelson, $\epsilon$ expansions for crumpled manifolds, Phys. Rev. Lett. 58, 1289 (1987).
  * Paczuski  _et al._ [1988] M. Paczuski, M. Kardar, and D. R. Nelson, Landau theory of the crumpling transition, Phys. Rev. Lett. 60, 2638 (1988).
  * Duplantier  _et al._ [1990] B. Duplantier, T. Hwa, and M. Kardar, Self-avoiding crumpled manifolds: Perturbative analysis and renormalizability, Phys. Rev. Lett. 64, 2022 (1990).
  * Lobkovsky  _et al._ [1995] A. Lobkovsky, S. Gentges, H. Li, D. Morse, and T. A. Witten, Scaling properties of stretching ridges in a crumpled elastic sheet, Science 270, 1482 (1995).
  * Lobkovsky and Witten [1997] A. E. Lobkovsky and T. A. Witten, Properties of ridges in elastic membranes, Phys. Rev. E 55, 1577 (1997).
  * Ben Amar and Pomeau [1997] M. Ben Amar and Y. Pomeau, Crumpled paper, Proceedings of the Royal Society of London. Series A: Mathematical, Physical and Engineering Sciences 453, 729 (1997).
  * Chaïeb and Melo [1997] S. Chaïeb and F. Melo, Experimental study of crease formation in an axially compressed sheet, Phys. Rev. E 56, 4736 (1997).
  * Cerda  _et al._ [1999] E. Cerda, S. Chaieb, F. Melo, and L. Mahadevan, Conical dislocations in crumpling, Nature 401, 46 (1999).
  * Mora and Boudaoud [2002] T. Mora and A. Boudaoud, Thin elastic plates: On the core of developable cones, EPL (Europhysics Letters) 59, 41 (2002).
  * Liang and Witten [2005] T. Liang and T. A. Witten, Crescent singularities in crumpled sheets, Phys. Rev. E 71, 016612 (2005).
  * Landau  _et al._ [1986] L. D. Landau, E. M. Lifshitz, A. M. Kosevich, and L. P. Pitaevskii,  _Theory of elasticity: volume 7_ , 3rd ed., Vol. 7 (Elsevier, 1986).
  * Gladden  _et al._ [2005] J. R. Gladden, N. Z. Handzy, A. Belmonte, and E. Villermaux, Dynamic buckling and fragmentation in brittle rods, Phys. Rev. Lett. 94, 035503 (2005).
  * Audoly and Neukirch [2005] B. Audoly and S. Neukirch, Fragmentation of rods by cascading cracks: Why spaghetti does not break in half, Phys. Rev. Lett. 95, 095505 (2005).
  * Nelson  _et al._ [2004] D. Nelson, T. Piran, and S. Weinberg,  _Statistical mechanics of membranes and surfaces_ , 2nd ed. (World Scientific, 2004).
  * Liang and Purohit [2016] X. Liang and P. K. Purohit, A fluctuating elastic plate and a cell model for lipid membranes, Journal of the Mechanics and Physics of Solids 90, 29 (2016).
  * Meyer  _et al._ [2007a] J. C. Meyer, A. K. Geim, M. I. Katsnelson, K. S. Novoselov, T. J. Booth, and S. Roth, The structure of suspended graphene sheets, Nature 446, 60 (2007a).
  * Meyer  _et al._ [2007b] J. C. Meyer, A. Geim, M. Katsnelson, K. Novoselov, D. Obergfell, S. Roth, C. Girit, and A. Zettl, On the roughness of single-and bi-layer graphene membranes, Solid State Communications 143, 101 (2007b).
  * Thompson-Flagg  _et al._ [2009] R. C. Thompson-Flagg, M. J. Moura, and M. Marder, Rippling of graphene, EPL (Europhysics Letters) 85, 46002 (2009).
  * Zang  _et al._ [2013] J. Zang, S. Ryu, N. Pugno, Q. Wang, Q. Tu, M. J. Buehler, and X. Zhao, Multifunctionality and control of the crumpling and unfolding of large-area graphene, Nature materials 12, 321 (2013).
  * Deng and Berry [2016] S. Deng and V. Berry, Wrinkled, rippled and crumpled graphene: an overview of formation mechanism, electronic properties, and applications, Materials Today 19, 197 (2016).
  * Ahmadpoor  _et al._ [2017] F. Ahmadpoor, P. Wang, R. Huang, and P. Sharma, Thermal fluctuations and effective bending stiffness of elastic thin sheets and graphene: A nonlinear analysis, Journal of the Mechanics and Physics of Solids 107, 294 (2017).
  * Düring  _et al._ [2006] G. Düring, C. Josserand, and S. Rica, Weak turbulence for a vibrating plate: Can one hear a kolmogorov spectrum?, Phys. Rev. Lett. 97, 025503 (2006).
  * Boudaoud  _et al._ [2008] A. Boudaoud, O. Cadot, B. Odille, and C. Touzé, Observation of wave turbulence in vibrating plates, Phys. Rev. Lett. 100, 234504 (2008).
  * Mordant [2008] N. Mordant, Are there waves in elastic wave turbulence?, Phys. Rev. Lett. 100, 234505 (2008).
  * Cadot  _et al._ [2008] O. Cadot, A. Boudaoud, and C. Touzé, Statistics of power injection in a plate set into chaotic vibration, The European Physical Journal B 66, 399 (2008).
  * Cobelli  _et al._ [2009] P. Cobelli, P. Petitjeans, A. Maurel, V. Pagneux, and N. Mordant, Space-time resolved wave turbulence in a vibrating plate, Phys. Rev. Lett. 103, 204301 (2009).
  * Humbert  _et al._ [2013] T. Humbert, O. Cadot, G. Düring, C. Josserand, S. Rica, and C. Touzé, Wave turbulence in vibrating plates: the effect of damping, EPL (Europhysics Letters) 102, 30002 (2013).
  * Düring  _et al._ [2015] G. Düring, C. Josserand, and S. Rica, Self-similar formation of an inverse cascade in vibrating elastic plates, Phys. Rev. E 91, 052916 (2015).
  * Düring  _et al._ [2017] G. Düring, C. Josserand, and S. Rica, Wave turbulence theory of elastic plates, Physica D: Nonlinear Phenomena 347, 42 (2017).
  * Düring  _et al._ [2019] G. Düring, C. Josserand, G. Krstulovic, and S. Rica, Strong turbulence for vibrating plates: Emergence of a kolmogorov spectrum, Physical Review Fluids 4, 064804 (2019).
  * Košmrlj and Nelson [2013] A. Košmrlj and D. R. Nelson, Mechanical properties of warped membranes, Phys. Rev. E 88, 012136 (2013).
  * Košmrlj and Nelson [2014] A. Košmrlj and D. R. Nelson, Thermal excitations of warped membranes, Phys. Rev. E 89, 022126 (2014).
  * Košmrlj and Nelson [2016] A. Košmrlj and D. R. Nelson, Response of thermalized ribbons to pulling and bending, Phys. Rev. B 93, 125431 (2016).
  * Košmrlj and Nelson [2017] A. Košmrlj and D. R. Nelson, Statistical mechanics of thin spherical shells, Phys. Rev. X 7, 011002 (2017).
  * Shankar and Nelson [2021] S. Shankar and D. R. Nelson, Thermalized buckling of isotropically compressed thin sheets, Phys. Rev. E 104, 054141 (2021).
  * Morshedifard  _et al._ [2021] A. Morshedifard, M. Ruiz-García, M. J. A. Qomi, and A. Košmrlj, Buckling of thermalized elastic sheets, Journal of the Mechanics and Physics of Solids 149, 104296 (2021).
  * Kramer and Lobkovsky [1996] E. M. Kramer and A. E. Lobkovsky, Universal power law in the noise from a crumpled elastic sheet, Phys. Rev. E 53, 1465 (1996).
  * Houle and Sethna [1996] P. A. Houle and J. P. Sethna, Acoustic emission from crumpling paper, Phys. Rev. E 54, 278 (1996).
  * Rad  _et al._ [2019] V. F. Rad, E. E. Ramírez-Miquet, H. Cabrera, M. Habibi, and A.-R. Moradi, Speckle pattern analysis of crumpled papers, Applied optics 58, 6549 (2019).
  * Nelson and Peliti [1987] D. Nelson and L. Peliti, Fluctuations in membranes with crystalline and hexatic order, J. Phys. France 48, 1085 (1987).
  * Kleinert and Schulte-Frohlinde [2001] H. Kleinert and V. Schulte-Frohlinde,  _Critical properties of phi4-theories_ (World Scientific, 2001).
  * Hohenberg and Halperin [1977] P. C. Hohenberg and B. I. Halperin, Theory of dynamic critical phenomena, Rev. Mod. Phys. 49, 435 (1977).
  * McComb [2003] W. D. McComb,  _Renormalization methods: a guide for beginners_ (OUP Oxford, 2003).
  * Schwartz and Edwards [1992] M. Schwartz and S. Edwards, Nonlinear deposition: a new approach, EPL (Europhysics Letters) 20, 301 (1992).
  * Schwartz and Edwards [1998] M. Schwartz and S. F. Edwards, Peierls-boltzmann equation for ballistic deposition, Phys. Rev. E 57, 5730 (1998).
  * Katzav and Schwartz [1999] E. Katzav and M. Schwartz, Self-consistent expansion for the kardar-parisi-zhang equation with correlated noise, Phys. Rev. E 60, 5677 (1999).
  * Katzav and Schwartz [2002] E. Katzav and M. Schwartz, Existence of the upper critical dimension of the kardar–parisi–zhang equation, Physica A: Statistical Mechanics and its Applications 309, 69 (2002).
  * Schwartz and Edwards [2002] M. Schwartz and S. Edwards, Stretched exponential in non-linear stochastic field theories, Physica A: Statistical Mechanics and its Applications 312, 363 (2002).
  * Katzav [2002] E. Katzav, Self-consistent expansion for the molecular beam epitaxy equation, Phys. Rev. E 65, 032103 (2002).
  * Katzav [2003] E. Katzav, Self-consistent expansion results for the nonlocal kardar-parisi-zhang equation, Phys. Rev. E 68, 046113 (2003).
  * Katzav and Schwartz [2004a] E. Katzav and M. Schwartz, Kardar-parisi-zhang equation with temporally correlated noise: A self-consistent approach, Phys. Rev. E 70, 011601 (2004a).
  * Katzav and Schwartz [2004b] E. Katzav and M. Schwartz, Numerical evidence for stretched exponential relaxations in the kardar-parisi-zhang equation, Phys. Rev. E 69, 052603 (2004b).
  * Katzav and Adda-Bedia [2006] E. Katzav and M. Adda-Bedia, Roughness of tensile crack fronts in heterogenous materials, EPL (Europhysics Letters) 76, 450 (2006).
  * Katzav  _et al._ [2007] E. Katzav, M. Adda-Bedia, M. Ben Amar, and A. Boudaoud, Roughness of moving elastic lines: Crack and wetting fronts, Phys. Rev. E 76, 051601 (2007).
  * Edwards and Schwartz [2002] S. F. Edwards and M. Schwartz, Lagrangian statistical mechanics applied to non-linear stochastic field equations, Physica A: Statistical Mechanics and its Applications 303, 357 (2002).
  * Schwartz and Katzav [2008] M. Schwartz and E. Katzav, The ideas behind self-consistent expansion, Journal of Statistical Mechanics: Theory and Experiment 2008, P04023 (2008).
  * Remez and Goldstein [2018] B. Remez and M. Goldstein, From divergent perturbation theory to an exponentially convergent self-consistent expansion, Phys. Rev. D 98, 056017 (2018).
  * Le Doussal and Radzihovsky [1992] P. Le Doussal and L. Radzihovsky, Self-consistent theory of polymerized membranes, Physical review letters 69, 1209 (1992).
  * Le Doussal and Radzihovsky [2018] P. Le Doussal and L. Radzihovsky, Anomalous elasticity, fluctuations and disorder in elastic membranes, Annals of Physics 392, 340 (2018).
  * Risken [1996] H. Risken,  _The Fokker-Planck Equation_ , 2nd ed. (Springer, 1996) pp. 63–95.
  * Balescu [1975] R. Balescu,  _Equilibrium and Non-Equilibrium Statistical Mechanics_, A Wiley interscience publication (Wiley, 1975).
  * Plischke and Bergersen [2006] M. Plischke and B. Bergersen,  _Equilibrium statistical physics_ , 3rd ed. (World Scientific, 2006).
  * Kardar [2007] M. Kardar,  _Statistical physics of particles_ (Cambridge University Press, 2007).
  * Isserlis [1916] L. Isserlis, On certain probable errors and correlation coefficients of multiple frequency distributions with skew regression, Biometrika 11, 185 (1916).
  * Isserlis [1918] L. Isserlis, On a formula for the product-moment coefficient of any order of a normal frequency distribution in any number of variables, Biometrika 12, 134 (1918).
  * Note [1] We find $f(0^{\circ})=-0.40$, $f(15^{\circ})=1.35$, $f(30^{\circ})=2.75$, $f(45^{\circ})=3.40$.
  * Kac [1966] M. Kac, Can one hear the shape of a drum?, The american mathematical monthly 73, 1 (1966).
  * Paulose  _et al._ [2012] J. Paulose, G. A. Vliegenthart, G. Gompper, and D. R. Nelson, Fluctuating shells under pressure, Proceedings of the National Academy of Sciences 109, 19551 (2012).
  * van Hemmen and Peterson [2005] J. L. van Hemmen and M. A. Peterson, Crumpling of curved sheets: Generalizing foeppl-von karman, arXiv preprint cond-mat/0503158 (2005).
  * Bender and Orszag [1999] C. M. Bender and S. A. Orszag,  _Advanced mathematical methods for scientists and engineers I_ (Springer-Verlag, 1999).
  * Baker and Graves-Morris [1996] G. A. Baker and P. Graves-Morris,  _Pade Approximants: Encyclopedia of Mathematics and It’s Applications, Vol. 59_ , 2nd ed. (Cambridge University Press, 1996).
  * Press  _et al._ [2007] W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery,  _Numerical recipes 3rd edition: The art of scientific computing_ (Cambridge university press, 2007).