# Thermally driven elastic membranes are quasi-linear across all scales

Chanania Steinbock, Eytan Katzav  Racah Institute of Physics, The Hebrew
University, Jerusalem 9190401, Israel<EMAIL_ADDRESS>
###### Abstract

We study the static and dynamic structure of thermally fluctuating elastic
thin sheets by investigating a model known as the overdamped dynamic Föppl-von
Kármán equation, in which the Föppl-von Kármán equation from elasticity theory
is driven by white noise. The resulting nonlinear equation is governed by a
single nondimensional coupling parameter $g$ where large and small values of
$g$ correspond to weak and strong nonlinear coupling respectively. By
analysing the weak coupling case with ordinary perturbation theory and the
strong coupling case with a self-consistent methodology known as the self-
consistent expansion, precise analytic predictions for the static and dynamic
structure factors are obtained. Importantly, the maximum frequency $n_{\max}$
supported by the system plays a role in determining which of three possible
classes such sheets belong to: (1) when $g\gg 1$, the system is mostly linear
with roughness exponent $\zeta=1$ and dynamic exponent $z=4$, (2) when $g\ll
2/n_{\max}$, the system is extremely nonlinear with roughness exponent
$\zeta=1/2$ and dynamic exponent $z=3$, (3) between these regimes, an
intermediate behaviour is obtained in which a crossover occurs such that the
nonlinear behaviour is observed for small frequencies while the linear
behaviour is observed for large frequencies, and thus the large frequency
linear tail is found to have a significant impact on the small frequency
behaviour of the sheet. Back-of-the-envelope calculations suggest that ultra-
thin materials such as graphene lie in this intermediate regime. Despite the
existence of these three distinct behaviours, the decay rate of the dynamic
structure factor is related to the static structure factor as if the system
were completely linear. This quasi-linearity occurs regardless of the size of
$g$ and at all length scales. Numerical simulations confirm the existence of
the three classes of behaviour and the quasi-linearity of all classes.

††:

Keywords: Thin sheets, out-of-equilibrium dynamics, Föppl-von Kármán
equations, Family-Vicsek scaling, self-consistent expansion, roughness

##  1 Introduction

Thin matter appears at all scales in nature, ranging from everyday sheets of
paper to microscopic cell membranes to geological tectonic plates.
Accordingly, the theory of how thin sheets behave in the presence of thermal
or uncorrelated noise has extremely wide applicability and indeed, the
relatively recent synthesis of ultra-thin materials such as graphene [1] has
made the theory more relevant than ever.

The early development of such a theory [2] began by considering so-called
“tethered surfaces” which were modelled by continuous elastic surfaces [3, 4].
It was found that though such sheets enter a crumpled phase at large
temperatures, for sufficiently low temperatures, the thermal fluctuations
actually stabilise an undulating flat phase [5, 6]. Using renormalisation
techniques and primarily focused on the static features of such surfaces, it
was found [7, 8, 9, 10] that the thermal fluctuations effectively renormalise
the sheet’s bending rigidity with the effect that its static structure factor
becomes a length scale-dependent universal power-law in Fourier space. Similar
results were subsequently obtained for related systems such as surfaces with
intrinsic curvature, alternate geometries, inhomogeneities or exposure to
tension [11, 12, 13, 14, 15, 16, 17] and were also obtained with other methods
such as perturbative expansions in the membrane dimension [18], variational
perturbation theory [19] or phase field models [20, 21]. Simultaneously, these
theoretical predictions were investigated in experimental settings [22, 23],
including by examination of randomly crumpled surfaces [24, 25, 26, 27, 28],
and using a variety of numerical techniques such as finite-difference methods
[29, 11, 30], Monte Carlo simulations [31, 32, 33, 34, 35] and molecular
dynamics simulations [36, 37, 38, 19].

While the static properties of thermally fluctuating thin sheets have been
studied extensively, far less attention has been devoted to their dynamic
character. Early studies of the dynamic properties of tethered membranes
neglected their elastic nature [3, 4] or focused on surfaces coupled to
randomly perturbing fluids [39]. In recent years, detailed numerical studies
[40, 21] focused on the dynamics have begun to be performed but comprehensive
theoretical descriptions remain lacking.

In this paper, we study both the static and dynamic structure of thermally
fluctuating thin sheets and find them to be intimately related. Recently, we
developed a model for athermally fluctuating thin sheets with the goal of
understanding the behaviour of fluctuating sheets immediately prior to
deliberate crumpling [41, 42]. Adapting this model to provide an out-of-
equilibrium model for a thermally fluctuating thin sheet, we use a self-
consistent methodology (the self-consistent expansion or SCE) to determine the
sheet’s static and dynamic structure factors and find among other things a
previously underappreciated “quasi-linearity” inherent to the system. In
particular, for the out-of-plane displacement field $\xi(\vec{r},t)$, the
sheet’s roughness exponent $\zeta$ and dynamic exponent $z$, can be defined as
the characteristic exponents of the evolution of the variance of the
displacement field

$\left\langle\frac{1}{L^{2}}\int_{[L\times
L]}d\vec{r}\,\left[\xi\left(\vec{r},t\right)-\xi\left(\vec{r},0\right)\right]^{2}\right\rangle=L^{2\zeta}f\left(\frac{t}{L^{z}}\right)\,,$
(1)

where $L$ denotes the linear size of the sheet, $f(t/L^{z})$ denotes some
scaling function and the averaging $\langle\cdot\rangle$ denotes noise
averaging. Such dynamic scaling was extensively studied in the context of
growing surfaces and diffusion-limited aggregation where it is known as
Family-Vicsek scaling [43, 44, 45, 46]. We explicitly show that even in the
presence of strong nonlinear coupling, the apparent roughness and dynamic
exponents are related by the relationship

$z=2\zeta+2\,,$ (2)

which is also the expected relationship from a purely linear theory [47, 48].
Though such a result was previously argued for using scaling arguments [40],
here, we explicitly show that this quasi-linearity extends well beyond the
large-scale long-wavelength behaviour and is in fact maintained at all scales.

The paper is structured as follows. In section 2, we outline the Langevin
equation used to describe a thermally fluctuating thin sheet and frame it in
terms of a Fokker-Planck equation. In section 3, we analyse this model for the
case of weak nonlinear coupling using ordinary perturbation theory and thereby
obtain the static and dynamic structure factors in sections 3.1 and 3.2
respectively. In section 4, we turn to the case of strong nonlinear coupling
and, using the self-consistent expansion, derive static and dynamic structure
factors in sections 4.1 and 4.2 respectively. These analytic results are
compared to numerical simulations in section 5 and the results are discussed
in section 6.

##  2 Langevin and Fokker-Planck Equations for Fluctuating Sheets

We model fluctuating elastic thin sheets using the dynamic Föppl-von Kármán
equations, originally investigated in the context of wave turbulence [49, 50,
51, 52, 53, 54, 55, 56, 57, 58, 59, 60] and recently further developed in [41,
42] to understand the structure of thin sheets driven by fluctuations
immediately prior to crumpling. The basic idea is to consider a mostly flat
sheet and simply write a Newton force equation for each mass element of the
sheet, that is, if $\xi\left(\vec{r},t\right)$ denotes the vertical
displacement of the sheet at some point relative to a reference plane with
coordinates $\vec{r}=(x,y)$, then Newton’s second law can generically be
written as

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

where $h$ and $\rho$ denote the thickness and density of the sheet
respectively. The forces $P_{\mathrm{driving}}$ and $P_{\mathrm{damping}}$
correspond to external pressures which inject energy into the system and
dissipate it respectively while $P_{\mathrm{elastic}}$ accounts for the
restorative elastic character of the sheet which favours a flat undeformed
configuration. Following [49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 41, 42,
60], this restorative elastic pressure is taken to be determined by the Föppl-
von Kármán equations [61]

$\displaystyle P_{\mathrm{elastic}}$
$\displaystyle=\kappa\nabla^{4}\xi-h\mathcal{L}(\xi,\chi)\,,$ (4)
$\displaystyle 0$
$\displaystyle=\nabla^{4}\chi+\frac{E}{2}\mathcal{L}(\xi,\xi)\,,$ (5)

where

$\kappa=\frac{Eh^{3}}{12\left(1-\nu^{2}\right)}$ (6)

denotes the bending or flexural rigidity of the sheet, $E$ and $\nu$ denote
the sheet’s Young’s modulus and Poisson ratio respectively,
$\chi\left(\vec{r},t\right)$ denotes the Airy stress potential of the sheet
and $\mathcal{L}(f,g)$ denotes a bilinear symmetric operator which is given in
Cartesian coordinates by

$\mathcal{L}\left(f,g\right)=\frac{\partial^{2}f}{\partial
x^{2}}\frac{\partial^{2}g}{\partial y^{2}}+\frac{\partial^{2}f}{\partial
y^{2}}\frac{\partial^{2}g}{\partial x^{2}}-2\frac{\partial^{2}f}{\partial
x\partial y}\frac{\partial^{2}g}{\partial x\partial y}\,.$ (7)

The Föppl-von Kármán elastic term describes the reduction of a three-
dimensional linear-elastic thin body to a two-dimensional elastic sheet. As a
three-dimensional body, the sheet is assumed to be purely elastic. In the two-
dimensional reduction, the nonlinear terms introduced by the bilinear operator
$\mathcal{L}$ correspond to a geometric nonlinearity in which the stretching
of the material modifies the local Gaussian curvature of the sheet. This in
turn is known to have the effect of focusing stresses in the sheet thus giving
rise to the appearance of singular structures in crumpled sheets such as
$d$-cones and ridges [62, 63, 64, 65].

Various driving and damping pressures can and have been considered previously.
In the context of wave turbulence, little consideration is given to the
particular sources and sinks of energy beyond energy injection occurring at
large scales and dissipation occurring at small scales as this is sufficient
to observe an energy cascade reminiscent of hydrodynamic turbulence. A more
concrete approach is considered in [41, 42] in which the driving pressure is
taken to be conserved white noise which conserves linear momentum, while the
damping pressure is taken to be ordinary fluid friction with coefficient of
friction $\alpha$

$P_{\textrm{damping}}=\alpha\frac{\partial\xi}{\partial t}\,.$ (8)

In this work, we adapt this latter approach by considering the same type of
damping while the system is driven out-of-equilibrium by thermal noise at
temperature $T$, that is, we consider driving pressures
$P_{\textrm{driving}}=\eta\left(\vec{r},t\right)$ where

$\displaystyle\left\langle\eta\left(\vec{r},t\right)\right\rangle$
$\displaystyle=0\,,$ (9)
$\displaystyle\left\langle\eta\left(\vec{r},t\right)\eta\left(\vec{r}\,^{\prime},t^{\prime}\right)\right\rangle$
$\displaystyle=2\alpha
k_{\mathrm{B}}T\delta\left(t-t^{\prime}\right)\delta\left(\vec{r}-\vec{r}\,^{\prime}\right)\,,$
(10)

and $k_{\mathrm{B}}$ is the Boltzmann constant. Unlike the conserved noise
considered in [41, 42] which is an appropriate model when the noise originates
from forces internal to the sheet, as would occur for example with active
sheets [66], the thermal noise considered here is an appropriate model for a
sheet driven by external forces such as coupling to a thermal bath, for
example, a surrounding fluid.

For a sheet with dimensions $L\times L$, the vertical displacement
$\xi(\vec{r},t)$, Airy stress potential $\chi(\vec{r},t)$ and noise
$\eta(\vec{r},t)$ can be written in terms of their of Fourier components
$\tilde{\xi}_{\vec{n}}(t)$, $\tilde{\chi}_{\vec{n}}(t)$ and
$\tilde{\eta}_{\vec{n}}(t)$ as

$\displaystyle\xi(\vec{r},t)$
$\displaystyle=\sum_{\vec{n}\in\mathbb{Z}^{2}}\tilde{\xi}_{\vec{n}}(t)e^{i\frac{2\pi}{L}\vec{n}\cdot\vec{r}}\,,$
(11) $\displaystyle\tilde{\xi}_{\vec{n}}(t)$
$\displaystyle=\frac{1}{L^{2}}\int_{[L\times
L]}d\vec{r}\,\xi(\vec{r},t)e^{-i\frac{2\pi}{L}\vec{n}\cdot\vec{r}}$ (12)

and similarly for $\chi(\vec{r},t)$ and $\eta(\vec{r},t)$. Writing the Föppl-
von Kármán equations in terms of these Fourier components and substituting the
above relationships into equation (3) allows us to write a single nonlinear
Langevin equation for the Fourier components $\tilde{\xi}_{\vec{n}}(t)$ of the
vertical displacement. For the sake of convenience, we will only consider the
overdamped limit here such that the left-hand side of equation (3) can be
neglected and thus we obtain the nonlinear Langevin equation

$\alpha\frac{\partial\tilde{\xi}_{\vec{n}}}{\partial
t}=-\frac{\left(2\pi\right)^{4}\kappa}{L^{4}}\left|\vec{n}\right|^{4}\tilde{\xi}_{\vec{n}}-\frac{\left(2\pi\right)^{4}Eh}{2L^{4}}\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}}$
(13)

where the kernel of the sum,
$V_{\vec{n},\vec{\ell}_{1},\vec{\ell}_{2},\vec{\ell}_{3}}$, is 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}}$
(14)

and each sum is taken over all lattice points of $\mathbb{Z}^{2}$ excluding
the point $\vec{\ell}_{1}=\vec{n}$. The variance of the noise in Fourier space
is given by

$\left\langle\tilde{\eta}_{\vec{n}}\left(t\right)\tilde{\eta}_{\vec{n}^{\prime}}\left(t^{\prime}\right)\right\rangle=\frac{2\alpha
k_{\mathrm{B}}T}{L^{2}}\delta\left(t-t^{\prime}\right)\delta_{\vec{n},-\vec{n}^{\prime}}\,.$
(15)

Formally, equation (13) is a nonlinear Langevin equation with a quartic
interaction. Accordingly, it can be thought of as a type of $\phi^{4}$-field
theory [67] with a spatially varying kernel
$V_{\vec{n},\vec{\ell}_{1},\vec{\ell}_{2},\vec{\ell}_{3}}$. Importantly, the
behaviour of the zeroth mode $\tilde{\xi}_{0}$ is completely decoupled from
that of the other modes since subbing $\vec{n}=0$ into equation (13) gives

$\alpha\frac{\partial\tilde{\xi}_{0}}{\partial t}=\tilde{\eta}_{0}(t)$ (16)

while the kernel $V_{\vec{n},\vec{\ell}_{1},\vec{\ell}_{2},\vec{\ell}_{3}}$
vanishes whenever $\vec{\ell}_{1}$, $\vec{\ell}_{2}$ or $\vec{\ell}_{3}$ equal
0 thus ensuring that for all $\vec{n}\neq 0$, the mode $\tilde{\xi}_{0}$
provides no contribution to the Langevin equation. This has the strong
implication that the mean position of the sheet, which is given by the zeroth
mode $\tilde{\xi}_{0}$, simply performs a random walk in space with no effect
on the structure of the sheet and the sums over $\vec{\ell}_{1}$,
$\vec{\ell}_{2}$ or $\vec{\ell}_{3}$ can be taken as skipping over the origin.

The structure of a sheet undergoing thermal fluctuations described by equation
(13) is most easily approached by studying its corresponding Fokker-Planck
equation for the probability $P=P(\\{\tilde{\xi}_{\vec{n}}(t)\\},t)$ of
finding the sheet with Fourier components $\\{\tilde{\xi}_{\vec{n}}(t)\\}$ at
time $t$ [68]

$\frac{\partial P}{\partial t}=\frac{k_{\mathrm{B}}T}{\alpha
L^{2}}\sum_{\vec{n}}\frac{\partial^{2}P}{\partial\tilde{\xi}_{\vec{n}}\partial\tilde{\xi}_{-\vec{n}}}+\frac{\left(2\pi\right)^{4}\kappa}{\alpha
L^{4}}\sum_{\vec{n}}\left|\vec{n}\right|^{4}\frac{\partial}{\partial\tilde{\xi}_{\vec{n}}}\left[\tilde{\xi}_{\vec{n}}P\right]\\\
+\frac{\left(2\pi\right)^{4}Eh}{2\alpha
L^{4}}\sum_{\vec{n}}\frac{\partial}{\partial\tilde{\xi}_{\vec{n}}}\Biggl{[}P\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}}\tilde{\xi}_{\vec{\ell}_{1}}\tilde{\xi}_{\vec{\ell}_{2}}\tilde{\xi}_{\vec{\ell}_{3}}\Biggr{]}\,.$
(17)

Multiplying this equation by any function
$\mathbb{F}(\\{\tilde{\xi}_{\vec{n}}(t)\\})$ of the Fourier components and
integrating over all components $\tilde{\xi}_{\vec{n}}(t)$ gives after some
integration by parts

$\frac{\partial\left\langle\mathbb{F}\right\rangle}{\partial
t}=\frac{k_{\mathrm{B}}T}{\alpha
L^{2}}\sum_{\vec{n}}\left\langle\frac{\partial^{2}\mathbb{F}}{\partial\tilde{\xi}_{\vec{n}}\partial\tilde{\xi}_{-\vec{n}}}\right\rangle-\frac{\left(2\pi\right)^{4}\kappa}{\alpha
L^{4}}\sum_{\vec{n}}\left|\vec{n}\right|^{4}\left\langle\frac{\partial\mathbb{F}}{\partial\tilde{\xi}_{\vec{n}}}\tilde{\xi}_{\vec{n}}\right\rangle\\\
-\frac{\left(2\pi\right)^{4}Eh}{2\alpha
L^{4}}\sum_{\vec{n}}\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\frac{\partial\mathbb{F}}{\partial\tilde{\xi}_{\vec{n}}}\tilde{\xi}_{\vec{\ell}_{1}}\tilde{\xi}_{\vec{\ell}_{2}}\tilde{\xi}_{\vec{\ell}_{3}}\right\rangle\,,$
(18)

where we have defined the averages

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

The first term on the right-hand side of equations (17) and (18) corresponds
to diffusion through the phase space of the fluctuating sheet due to the
thermal noise while the second term corresponds to the limitation of this
diffusion due to the restorative effect of the sheet’s linear elasticity. The
final term of these equations modifies this restorative effect by accounting
for the geometric nonlinearity of a deformed sheet as described by the
nonlinear portion of the Föppl-von Kármán equations (equations (4) and (5))
and its corresponding kernel in Fourier space
$V_{\vec{n},\vec{\ell}_{1},\vec{\ell}_{2},\vec{\ell}_{3}}$. In instances where
the effect of the geometric nonlinearity is weak, the final term can be
neglected at zeroth order and a linear theory is obtained which in principle
can be perturbed around. In contrast, when the effect of the nonlinearity is
large, the linear theory fails to provide any sort of reasonable approximation
of the sheet and more sophisticated techniques of analysis must be employed.
In the next sections, we will argue that both of these regimes are relevant,
interesting and worth addressing. We will begin by analysing the behaviour of
the sheet in its linear regime before moving onto its nonlinear regime.

##  3 The Linear Regime

To identify the linear regime, it is helpful to begin by nondimensionalising
our equations. Examining the physical parameters which enter our model, we
find that one can define the dimensionless quantity

$g=\frac{\left(2\pi\right)^{2}\kappa}{\sqrt{2k_{\mathrm{B}}TEhL^{2}}}$ (20)

where the size of $g$ will determine whether we are in the linear or nonlinear
regime. If we define dimensionless scalings of the time $\hat{t}$ and vertical
displacement $\hat{\xi}$ as follows,

$\displaystyle\hat{t}$
$\displaystyle=\frac{\left(2\pi\right)^{4}\kappa}{\alpha L^{4}}t\,,$ (21)
$\displaystyle\hat{\xi}_{\vec{n}}\left(\hat{t}\,\right)$
$\displaystyle=\left[\frac{\left(2\pi\right)^{4}\kappa}{2k_{\mathrm{B}}TL^{2}}\right]^{1/2}\tilde{\xi}_{\vec{n}}\left(t\right)\,,$
(22)

then equation (18) can be written as

$\frac{\partial\left\langle\mathbb{F}\right\rangle}{\partial\hat{t}}=\frac{1}{2}\sum_{\vec{n}}\left\langle\frac{\partial^{2}\mathbb{F}}{\partial\hat{\xi}_{\vec{n}}\partial\hat{\xi}_{-\vec{n}}}\right\rangle-\sum_{\vec{n}}\left|\vec{n}\right|^{4}\left\langle\frac{\partial\mathbb{F}}{\partial\hat{\xi}_{\vec{n}}}\hat{\xi}_{\vec{n}}\right\rangle\\\
-\frac{1}{2}\frac{1}{g^{2}}\sum_{\vec{n}}\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\frac{\partial\mathbb{F}}{\partial\hat{\xi}_{\vec{n}}}\hat{\xi}_{\vec{\ell}_{1}}\hat{\xi}_{\vec{\ell}_{2}}\hat{\xi}_{\vec{\ell}_{3}}\right\rangle\,,$
(23)

which contains just the single dimensionless parameter $g$ and thus we find
that we will be in the linear regime if $g\gg 1$.

To understand the physical significance of the dimensionless parameter $g$, it
is helpful to recast equation (20) in terms which make the dependence on the
physical size explicit. In particular, it is known from equation (6) that the
flexural rigidity $\kappa$ grows with the sheet thickness as $\sim h^{3}$.
Accordingly, in terms of the Young’s modulus $E$ and Poisson ratio $\nu$,
which are both independent of sheet thickness, $g$ can be written as

$g=\frac{\left(2\pi\right)^{2}}{12(1-\nu^{2})}\sqrt{\frac{Eh^{5}}{2k_{\mathrm{B}}TL^{2}}}\,.$
(24)

This expression makes it clear that $g$ grows with sheet thickness and thus
nonlinear effects will become more important for thinner sheets.

Physically, the condition $g\gg 1$ occurs for macroscopic sheets at ordinary
temperatures. For instance, a sheet of aluminium with length $L\sim 20\textrm{
cm}$ to a side, would have a typical thickness around $h\sim 0.02\textrm{
mm}$, a Young’s modulus around $E\sim 70\textrm{ GPa}$, a Poisson ratio of
$\nu\sim 0.35$ and thus a bending rigidity $\kappa\sim 5\times 10^{-5}\textrm{
N m}$. Accordingly, at temperatures of say $T\sim 300\textrm{ K}$, we would
have $g\sim 10^{5}$ which is indeed much much larger than unity.

For such sheets, the Föppl-von Kármán geometric nonlinearity can be completely
neglected and if we so desire, a perturbative expansion can be carried out
around the linear part. To this end, let
$\left\langle\mathbb{F}\right\rangle^{(m)}$ denote an $m^{\textrm{th}}$ order
approximation for any quantity $\left\langle\mathbb{F}\right\rangle$, that is,
$\left\langle\mathbb{F}\right\rangle^{(m)}$ denotes an expansion containing
all corrections up to order $m$. Then the following iterative scheme

$\frac{\partial\left\langle\mathbb{F}\right\rangle^{\left(m\right)}}{\partial\hat{t}}=\frac{1}{2}\sum_{\vec{n}}\left\langle\frac{\partial^{2}\mathbb{F}}{\partial\hat{\xi}_{\vec{n}}\partial\hat{\xi}_{-\vec{n}}}\right\rangle^{\left(m\right)}-\sum_{\vec{n}}\left|\vec{n}\right|^{4}\left\langle\frac{\partial\mathbb{F}}{\partial\hat{\xi}_{\vec{n}}}\hat{\xi}_{\vec{n}}\right\rangle^{\left(m\right)}\\\
-\frac{1}{2}\frac{1}{g^{2}}\sum_{\vec{n}}\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\frac{\partial\mathbb{F}}{\partial\hat{\xi}_{\vec{n}}}\hat{\xi}_{\vec{\ell}_{1}}\hat{\xi}_{\vec{\ell}_{2}}\hat{\xi}_{\vec{\ell}_{3}}\right\rangle^{\left(m-1\right)}$
(25)

can be used to calculate any expectation value
$\left\langle\mathbb{F}\right\rangle$ order by order. In this notation,
$\left\langle\mathbb{F}\right\rangle^{(0)}$ corresponds to the solution to the
purely linear theory and since the small parameter of equation (25) is
$1/g^{2}$, $\left\langle\mathbb{F}\right\rangle^{(m)}$ will contain a term
proportional to $1/g^{2m}$. In the following subsections, we demonstrate this
approach by calculating the static and dynamic structure factors up to first
order in the linear regime.

###  3.1 The Static Structure Factor

The static structure factor
$\langle\hat{\xi}_{\vec{n}}\hat{\xi}_{-\vec{n}}\rangle$ can be calculated at
zeroth order by subbing
$\mathbb{F}=\hat{\xi}_{\vec{n}_{1}}\hat{\xi}_{\vec{n}_{2}}$ and $m=0$ into
equation (25)

$0=\frac{1}{2}\sum_{\vec{n}}\left\langle\frac{\partial^{2}\left(\hat{\xi}_{\vec{n}_{1}}\hat{\xi}_{\vec{n}_{2}}\right)}{\partial\hat{\xi}_{\vec{n}}\partial\hat{\xi}_{-\vec{n}}}\right\rangle^{\left(0\right)}-\sum_{\vec{n}}\left|\vec{n}\right|^{4}\left\langle\frac{\partial\left(\hat{\xi}_{\vec{n}_{1}}\hat{\xi}_{\vec{n}_{2}}\right)}{\partial\hat{\xi}_{\vec{n}}}\hat{\xi}_{\vec{n}}\right\rangle^{\left(0\right)}\,.$
(26)

Note that the left-hand side has been set to 0 since in steady-state, the
equal-time two-point function
$\langle\hat{\xi}_{\vec{n}_{1}}\hat{\xi}_{\vec{n}_{2}}\rangle$ does not depend
on time. Calculating the derivatives and simplifying the sums, we obtain the
static two-point function

$\left\langle\hat{\xi}_{\vec{n}_{1}}\hat{\xi}_{\vec{n}_{2}}\right\rangle^{\left(0\right)}=\frac{\delta_{\vec{n}_{1},-\vec{n}_{2}}}{\left|\vec{n}_{1}\right|^{4}+\left|\vec{n}_{2}\right|^{4}}=\frac{\delta_{\vec{n}_{1},-\vec{n}_{2}}}{2\left|\vec{n}_{1}\right|^{4}}\,.$
(27)

If we wish to consider macroscopic sheets being driven by white noise of
arbitrary amplitude, i.e. by effective temperatures much larger than room
temperature and thus with smaller corresponding values of $g$, corrections to
this expression may become relevant. To this end, subbing
$\mathbb{F}=\hat{\xi}_{\vec{n}_{1}}\hat{\xi}_{\vec{n}_{2}}$ and $m=1$ into
equation (25) gives, after some simplification, the equation

$0=\delta_{\vec{n}_{1},-\vec{n}_{2}}-\left(\left|\vec{n}_{1}\right|^{4}+\left|\vec{n}_{2}\right|^{4}\right)\left\langle\hat{\xi}_{\vec{n}_{1}}\hat{\xi}_{\vec{n}_{2}}\right\rangle^{\left(1\right)}\\\
-\frac{1}{2}\frac{1}{g^{2}}\Biggl{[}\,\sum_{\begin{subarray}{c}\vec{\ell}_{1}\neq\vec{n}_{1}\\\
\vec{\ell}_{2},\vec{\ell}_{3}\end{subarray}}V_{\vec{n}_{1},\vec{\ell}_{1},\vec{\ell}_{2},\vec{\ell}_{3}}\left\langle\hat{\xi}_{\vec{n}_{2}}\hat{\xi}_{\vec{\ell}_{1}}\hat{\xi}_{\vec{\ell}_{2}}\hat{\xi}_{\vec{\ell}_{3}}\right\rangle^{\left(0\right)}+\sum_{\begin{subarray}{c}\vec{\ell}_{1}\neq\vec{n}_{2}\\\
\vec{\ell}_{2},\vec{\ell}_{3}\end{subarray}}V_{\vec{n}_{2},\vec{\ell}_{1},\vec{\ell}_{2},\vec{\ell}_{3}}\left\langle\hat{\xi}_{\vec{n}_{1}}\hat{\xi}_{\vec{\ell}_{1}}\hat{\xi}_{\vec{\ell}_{2}}\hat{\xi}_{\vec{\ell}_{3}}\right\rangle^{\left(0\right)}\Biggr{]}\,.$
(28)

To proceed, the zeroth order equal-time four-point function
$\langle\hat{\xi}_{\vec{n}_{1}}\hat{\xi}_{\vec{\ell}_{1}}\hat{\xi}_{\vec{\ell}_{2}}\hat{\xi}_{\vec{\ell}_{3}}\rangle^{\left(0\right)}$
is required. This can be obtained by subbing
$\mathbb{F}=\hat{\xi}_{\vec{n}_{1}}\hat{\xi}_{\vec{\ell}_{1}}\hat{\xi}_{\vec{\ell}_{2}}\hat{\xi}_{\vec{\ell}_{3}}$
and $m=0$ into equation (25), which after some straightforward manipulation
results in

$\left\langle\hat{\xi}_{\vec{n}_{1}}\hat{\xi}_{\vec{\ell}_{1}}\hat{\xi}_{\vec{\ell}_{2}}\hat{\xi}_{\vec{\ell}_{3}}\right\rangle^{\left(0\right)}=\left\langle\hat{\xi}_{\vec{n}_{1}}\hat{\xi}_{\vec{\ell}_{1}}\right\rangle^{\left(0\right)}\left\langle\hat{\xi}_{\vec{\ell}_{2}}\hat{\xi}_{\vec{\ell}_{3}}\right\rangle^{\left(0\right)}\\\
+\left\langle\hat{\xi}_{\vec{n}_{1}}\hat{\xi}_{\vec{\ell}_{2}}\right\rangle^{\left(0\right)}\left\langle\hat{\xi}_{\vec{\ell}_{1}}\hat{\xi}_{\vec{\ell}_{3}}\right\rangle^{\left(0\right)}+\left\langle\hat{\xi}_{\vec{n}_{1}}\hat{\xi}_{\vec{\ell}_{3}}\right\rangle^{\left(0\right)}\left\langle\hat{\xi}_{\vec{\ell}_{1}}\hat{\xi}_{\vec{\ell}_{2}}\right\rangle^{\left(0\right)}$
(29)

This is of course none other than Isserlis’ theorem [69, 70], also known as
Wick’s theorem [71], for the four-point function of a multivariate Gaussian
distribution. Accordingly, combining this result and the definition of
$V_{\vec{n}_{2},\vec{\ell}_{1},\vec{\ell}_{2},\vec{\ell}_{3}}$ given by
equation (14) allows us to obtain

$\sum_{\begin{subarray}{c}\vec{\ell}_{1}\neq\vec{n}_{2}\\\
\vec{\ell}_{2},\vec{\ell}_{3}\end{subarray}}V_{\vec{n}_{2},\vec{\ell}_{1},\vec{\ell}_{2},\vec{\ell}_{3}}\left\langle\hat{\xi}_{\vec{n}_{1}}\hat{\xi}_{\vec{\ell}_{1}}\hat{\xi}_{\vec{\ell}_{2}}\hat{\xi}_{\vec{\ell}_{3}}\right\rangle^{\left(0\right)}=\left\langle\hat{\xi}_{\vec{n}_{1}}\hat{\xi}_{\vec{n}_{2}}\right\rangle^{\left(0\right)}\sum_{\vec{\ell}_{1}\neq\vec{n}_{2}}\frac{\left|\vec{n}_{2}\times\vec{\ell}_{1}\right|^{4}}{\left|\vec{\ell}_{1}\right|^{4}\left|\vec{n}_{2}-\vec{\ell}_{1}\right|^{4}}$
(30)

and thus equation (28) can be rearranged to give the static structure factor
as

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

Though the summand in this expression is ill defined when $\vec{\ell}=(0,0)$,
recall from section 2 that no harm is done when the sums skip over the origin
since the zeroth mode $\tilde{\xi}_{0}$ is completely decoupled from the other
modes anyway. The sum can be quickly evaluated by approximating it as the
integral

$\sum_{\vec{\ell}\neq
0,\vec{n}}\frac{\left|\vec{n}\times\vec{\ell}\,\right|^{4}}{\left|\vec{\ell}\,\right|^{4}\left|\vec{n}-\vec{\ell}\,\right|^{4}}\approx\int
d^{2}\ell\,\frac{\left|\vec{n}\times\vec{\ell}\,\right|^{4}}{\left|\vec{\ell}\,\right|^{4}\left|\vec{n}-\vec{\ell}\,\right|^{4}}\,.$
(32)

Noting that
$|\vec{n}\times\vec{\ell}\,|^{4}=|\vec{n}|^{4}|\vec{\ell}\,|^{4}\sin^{4}\theta$
where $\theta$ is the angle between $\vec{n}$ and $\vec{\ell}$ allows this
integral to be calculated in polar coordinates. Writing the integral as

$\int
d^{2}\ell\,\frac{\left|\vec{n}\times\vec{\ell}\,\right|^{4}}{\left|\vec{\ell}\,\right|^{4}\left|\vec{n}-\vec{\ell}\,\right|^{4}}=\\\
n^{4}\int_{0}^{\infty}d\ell\,\ell\int_{0}^{2\pi}d\theta\frac{\sin^{4}\theta}{\left(n^{2}+\ell^{2}-2n\ell\cos\theta\right)^{2}}$
(33)

and noting that the angular integral contributes

$\int_{0}^{2\pi}d\theta\,\frac{\sin^{4}\theta}{\left(n^{2}+\ell^{2}-2n\ell\cos\theta\right)^{2}}=\frac{3\pi}{4}\begin{cases}n^{-4}&\ell\leq
n\\\ \ell^{-4}&\ell\geq n\end{cases}\,,$ (34)

we find that

$\int
d^{2}\ell\,\frac{\left|\vec{n}\times\vec{\ell}\,\right|^{4}}{\left|\vec{\ell}\,\right|^{4}\left|\vec{n}-\vec{\ell}\,\right|^{4}}=\frac{3\pi}{4}n^{2}\,.$
(35)

Accordingly, for large $g$, the static structure factor
$\langle\hat{\xi}_{\vec{n}}\hat{\xi}_{-\vec{n}}\rangle$ is given up to first
order by

$\left\langle\hat{\xi}_{\vec{n}}\hat{\xi}_{-\vec{n}}\right\rangle^{\left(1\right)}=\frac{1}{2n^{4}}\left(1-\frac{3\pi}{8g^{2}}\frac{1}{n^{2}}\right)\,.$
(36)

Equivalently, we can rewrite this as the straightforward Padé approximant [72]

$\left\langle\hat{\xi}_{\vec{n}}\hat{\xi}_{-\vec{n}}\right\rangle^{\left(1\right)}=\frac{1}{2n^{4}\left(1+\frac{3\pi}{8g^{2}}\frac{1}{n^{2}}\right)}\,,$
(37)

which makes it clear that the first order correction is simply a
renormalisation of the linear coupling $n^{4}$ by
$n^{4}\left(1+\frac{3\pi}{8g^{2}}\frac{1}{n^{2}}\right)$. Calculating higher
order corrections is of course possible though tedious.

Returning to dimensional quantities, we find that

$\left\langle\tilde{\xi}_{\vec{n}}\tilde{\xi}_{-\vec{n}}\right\rangle^{\left(1\right)}=\frac{k_{\mathrm{B}}TL^{2}}{\left(2\pi\right)^{4}\kappa
n^{4}\left(1+\frac{3\pi}{8\left(2\pi\right)^{4}}\frac{2k_{\mathrm{B}}TEhL^{2}}{\kappa^{2}}\frac{1}{n^{2}}\right)}\,.$
(38)

The scaling picture described by equation (1) can be translated to Fourier
space [73] where it implies that the structure factor decays as

$\left\langle\tilde{\xi}_{\vec{n}}\tilde{\xi}_{-\vec{n}}\right\rangle\sim\frac{1}{n^{d+2\zeta}}\,,$
(39)

where $\zeta$ denotes the roughness exponent of the sheet and $d$ its
dimension. We have $d=2$ and thus comparing this definition of the roughness
exponent with the result presented by equation (38), it is clear that the
roughness exponent of the linear regime is simply $\zeta=1$ though this only
characterises the very large scale behaviour of the sheet.

Finally, when addressing the nonlinear regime, it will become apparent that,
under certain conditions, systems which have an upper frequency cut-off can
behave differently from those which don’t. Accordingly, it is worth
considering the effect of an upper cut-off already at this stage.
Conceptually, suppose that at sufficiently fine scales, the system has some
minimum length scale $a\ll L$ beneath which the Föppl-von Kármán equations no
longer apply. This introduces a new dimensionless quantity into our system,
$n_{\mathrm{max}}=L/(2a)\gg 1$, which enters our calculations as an upper
bound on the sum and integral given in equation (32). Replacing the $\infty$
in the upper-bound of equation (33) with the upper frequency cut-off
$n_{\mathrm{max}}$, we ultimately find that

$\left\langle\hat{\xi}_{\vec{n}}\hat{\xi}_{-\vec{n}}\right\rangle^{\left(1\right)}=\frac{1}{2n^{4}\left(1+\frac{3\pi}{8g^{2}}\frac{1}{n^{2}}\left[1-\frac{1}{2}\left(\frac{n}{n_{\mathrm{max}}}\right)^{2}\right]\right)}\,.$
(40)

That is, the existence of an upper frequency cut-off has the effect of
slightly modifying the first correction to the renormalised linear coupling.
Notice that in the limit of $n_{\mathrm{max}}\rightarrow\infty$, we simply
return to the result described by equation (37) and thus the upper cut-off has
little effect in the linear regime. Furthermore, for small $n$ (i.e. $n\sim
1$), the effect of the upper cut-off is extremely small and thus the
expression given by equation (37) is perfectly sufficient, while for large $n$
(i.e. $n\sim n_{\mathrm{max}}$), the effect of the upper cut-off on the first
order correction might indeed be large but since the entire first order
correction becomes increasingly small for large $n$, the effect of the upper
cut-off is again negligible. We thus find that in the linear regime, an upper
frequency cut-off can be safely neglected.

###  3.2 The Dynamic Structure Factor

Equation (25) can also be used to calculate dynamic quantities such as the
dynamic structure factor
$\langle\hat{\xi}_{\vec{n}}(0)\hat{\xi}_{-\vec{n}}(\hat{t}\,)\rangle$, though
first, we should clarify a notational ambiguity. The derivatives in equation
(25) are with respect to $\hat{\xi}_{\vec{n}}=\hat{\xi}_{\vec{n}}(\hat{t}\,)$,
that is, for the sake of conciseness, the time dependence has been implicitly
understood to exist, even if it is not noted explicitly. Until now, this made
little difference as the static quantities considered were independent of time
such that
$\langle\hat{\xi}_{\vec{n}}(\hat{t}\,)\hat{\xi}_{-\vec{n}}(\hat{t}\,)\rangle=\langle\hat{\xi}_{\vec{n}}(0)\hat{\xi}_{-\vec{n}}(0)\rangle$
and similarly for the static four-point function. For time-dependent
quantities however, the time dependence cannot be ignored thus while we will
continue to use this notation, we will need to bear in mind that if
unspecified, $\hat{\xi}_{\vec{n}}$ is merely shorthand for
$\hat{\xi}_{\vec{n}}(\hat{t}\,)$. If some other time such as
$\hat{\xi}_{\vec{n}}(0)$ is intended, this will always be stated explicitly.

Now, subbing in
$\mathbb{F}=\hat{\xi}_{\vec{n}_{1}}\left(0\right)\hat{\xi}_{\vec{n}_{2}}(\hat{t}\,)$
and $m=0$ into equation (25) results in the ODE

$\frac{\partial}{\partial\hat{t}}\left\langle\hat{\xi}_{\vec{n}_{1}}\left(0\right)\hat{\xi}_{\vec{n}_{2}}(\hat{t}\,)\right\rangle^{\left(0\right)}=-\left|\vec{n}_{2}\right|^{4}\left\langle\hat{\xi}_{\vec{n}_{1}}\left(0\right)\hat{\xi}_{\vec{n}_{2}}(\hat{t}\,)\right\rangle^{\left(0\right)}\,,$
(41)

which is simply solved by

$\left\langle\hat{\xi}_{\vec{n}_{1}}\left(0\right)\hat{\xi}_{\vec{n}_{2}}(\hat{t}\,)\right\rangle^{\left(0\right)}=\left\langle\hat{\xi}_{\vec{n}_{1}}\hat{\xi}_{\vec{n}_{2}}\right\rangle^{\left(0\right)}e^{-|\vec{n}_{2}|^{4}\hat{t}}\,.$
(42)

This exponential decay is of course the expected result for a completely
linear theory. As with the static structure factor, corrections can be
obtained by subbing
$\mathbb{F}=\hat{\xi}_{\vec{n}_{1}}\left(0\right)\hat{\xi}_{\vec{n}_{2}}(\hat{t}\,)$
with $m=1$ into equation (25)

$\frac{\partial}{\partial\hat{t}}\left\langle\hat{\xi}_{\vec{n}_{1}}\left(0\right)\hat{\xi}_{\vec{n}_{2}}(\hat{t}\,)\right\rangle^{\left(1\right)}=-\left|\vec{n}_{2}\right|^{4}\left\langle\hat{\xi}_{\vec{n}_{1}}\left(0\right)\hat{\xi}_{\vec{n}_{2}}\right\rangle^{\left(1\right)}\\\
-\frac{1}{2}\frac{1}{g^{2}}\sum_{\begin{subarray}{c}\vec{\ell}_{1}\neq\vec{n}_{2}\\\
\vec{\ell}_{2},\vec{\ell}_{3}\end{subarray}}V_{\vec{n}_{2},\vec{\ell}_{1},\vec{\ell}_{2},\vec{\ell}_{3}}\left\langle\hat{\xi}_{\vec{n}_{1}}\left(0\right)\hat{\xi}_{\vec{\ell}_{1}}\hat{\xi}_{\vec{\ell}_{2}}\hat{\xi}_{\vec{\ell}_{3}}\right\rangle^{\left(0\right)}$
(43)

and as before, we find that we need the zeroth order time-dependent four-point
function to proceed. Subbing
$\mathbb{F}=\hat{\xi}_{\vec{n}_{1}}\left(0\right)\hat{\xi}_{\vec{\ell}_{1}}\hat{\xi}_{\vec{\ell}_{2}}\hat{\xi}_{\vec{\ell}_{3}}$
and $m=0$ into equation (25) gives the ODE

$\frac{\partial}{\partial\hat{t}}\left\langle\hat{\xi}_{\vec{n}_{1}}\left(0\right)\hat{\xi}_{\vec{\ell}_{1}}\hat{\xi}_{\vec{\ell}_{2}}\hat{\xi}_{\vec{\ell}_{3}}\right\rangle^{\left(0\right)}=-\left(\left|\vec{\ell}_{1}\right|^{4}+\left|\vec{\ell}_{2}\right|^{4}+\left|\vec{\ell}_{3}\right|^{4}\right)\left\langle\hat{\xi}_{\vec{n}_{1}}\left(0\right)\hat{\xi}_{\vec{\ell}_{1}}\hat{\xi}_{\vec{\ell}_{2}}\hat{\xi}_{\vec{\ell}_{3}}\right\rangle^{\left(0\right)}\\\
+\delta_{\vec{\ell}_{1},-\vec{\ell}_{2}}\left\langle\hat{\xi}_{\vec{n}_{1}}\left(0\right)\hat{\xi}_{\vec{\ell}_{3}}\right\rangle^{\left(0\right)}+\delta_{\vec{\ell}_{1},-\vec{\ell}_{3}}\left\langle\hat{\xi}_{\vec{n}_{1}}\left(0\right)\hat{\xi}_{\vec{\ell}_{2}}\right\rangle^{\left(0\right)}+\delta_{\vec{\ell}_{2},-\vec{\ell}_{3}}\left\langle\hat{\xi}_{\vec{n}_{1}}\left(0\right)\hat{\xi}_{\vec{\ell}_{1}}\right\rangle^{\left(0\right)}\,.$
(44)

Focusing on any of the last three terms of this equation, note that equations
(27) and (42) can be used to write

$\delta_{\vec{\ell}_{1},-\vec{\ell}_{2}}\left\langle\hat{\xi}_{\vec{n}_{1}}\left(0\right)\hat{\xi}_{\vec{\ell}_{3}}\right\rangle^{\left(0\right)}=\\\
\left(\left|\vec{\ell}_{1}\right|^{4}+\left|\vec{\ell}_{2}\right|^{4}+\left|\vec{\ell}_{3}\right|^{4}-\left|\vec{n}_{1}\right|^{4}\right)\left\langle\hat{\xi}_{\vec{\ell}_{1}}\hat{\xi}_{\vec{\ell}_{2}}\right\rangle^{\left(0\right)}\left\langle\hat{\xi}_{\vec{n}_{1}}\hat{\xi}_{\vec{\ell}_{3}}\right\rangle^{\left(0\right)}e^{-|\vec{n}_{1}|^{4}\hat{t}}\,,$
(45)

where we have added and subtracted $\left|\vec{\ell}_{3}\right|^{4}$ and
$\left|\vec{n}_{1}\right|^{4}$ since the two point function
$\left\langle\hat{\xi}_{\vec{n}_{1}}\hat{\xi}_{\vec{\ell}_{3}}\right\rangle^{\left(0\right)}$
anyway vanishes unless $\vec{\ell}_{3}=-\vec{n}_{1}$. Accordingly, upon
comparison with equation (29), the sum of the last three terms of equation
(44) can be written as

$\delta_{\vec{\ell}_{1},-\vec{\ell}_{2}}\left\langle\hat{\xi}_{\vec{n}_{1}}\left(0\right)\hat{\xi}_{\vec{\ell}_{3}}\right\rangle^{\left(0\right)}+\delta_{\vec{\ell}_{1},-\vec{\ell}_{3}}\left\langle\hat{\xi}_{\vec{n}_{1}}\left(0\right)\hat{\xi}_{\vec{\ell}_{2}}\right\rangle^{\left(0\right)}+\delta_{\vec{\ell}_{2},-\vec{\ell}_{3}}\left\langle\hat{\xi}_{\vec{n}_{1}}\left(0\right)\hat{\xi}_{\vec{\ell}_{1}}\right\rangle^{\left(0\right)}=\\\
\left(\left|\vec{\ell}_{1}\right|^{4}+\left|\vec{\ell}_{2}\right|^{4}+\left|\vec{\ell}_{3}\right|^{4}-\left|\vec{n}_{1}\right|^{4}\right)\left\langle\hat{\xi}_{\vec{n}_{1}}\hat{\xi}_{\vec{\ell}_{1}}\hat{\xi}_{\vec{\ell}_{2}}\hat{\xi}_{\vec{\ell}_{3}}\right\rangle^{\left(0\right)}e^{-|\vec{n}_{1}|^{4}\hat{t}}\,,$
(46)

and thus it becomes easy to check that the solution to equation (44) is
unsurprisingly

$\left\langle\hat{\xi}_{\vec{n}_{1}}\left(0\right)\hat{\xi}_{\vec{\ell}_{1}}\hat{\xi}_{\vec{\ell}_{2}}\hat{\xi}_{\vec{\ell}_{3}}\right\rangle^{\left(0\right)}=\left\langle\hat{\xi}_{\vec{n}_{1}}\hat{\xi}_{\vec{\ell}_{1}}\hat{\xi}_{\vec{\ell}_{2}}\hat{\xi}_{\vec{\ell}_{3}}\right\rangle^{\left(0\right)}e^{-|\vec{n}_{1}|^{4}\hat{t}}\,.$
(47)

Returning to equation (43) and making use of equation (30), we obtain the
following nonhomogeneous ODE for the first order time-dependent two-point
function

$\frac{\partial}{\partial\hat{t}}\left\langle\hat{\xi}_{\vec{n}_{1}}\left(0\right)\hat{\xi}_{\vec{n}_{2}}(\hat{t}\,)\right\rangle^{\left(1\right)}=-\left|\vec{n}_{2}\right|^{4}\left\langle\hat{\xi}_{\vec{n}_{1}}\left(0\right)\hat{\xi}_{\vec{n}_{2}}(\hat{t}\,)\right\rangle^{\left(1\right)}\\\
-\frac{3\pi}{8g^{2}}\left|\vec{n}_{2}\right|^{2}\left\langle\hat{\xi}_{\vec{n}_{1}}\hat{\xi}_{\vec{n}_{2}}\right\rangle^{\left(0\right)}e^{-|\vec{n}_{2}|^{4}\hat{t}}\,.$
(48)

Here, we have freely replaced the decay rate of the exponential
$|\vec{n}_{1}|^{4}$ with $|\vec{n}_{2}|^{4}$ since the static two-point
function ensures that this term vanishes anyway if
$\vec{n}_{1}\neq-\vec{n}_{2}$. We have also replaced the sum, now given by
equation (30), with its approximate value given by equation (35).

Equation (48) is easily solved but the fact that its nonhomogeneous part
decays with the decay rate of the homogeneous part results in the solution
containing a secular term

$\left\langle\hat{\xi}_{\vec{n}_{1}}\left(0\right)\hat{\xi}_{\vec{n}_{2}}(\hat{t}\,)\right\rangle^{\left(1\right)}=\\\
\left[\left\langle\hat{\xi}_{\vec{n}_{1}}\hat{\xi}_{\vec{n}_{2}}\right\rangle^{\left(1\right)}-\frac{3\pi}{8g^{2}}\left|\vec{n}_{2}\right|^{2}\left\langle\hat{\xi}_{\vec{n}_{1}}\hat{\xi}_{\vec{n}_{2}}\right\rangle^{\left(0\right)}\hat{t}\,\right]e^{-|\vec{n}_{2}|^{4}\hat{t}}\,.$
(49)

This secular term is an unphysical artefact of the perturbation theory which
limits the solution to short times. If we wish to obtain a first order
correction which is valid for all times, we need to perform a more
sophisticated perturbation theory, such as the Poincaré-Lindstedt method [74,
75, 76], in which the decay rate of equation (25) is also expanded order by
order. While such an analysis is possible, it is faster to simply note that
equation (49) is merely the large $g$ expansion of

$\left\langle\hat{\xi}_{\vec{n}_{1}}\left(0\right)\hat{\xi}_{\vec{n}_{2}}(\hat{t}\,)\right\rangle^{\left(1\right)}=\\\
\left\langle\hat{\xi}_{\vec{n}_{1}}\hat{\xi}_{\vec{n}_{2}}\right\rangle^{\left(1\right)}e^{-\left|\vec{n}_{2}\right|^{4}\left(1+\frac{3\pi}{8g^{2}}\frac{1}{\left|\vec{n}_{2}\right|^{2}}\right)\hat{t}}$
(50)

where we have used in the exponential the renormalised linear coupling implied
by equation (37). Accordingly, we have found that in the linear regime, the
weakly coupled nonlinearity merely has the effect of renormalising the linear
coupling $n^{4}\rightarrow
n^{4}\left(1+\frac{3\pi}{8g^{2}}\frac{1}{n^{2}}\right)$, both for static and
dynamic quantities.

In dimensional terms, we obtain

$\left\langle\tilde{\xi}_{\vec{n}_{1}}\left(0\right)\tilde{\xi}_{\vec{n}_{2}}(t\,)\right\rangle^{\left(1\right)}=\\\
\left\langle\tilde{\xi}_{\vec{n}_{1}}\tilde{\xi}_{\vec{n}_{2}}\right\rangle^{\left(1\right)}e^{-\left(2\pi\right)^{4}\frac{\kappa\left|\vec{n}_{2}\right|^{4}}{\alpha
L^{4}}\left(1+\frac{3\pi}{8g^{2}}\frac{1}{\left|\vec{n}_{2}\right|^{2}}\right)t}\,.$
(51)

In Fourier space, the dynamical scaling described by equation (1) takes the
form [73]

$\left\langle\tilde{\xi}_{\vec{n}}\left(0\right)\tilde{\xi}_{-\vec{n}}\left(t\right)\right\rangle=\left\langle\tilde{\xi}_{\vec{n}}\tilde{\xi}_{-\vec{n}}\right\rangle\Phi\left(n^{z}t\right)$
(52)

where $z$ denotes the dynamic exponent and $\Phi(u)$ denotes a scaling
function. Direct comparison with equation (51) immediately shows that for our
system, $z=4$ and $\Phi(u)$ exhibits an exponentially decay.

Indeed, for any linear system, the roughness exponent $\zeta$ and dynamic
exponent $z$ will be related by

$z=2\zeta+d\,,$ (53)

where $d$ denotes the dimension of the system [47, 48]. Accordingly, for our
case of $d=2$, it is unsurprising to find that the values $\zeta=1$ and $z=4$
satisfy this relationship. Importantly however, the fact that this linear
character extends beyond the large-scale behaviour described by the exponents
to subsequent corrections is nontrivial. Indeed, in the next section, we will
unexpectedly find that this linear character over all scales remains even in
the context of strong nonlinear coupling and is not merely an artefact of weak
coupling.

##  4 The Nonlinear Regime

While ordinary everyday systems are expected to reside in the linear regime of
$g\gg 1$, it is not difficult to imagine circumstances where $g\ll 1$ and thus
the nonlinear term will dominate equation (23). In particular, since equation
(24) shows that $g\sim h^{5/2}$, we expect $g$ to be small for sufficiently
thin materials such as graphene or biological membranes. Indeed, single layer
graphene will have the thickness of a single atom, i.e. $h\sim 1\textrm{
\AA}$, and thus a sheet of size $L\sim 250\textrm{ nm}$ [77] with Young’s
modulus $E\sim 2\textrm{ TPa}$ [78] and bending rigidity $\kappa\sim 2\times
10^{-19}\textrm{ N m}$ [79] at temperatures of say $T\sim 300\textrm{ K}$ will
have $g\sim 0.02$ which is indeed much smaller than unity. Similarly,
polymerised lipid bilayers surrounding cells or organelles ($L\sim 1\textrm{
$\mathrm{\SIUnitSymbolMicro m}$}-100\textrm{ $\mathrm{\SIUnitSymbolMicro
m}$}$) with thickness $h\sim 5\textrm{ nm}$ and Young’s modulus $E\sim
20\textrm{ MPa}$ [80] will have values of $g$ ranging from $\sim 0.003$ to
$\sim 0.3$. A theory which can also describe the small $g$ nonlinear regime is
thus certainly of value.

To investigate this regime, it will be helpful to rescale equation (23) as
follows. If we define a new dimensionless time $\bar{t}$ and vertical
displacement $\bar{\xi}_{\vec{n}}$ as

$\bar{t}=\frac{1}{g}\hat{t}=\left[\frac{\left(2\pi\right)^{4}2k_{\mathrm{B}}TEh}{\alpha^{2}L^{6}}\right]^{1/2}t$
(54)

and

$\bar{\xi}_{\vec{n}}(\bar{t}\,)=\frac{1}{g^{1/2}}\hat{\xi}_{\vec{n}}(\hat{t}\,)=\left[\frac{\left(2\pi\right)^{4}Eh}{2k_{\mathrm{B}}TL^{2}}\right]^{1/4}\tilde{\xi}_{\vec{n}}\left(t\right)\,,$
(55)

then equation (23) can be rewritten as

$\frac{\partial\left\langle\mathbb{F}\right\rangle}{\partial\bar{t}}=\frac{1}{2}\sum_{\vec{n}}\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_{\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\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\,,$
(56)

from which it is clear that $g$ is coupled to the linear part. Accordingly,
for $g\ll 1$ attempting to perturb the equation around the linear part is
singular and thus more sophisticated techniques are required if we wish to
understand this regime.

The Self-Consistent Expansion (SCE) is such a technique. Originally developed
to investigate the KPZ equation [81, 82, 83, 84, 85, 86, 87, 73, 88], the SCE
has been successfully used to study a variety of out-of-equilibrium systems
including turbulence, fracture and wetting fronts [89, 90, 91]. More recently
it was shown that the SCE is able to produce excellent results at low order
for the zero-dimensional $\phi^{4}$-theory and even exact convergence at high
orders [92, 93]. On this basis, the SCE was used in [41] and [42] to study the
related problem of a sheet fluctuating under conserved noise immediately prior
to crumpling and following this success, here too we will apply the SCE to
equation (56).

The idea behind the SCE is that we can expand equation (56) around any linear
theory, not just the one “handed to us”. That is, consider some arbitrary
linear coupling $\Gamma_{|\vec{n}|}$. We can identically insert this coupling
into equation (56) by writing

$\frac{\partial\left\langle\mathbb{F}\right\rangle}{\partial\bar{t}}=\frac{1}{2}\sum_{\vec{n}}\left\langle\frac{\partial^{2}\mathbb{F}}{\partial\bar{\xi}_{\vec{n}}\partial\bar{\xi}_{-\vec{n}}}\right\rangle-\sum_{\vec{n}}\Gamma_{|\vec{n}|}\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_{|\vec{n}|}\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_{\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\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\,.$
(57)

The goal is to choose the linear coupling $\Gamma_{|\vec{n}|}$ such that the
expansion around this linear theory is close to the behaviour of the nonlinear
system though this task will be deferred to later. For now, we can simply
perform an ordinary perturbative expansion in which the nonlinear term and
“residual” linear term are corrections, that is, we can write

$\frac{\partial\left\langle\mathbb{F}\right\rangle^{\left(m\right)}}{\partial\bar{t}}=\frac{1}{2}\sum_{\vec{n}}\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_{|\vec{n}|}\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_{|\vec{n}|}\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_{\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\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)}\,.$
(58)

In the following subsections, we will analyse the static and dynamic structure
factors which emerge from such a perturbative expansion.

###  4.1 The Static Structure Factor

As with the linear regime, the equal-time two-point and four-point functions
can be obtained by subbing in
$\mathbb{F}=\bar{\xi}_{\vec{n}_{1}}\bar{\xi}_{\vec{n}_{2}}$ and
$\mathbb{F}=\bar{\xi}_{\vec{n}_{1}}\bar{\xi}_{\vec{\ell}_{1}}\bar{\xi}_{\vec{\ell}_{2}}\bar{\xi}_{\vec{\ell}_{3}}$
and similarly to before, at zeroth order, we obtain

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

and Isserlis’ theorem

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

Proceeding to first order, subbing
$\mathbb{F}=\bar{\xi}_{\vec{n}_{1}}\bar{\xi}_{\vec{n}_{2}}$ and $m=1$ into
equation (58) gives

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

which after simplifying the sums by subbing in equation (60) and rearranging
gives

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

We now argue that the appropriate linear theory to expand around is the one
for which subsequent corrections to the linear theory are minimised, that is,
we self-consistently desire that

$\left\langle\bar{\xi}_{\vec{n}_{1}}\bar{\xi}_{\vec{n}_{2}}\right\rangle^{\left(1\right)}=\left\langle\bar{\xi}_{\vec{n}_{1}}\bar{\xi}_{\vec{n}_{2}}\right\rangle^{\left(0\right)}\,.$
(63)

Clearly this will be the case if $\Gamma_{|\vec{n}|}$ obeys the discrete
integral equation

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

This equation defines the linear coupling which we have desired to expand
around and though such equations are generically difficult to solve (indeed,
there may not exist an exact solution to this equation), we can relatively
easily obtain corresponding approximate asymptotic behaviours of
$\Gamma_{|\vec{n}|}$ by considering the limit of a continuous spectrum, though
this may result in missing solutions which also contain a dependence on the
direction of $\vec{n}$. The discrete frequencies $\vec{n}$ of our system will
appear continuous in the limit where the size of the sheet $L$ is taken to
infinity thus it is natural to define the scaled frequencies
$\vec{q}=\vec{n}/L$ and $\vec{k}=\vec{\ell}/L$. At the same time, it will be
helpful to scale $\Gamma_{|\vec{n}|}$ by some amount to keep it finite and so
we define $\Gamma(q)=\Gamma_{|\vec{n}|}/L^{\nu}$ where $\nu$ is an appropriate
scaling exponent which will become clear shortly. Then equation (64) becomes

$\Gamma\left(q\right)=L^{4-\nu}gq^{4}+\frac{1}{2}L^{6-2\nu}\sum_{\vec{k}\neq\vec{q}}\frac{1}{L^{2}}\frac{\left|\vec{q}\times\vec{k}\,\right|^{4}}{\Gamma\left(k\right)\left|\vec{q}-\vec{k}\,\right|^{4}}\,.$
(65)

Setting $\nu=3$ and taking the limit $L\rightarrow\infty$, we obtain the
continuous integral equation for $\Gamma(q)$

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

It is worth noting that the coefficient $Lg$ remains constant as
$L\rightarrow\infty$ since $g\sim 1/L$ and thus this continuous integral
equation is a well-defined continuous analogue of our discrete integral
equation. We can write the integral in polar coordinates as

$\Gamma\left(q\right)=(Lg)q^{4}\\\
+\frac{q^{4}}{2}\int_{0}^{q_{\mathrm{max}}}dk\,\frac{k^{5}}{\Gamma\left(k\right)}\int_{0}^{2\pi}d\theta\,\frac{\sin^{4}\theta}{\left(q^{2}+k^{2}-2qk\cos\theta\right)^{2}}\,,$
(67)

where $q_{\mathrm{max}}=n_{\mathrm{max}}/L$ denotes a potential upper
frequency cut-off. As mentioned towards the end of Section 3.1, if the system
has some finite resolution scale $a\ll L$, an upper frequency cut-off
$n_{\mathrm{max}}=L/(2a)\gg 1$ is introduced. It will become clear that in the
nonlinear regime, the presence of a finite upper frequency cut-off can result
in very different behaviour and thus we already allow for such a possibility
at this stage. The case with no upper cut-off can always be obtained later by
simply taking the limit $q_{\mathrm{max}}\rightarrow\infty$.

The angular integral has already been encountered and is given by equation
(34) thus we can obtain the following integral equation for $\Gamma(q)$

$\Gamma\left(q\right)=(Lg)q^{4}\\\
+\frac{3\pi}{8}\left[\int_{0}^{q}dk\,\frac{k^{5}}{\Gamma\left(k\right)}+q^{4}\int_{q}^{q_{\mathrm{max}}}dk\,\frac{k}{\Gamma\left(k\right)}\right]$
(68)

and with this equation in hand, we are ready to determine the asymptotic
behaviours of $\Gamma(q)$.

Beginning with large $q$, suppose that $\Gamma(q)$ has an asymptotic expansion

$\Gamma(q)\sim A_{\infty}q^{\nu_{\infty}}\,,$ (69)

where $A_{\infty}$ and $\nu_{\infty}$ determine the leading order behaviour of
$\Gamma(q)$ and higher order corrections are neglected. We can write equation
(68) as

$\Gamma\left(q\right)=(Lg)q^{4}+\frac{3\pi}{8}\int_{q}^{q_{\mathrm{max}}}dk\,\frac{k\left(q^{4}-k^{4}\right)}{\Gamma\left(k\right)}+C_{\infty}\,,$
(70)

where

$C_{\infty}=\frac{3\pi}{8}\int_{0}^{q_{\mathrm{max}}}dk\,\frac{k^{5}}{\Gamma\left(k\right)}$
(71)

denotes a constant. Now subbing equation (69) into equation (70), we obtain at
leading order that

$A_{\infty}q^{\nu_{\infty}}=\left[(Lg)-\frac{3\pi}{8A_{\infty}}\frac{q_{\mathrm{max}}^{2-\nu_{\infty}}}{\left(\nu_{\infty}-2\right)}\right]q^{4}+\frac{3\pi}{8A_{\infty}}\frac{4}{\left(6-\nu_{\infty}\right)\left(\nu_{\infty}-2\right)}q^{6-\nu_{\infty}}\\\
+C_{\infty}-\frac{3\pi}{8A_{\infty}^{\left(0\right)}}\frac{q_{\mathrm{max}}^{6-\nu_{\infty}}}{\left(6-\nu_{\infty}\right)}\,,$
(72)

where we have implicitly assumed that $\nu_{\infty}\neq 2,6$. (Assuming
$\nu_{\infty}=2,6$ results in logarithmic corrections on the right-hand side
which cannot be balanced and thus this assumption is justified.) Examination
of this equation reveals that there are two choices for $\nu_{\infty}$ which
can satisfy the equation at leading order, either $\nu_{\infty}=4$ or
$\nu_{\infty}=3$. If $\nu_{\infty}=4$, then we must also have

$A_{\infty}=(Lg)-\frac{3\pi}{8A_{\infty}}\frac{1}{2q_{\mathrm{max}}^{2}}\,,$
(73)

which can be solved for $A_{\infty}$ to give

$A_{\infty}=\frac{(Lg)}{2}\left[1+\sqrt{1-\frac{3\pi}{4}\frac{1}{(Lgq_{\mathrm{max}})^{2}}}\right]\,.$
(74)

We thus find that the first option for the large $q$ asymptotic behaviour of
$\Gamma(q)$ is

$\Gamma(q)\sim\frac{(Lg)}{2}\left[1+\sqrt{1-\frac{3\pi}{4}\frac{1}{(Lgq_{\mathrm{max}})^{2}}}\right]q^{4}\,.$
(75)

Recalling that $\Gamma(q)=\Gamma_{|\vec{n}|}/L^{3}$ and $q=n/L$, we find that
this asymptotic behaviour corresponds to

$\Gamma_{\left|\vec{n}\right|}\sim\frac{g}{2}\left[1+\sqrt{1-\frac{3\pi}{4}\frac{1}{(gn_{\mathrm{max}})^{2}}}\right]n^{4}\,,$
(76)

or when $gn_{\mathrm{max}}\gg 1$

$\Gamma_{\left|\vec{n}\right|}\sim
gn^{4}\left[1-\frac{3\pi}{16}\frac{1}{(gn_{\mathrm{max}})^{2}}\right]\,.$ (77)

Noting that

$\left\langle\hat{\xi}_{\vec{n}}\hat{\xi}_{-\vec{n}}\right\rangle^{\left(0\right)}=g\left\langle\bar{\xi}_{\vec{n}}\bar{\xi}_{-\vec{n}}\right\rangle^{\left(0\right)}=\frac{g}{2\Gamma_{|\vec{n}|}}\,,$
(78)

and comparing with equation (40), we thus find that at leading order, this
large frequency asymptotic behaviour completely agrees with the analysis
performed for the linear regime. It is important however to point out that
agreement merely requires that $gn_{\mathrm{max}}\gg 1$ rather than $g\gg 1$
and thus if $n_{\mathrm{max}}$ is sufficiently large, it is possible to obtain
a linear-like response for large frequencies, even if $g\ll 1$. If no upper
frequency cut-off exists whatsoever, this large frequency linear-like response
is actually obtained for any nonzero $g$.

At the same time, from equation (74), it is apparent that if

$gn_{\mathrm{max}}=Lgq_{\mathrm{max}}<\sqrt{\frac{3\pi}{4}}\approx 1.53...\,,$
(79)

then $A_{\infty}$ will cease to be real, i.e. $A_{\infty}\in\mathbb{C}$, which
cannot be the case. Accordingly, when $gn_{\mathrm{max}}$ is sufficiently
small, we find that the linear $n^{4}$ asymptotic behaviour becomes untenable.
Under these circumstances, the other possible asymptotic behaviour mentioned
above emerges, i.e. we must have $\nu_{\infty}=3$ in which case we simply find
that

$A_{\infty}=\sqrt{\frac{\pi}{2}}\,.$ (80)

This eventuality however can only be realised if the coefficient of the
$q^{4}$ term in equation (72) vanishes, that is we must also have

$Lgq_{\mathrm{max}}=\frac{3\pi}{8A_{\infty}}\,,$ (81)

i.e.

$gn_{\mathrm{max}}\approx\frac{3\sqrt{2\pi}}{8}\approx 0.94...\,.$ (82)

We thus find that for values of $gn_{\mathrm{max}}$ close to or smaller than
1, the linear-like quartic large frequency behaviour ceases to exist and a
cubic behaviour becomes dominant.

Moving onto the small $q$ behaviour, suppose that $\Gamma(q)$ has a small $q$
expansion

$\Gamma(q)\sim A_{0}q^{\nu_{0}}\,,$ (83)

where $A_{\infty}$ and $\nu_{0}$ determine the small $q$ leading order
behaviour. We can write equation (68) as

$\Gamma\left(q\right)=\left[(Lg)+C_{0}\right]q^{4}\\\
+\frac{3\pi}{8}\left[\int_{0}^{q}dk\,\frac{k^{5}}{\Gamma\left(k\right)}-q^{4}\int_{\lambda}^{q}dk\,\frac{k}{\Gamma\left(k\right)}\right]\,,$
(84)

where

$C_{0}=\frac{3\pi}{8}\int_{\lambda}^{q_{\mathrm{max}}}dk\,\frac{k}{\Gamma\left(k\right)}$
(85)

and $\lambda<q$ denotes some small frequency which can be taken to 0 but is
temporarily needed to prevent the integrals diverging. Subbing equation (83)
into equation (84) and expanding order by order, we obtain

$A_{0}q^{\nu_{0}}=\left[(Lg)+C_{0}-\frac{3\pi}{8A_{0}}\frac{\lambda^{2-\nu_{0}}}{\nu_{0}-2}\right]q^{4}\\\
+\frac{3\pi}{8A_{0}}\frac{4}{\left(6-\nu_{0}\right)\left(\nu_{0}-2\right)}q^{6-\nu_{0}}\,,$
(86)

which can only be satisfied at lowest order by $\nu_{0}=3$, in which case

$A_{0}=\sqrt{\frac{\pi}{2}}\,.$ (87)

This leading order cubic small frequency behaviour is obtained independent of
the value of $g$ and it is curious to observe that at the critical value
$gn_{\mathrm{max}}=3\sqrt{2\pi}/8$ that $A_{0}=A_{\infty}$. Indeed
substituting

$\Gamma(q)=\sqrt{\frac{\pi}{2}}q^{3}$ (88)

directly into equation (68) makes it clear that this is the exact solution to
the integral equation when
$Lgq_{\mathrm{max}}=gn_{\mathrm{max}}=3\sqrt{2\pi}/8$.

To summarise, we have identified the following three asymptotic behaviours of
$\Gamma(q)$

  * •

For large $q$, when $Lgq_{\mathrm{max}}\gg 1$, we have found that

$\Gamma(q)\sim\frac{(Lg)}{2}q^{4}\,.$ (89)

This behaviour perfectly corresponds to the behaviour we found in the linear
regime.

  * •

For large $q$, when $Lgq_{\mathrm{max}}\lesssim 1$, we have found that

$\Gamma(q)\sim\sqrt{\frac{\pi}{2}}q^{3}\,.$ (90)

  * •

For small $q$, regardless of the value of $g$, we have found that

$\Gamma(q)\sim\sqrt{\frac{\pi}{2}}q^{3}\,.$ (91)

Together, these asymptotics suggest that the overall behaviour can be
approximated by the combination

$\Gamma(q)\approx
C\left(Lg,q_{\mathrm{max}}\right)q^{3}+A\left(Lg,q_{\mathrm{max}}\right)q^{4}$
(92)

where $A(Lg,q_{\mathrm{max}})$ and $C\left(Lg,q_{\mathrm{max}}\right)$ denote
constants which depend on $Lg$ and $q_{\mathrm{max}}$, though
$C\sim\sqrt{\pi/2}\sim 1$ is expected to have a much weaker dependence on $g$
and thus will generally be treated as constant. In terms of the discrete
frequencies we began with, this translates to the functional form

$\Gamma_{\left|\vec{n}\right|}\approx
Cn^{3}+A\left(g,n_{\mathrm{max}}\right)n^{4}\,.$ (93)

We now note that the crossover between the $n^{3}$ and $n^{4}$ behaviours
occurs when
$\left|C\right|n^{3}\sim\left|A\left(g,n_{\mathrm{max}}\right)\right|n^{4}$
and thus we can define a crossover frequency $n_{c}$ as

$n_{c}\sim\left|\frac{C}{A\left(g,n_{\mathrm{max}}\right)}\right|\,.$ (94)

At the same time, $n_{c}$ must obey $1\leq n_{c}\leq n_{\mathrm{max}}$ thus we
expect that if $n_{c}$ is too small, a purely quartic behaviour will be
obtained (as was predicted for the linear regime) while if $n_{c}$ is too
large, a purely cubic behaviour will be obtained. For values of $g$ and
$n_{\mathrm{max}}$ which cause $n_{c}$ to fall between these two extremes, a
state of coexistence is obtained in which the cubic behaviour dominates the
small frequencies while the quartic behaviour dominates the large frequencies.

We finish this analysis by providing a useful heuristic for
$A(Lg,q_{\mathrm{max}})$. Differentiating equation (68) with respect to $q$,
one obtains

$\Gamma^{\prime}\left(q\right)=4(Lg)q^{3}+\frac{3\pi}{2}q^{3}\int_{q}^{q_{\mathrm{max}}}dk\,\frac{k}{\Gamma\left(k\right)}\,.$
(95)

Subbing in $q=q_{\mathrm{max}}$ thus gives

$\Gamma^{\prime}\left(q_{\mathrm{max}}\right)=4(Lg)q_{\mathrm{max}}^{3}\,.$
(96)

At the same time, differentiating equation (92) and subbing in
$q=q_{\mathrm{max}}$ gives

$\Gamma^{\prime}(q_{\mathrm{max}})\approx
3Cq_{\mathrm{max}}^{2}+4A\left(Lg,q_{\mathrm{max}}\right)q_{\mathrm{max}}^{3}\,.$
(97)

Accordingly, after equating, we obtain

$A\left(Lg,q_{\mathrm{max}}\right)\approx Lg-\frac{3C}{4q_{\mathrm{max}}}\,,$
(98)

or in terms of discrete frequencies

$A\left(g,n_{\mathrm{max}}\right)\approx g-\frac{3C}{4n_{\mathrm{max}}}\,.$
(99)

$A\left(g,n_{\mathrm{max}}\right)$ is thus found to grow linearly with $g$ and
crucially, for small $g$, it is found to vanish and even becomes negative. The
crossover frequency $n_{c}$ defined above in equation (94) occurs around

$n_{c}\sim\left|\frac{C}{gn_{\mathrm{max}}-\frac{3}{4}C}\right|n_{\mathrm{max}}\,,$
(100)

thus if we take $C\sim 1$, we expect a heavily dominated quartic behaviour if
$g\gtrsim 1$ and a heavily dominated cubic behaviour if $g\lesssim
2/n_{\mathrm{max}}$. For values between these extremes, neither the cubic
contribution nor the quartic contribution can be neglected.

To summarise, we have found that, up to first order, the static structure
factor in the nonlinear regime behaves like

$\left\langle\bar{\xi}_{\vec{n}}\bar{\xi}_{-\vec{n}}\right\rangle^{\left(1\right)}=\frac{1}{2\Gamma_{|\vec{n}|}}=\frac{1}{2n^{3}\left[C+A\left(g,n_{\mathrm{max}}\right)n\right]}\,,$
(101)

where $C\sim 1$ and $A\left(g,n_{\mathrm{max}}\right)$ is approximately given
by equation (99). In dimensional terms, this reads as

$\left\langle\tilde{\xi}_{\vec{n}}\tilde{\xi}_{-\vec{n}}\right\rangle^{\left(1\right)}=\frac{\sqrt{\left(2k_{\mathrm{B}}TL^{2}\right)/\left[\left(2\pi\right)^{4}Eh\right]}}{2n^{3}\left[C+A\left(g,n_{\mathrm{max}}\right)n\right]}\,.$
(102)

Note that the bending rigidity $\kappa$ does not enter the expression
explicitly except via the dependence of $A\left(g,n_{\mathrm{max}}\right)$ on
$g$. Applying the same logic as at the end of section 3.1, we find upon
comparison with equation (39) that the roughness exponent of the sheet in the
nonlinear regime is $\zeta=1/2$ which is found to be in agreement with results
presented in [7]. It is important however to emphasise that this roughness
exponent only characterises the sheet well if we are in the extremely strong
nonlinear regime of $g\lesssim 2/n_{\mathrm{max}}$. When one finds oneself in
the interval with $2/n_{\mathrm{max}}\lesssim g\lesssim 1$, the existence of
the “linear” quartic crossover cannot be neglected and thus the roughness
exponent alone is less informative. Indeed, since
$n_{\mathrm{max}}=L/\left(2a\right)$ where $a$ is the lattice length of the
material, the extremely strong nonlinear regime of $g\lesssim
2/n_{\mathrm{max}}$ will only be obtained if

$\frac{\left(2\pi\right)^{2}\kappa}{\sqrt{2k_{\mathrm{B}}TEh\left(4a\right)^{2}}}\lesssim
1\,.$ (103)

For graphene, for instance, with a measured and theoretical bond-length of
$a\sim 0.15\textrm{ nm}$ [79], the left-hand side of this expression is
actually roughly 8 and thus even graphene does not fall into this most extreme
regime.

###  4.2 The Dynamic Structure Factor

As in the linear regime, we can obtain time-dependent two-point and four point
functions by subbing
$\mathbb{F}=\bar{\xi}_{\vec{n}_{1}}\left(0\right)\bar{\xi}_{\vec{n}_{2}}(\bar{t}\,)$
and
$\mathbb{F}=\bar{\xi}_{\vec{n}_{1}}\left(0\right)\bar{\xi}_{\vec{\ell}_{1}}\bar{\xi}_{\vec{\ell}_{2}}\bar{\xi}_{\vec{\ell}_{3}}$
into equation (58). At zeroth order ($m=0$), in complete analogy with the
linear regime, this gives ODEs with solutions

$\left\langle\bar{\xi}_{\vec{n}_{1}}\left(0\right)\bar{\xi}_{\vec{n}_{2}}(\bar{t}\,)\right\rangle^{\left(0\right)}=\left\langle\bar{\xi}_{\vec{n}_{1}}\bar{\xi}_{\vec{n}_{2}}\right\rangle^{\left(0\right)}e^{-\Gamma_{|\vec{n}_{2}|}\bar{t}}$
(104)

and

$\left\langle\bar{\xi}_{\vec{n}_{1}}\left(0\right)\bar{\xi}_{\vec{\ell}_{1}}\bar{\xi}_{\vec{\ell}_{2}}\bar{\xi}_{\vec{\ell}_{3}}\right\rangle^{\left(0\right)}=\left\langle\bar{\xi}_{\vec{n}_{1}}\bar{\xi}_{\vec{\ell}_{1}}\bar{\xi}_{\vec{\ell}_{2}}\bar{\xi}_{\vec{\ell}_{3}}\right\rangle^{\left(0\right)}e^{-\Gamma_{|\vec{n}_{1}|}\bar{t}}\,,$
(105)

where the zeroth order equal-time two-point and four-point functions are given
by equations (59) and (60).

We can obtain the first order correction to the dynamic two-point function by
subbing $m=1$ and
$\mathbb{F}=\bar{\xi}_{\vec{n}_{1}}\left(0\right)\bar{\xi}_{\vec{n}_{2}}(\bar{t}\,)$
into equation (58). After simplification, we obtain the ODE

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

As pointed out previously, since the nonhomogeneous part of this equation
decays at the natural decay rate of the equation, the solution to this
equation will contain a nonphysical secular term unless the coefficient of the
nonhomogeneous part happens to vanish. This will occur if

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

but this is none other than the discrete integral equation for
$\Gamma_{|\vec{n}|}$ obtained in the previous subsection (equation (64))!
Accordingly, all the results obtained for $\Gamma_{|\vec{n}|}$ in the previous
subsection apply here as well and we thus find that the dynamic structure
factor is given up to first order by an exponential decay with decay rate
$\Gamma_{|\vec{n}|}$

$\left\langle\bar{\xi}_{\vec{n}}\left(0\right)\bar{\xi}_{-\vec{n}}(\bar{t}\,)\right\rangle^{\left(1\right)}=\left\langle\bar{\xi}_{\vec{n}}\bar{\xi}_{-\vec{n}}\right\rangle^{\left(1\right)}e^{-\Gamma_{|\vec{n}|}\bar{t}}\,.$
(108)

The above argument that secular terms should be eliminated has already been
alluded to in Section 3.2 when deriving the dynamic structure factor in the
linear regime. There it was observed that the ordinary application of
perturbation theory fails to eliminate the subsequent secular term and thus
more sophisticated techniques such as the Poincaré-Lindstedt method or
resummation via Padé approximants must be used. Here in contrast, the SCE
immediately provides a means of eliminating the secular term. Note that
eliminating the secular term is conceptually equivalent to imposing that the
zeroth order approximation be exact up to first order

$\left\langle\bar{\xi}_{\vec{n}}\left(0\right)\bar{\xi}_{-\vec{n}}(\bar{t}\,)\right\rangle^{\left(1\right)}=\left\langle\bar{\xi}_{\vec{n}}\left(0\right)\bar{\xi}_{-\vec{n}}(\bar{t}\,)\right\rangle^{\left(0\right)}$
(109)

and thus can also be construed as a self-consistent argument for
$\Gamma_{|\vec{n}|}$.

The fact that the exact same integral equation is obtained for
$\Gamma_{|\vec{n}|}$ from the analysis of the static and dynamic structure
factors was by no means guaranteed and indeed, it is known to not be the case
for the ordinary $\phi^{4}$-theory. Indeed, this fact suggests that the system
is in some sense quasi-linear, i.e. it behaves like a purely linear theory
with linear coupling given by $\Gamma_{|\vec{n}|}$. In dimensional terms,
equation (108) can be written as

$\left\langle\tilde{\xi}_{\vec{n}}\left(0\right)\tilde{\xi}_{-\vec{n}}(t)\right\rangle^{\left(1\right)}=\\\
\left\langle\tilde{\xi}_{\vec{n}}\tilde{\xi}_{-\vec{n}}\right\rangle^{\left(1\right)}e^{-\frac{\left(2\pi\right)^{2}\sqrt{2k_{\mathrm{B}}TEh}}{\alpha
L^{3}}n^{3}\left[C+A\left(g,n_{\mathrm{max}}\right)n\right]t}\,,$ (110)

which makes it clear that at large scales, the dynamic exponent is simply
$z=3$ which, together with our identification of $\zeta=1/2$ for the nonlinear
regime, indeed satisfies the relationship expected for a linear theory given
by equation (53). It is important to note however that this quasi-linearity
extends well beyond merely the large scale behaviour of the sheet
characterised by these exponents. The fact that $\Gamma_{|\vec{n}|}$ has the
same functional form for both the static and dynamic structures at all
frequency scales in fact suggests that this quasi-linearity occurs at all
scales.

##  5 Comparison with Simulations

Equation (13) is a straightforward nonlinear Langevin equation and can thus be
simulated using standard numerical techniques. First though, it will be
helpful to nondimensionalise the equation. We do so using the
nondimensionalisation used for the nonlinear analysis performed in Section 4,
that is, using equations (54) and (55) while scaling the noise by

$\bar{\eta}_{\vec{n}}\left(\bar{t}\right)=\left[\frac{L^{10}}{\left(2\pi\right)^{4}Eh\left(2k_{\mathrm{B}}T\right)^{3}}\right]^{1/4}\tilde{\eta}_{\vec{n}}\left(t\right)\,,$
(111)

we obtain the following dimensionless Langevin 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)\,,$
(112)

where the variance of the dimensionless noise is simply

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

Following [41] and [42], for various values of $g$, we simulate equation (112)
by forward iteration over a $41\times 41$ square lattice implying a maximum
frequency $n_{\mathrm{max}}=20$. We do this for both large and small values of
$g$ corresponding to the linear and nonlinear regimes respectively. To
determine the dimensionless time step $\delta\bar{t}$ and time duration
$\bar{T}$ of each simulation, we recall that equations (93) and (99) suggest
that modes of size $n$ decay at rate

$\Gamma_{n}\sim Cn^{3}+\left(g-\frac{3C}{4n_{\mathrm{max}}}\right)n^{4}\,,$
(114)

where $C\sim 1$ is relatively constant across values of $g$. Subbing in the
minimum frequency $n=1$ and maximum frequency $n=n_{\mathrm{max}}$, we obtain
two time scales

$\tau_{\mathrm{min}}\sim\frac{1}{\Gamma_{n_{\mathrm{max}}}}\sim\frac{4}{n_{\mathrm{max}}^{3}\left(1+4gn_{\mathrm{max}}\right)}$
(115)

and

$\tau_{\mathrm{max}}\sim\frac{1}{\Gamma_{1}}\sim\frac{1}{1+g}\,.$ (116)

The fastest decaying modes decay with characteristic time
$\tau_{\mathrm{min}}$ thus our time step $\delta\bar{t}$ must be smaller than
this if we wish to faithfully simulate the equation. Accordingly, for each
simulation, we select the time step $\delta\bar{t}$ to be equal to the largest
power of $2$ smaller than $\tau_{\mathrm{min}}$ (i.e.
$\delta\bar{t}=2^{\max(k)}<\tau_{\mathrm{min}}$). It is worth noting that in
practice, if the time step $\delta\bar{t}$ is larger than
$\tau_{\mathrm{min}}$, the simulation becomes unstable and thus breaks down
after a small number of steps. Accordingly, on a practical level, we can use
the stability of the simulation to gauge whether our time step is sufficiently
small.

The slowest decaying modes decay with characteristic time
$\tau_{\mathrm{max}}$ thus if we wish to observe their decay, we must run the
simulations over time intervals $\bar{T}$ larger than $\tau_{\mathrm{max}}$.
Accordingly, we select the time duration of each simulation to be equal to $4$
times the smallest power of $2$ larger than $\tau_{\mathrm{max}}$ (i.e.
$\bar{T}/4=2^{\min(k)}>\tau_{\mathrm{max}}$), the factor of $4$ being chosen
to ensure that we have at least several independent intervals in which the
decay is guaranteed to occur.

The aforementioned decision to use a $41\times 41$ lattice, i.e.
$n_{\mathrm{max}}=20$, accounts for the competing interests of spatial
resolution and temporal resolution in our simulation. Equation (115) makes it
clear that finer spatial resolution (i.e. larger values of $n_{\mathrm{max}}$)
requires substantially finer temporal resolution. Since we wish to be able to
run the simulations for long enough to obtain statistics for even the slowest
decaying modes, using a temporal resolution which is extremely fine is
computationally expensive. Running simulations with different values of
$n_{\mathrm{max}}$, it was observed that $n_{\mathrm{max}}=20$ handles this
trade-off satisfactorily, that is, the spatial resolution is sufficiently fine
for our theory to be applicable while the temporal resolution is not too fine
to make the simulations too onerous to run.

Finally, since it is time consuming to calculate the nonlinear quartic term of
equation (112) directly, we implement the pseudo-spectral method described in
[57] and used in [41] and [42] in which the term is calculated as the Fourier
transform of its real-space counterpart. This has the additional effect of
imposing periodic boundary conditions on the simulation. For each value of $g$
used, 10 simulations were run and the results averaged across simulations.

Figure 1: The static structure factor
$\left\langle\bar{\xi}_{\vec{n}}\bar{\xi}_{-\vec{n}}\right\rangle$ as a
function of $n$ for various values of $g$. The data points (symbols) denote
simulation data while the solid lines denote fits of equation (117) to this
data. The dashed and dotted lines are guidelines proportional to $n^{-3}$ and
$n^{-4}$ respectively. For large values of $g$, the $n^{-4}$ behaviour
dominates, confirming the prediction of the linear theory while for small
values of $g$, the $n^{-3}$ behaviour dominates, confirming the prediction of
the nonlinear theory. For intermediate values, a crossover between $n^{-3}$
behaviour for small $n$ and $n^{-4}$ behaviour for large $n$ is observed.

Figure 1 shows the static structure factor for several values of $g$ ranging
from significantly smaller than unity to significantly larger than unity. The
points (symbols) denote the simulation data while the solid lines denote fits
to the data of the equation

$\left\langle\bar{\xi}_{\vec{n}}\bar{\xi}_{-\vec{n}}\right\rangle=\frac{1}{2\left(Cn^{3}+An^{4}\right)}\,,$
(117)

where $A$ and $C$ are treated as fitting parameters. The dashed and dotted
lines are guidelines proportional to $n^{-3}$ and $n^{-4}$ respectively. As
can be seen, as $g$ increases, the simulation data moves from being roughly
parallel to the $n^{-3}$ guideline, predicted for the theory of the nonlinear
regime, to being roughly parallel to the $n^{-4}$ guideline, predicted for the
theory of the linear regime. For intermediate values of $g$, a clear crossover
between $n^{-3}$ behaviour for small $n$ and $n^{-4}$ behaviour for large $n$
is observed. The fact that the fits are excellent across all values of $n$ for
all values of $g$ confirms that equation (117) provides the correct functional
form of the static structure factor. The fit parameters $C(g)$ and $A(g)$ are
plotted in figure 5 and will be discussed in the continuation.

Before moving to the dynamic structure factor, it is worth observing that it
is easy to naively fit power-laws $n^{-\nu}$ to the data presented in figure 1
and this will give values of $\nu$ between 3 and 4. Indeed, because the change
from $n^{-3}$ behaviour to $n^{-4}$ behaviour is visually quite subtle, this
will even be the case when one only attempts to fit the small $n$ values.
Previous research has indeed observed a variety of anomalous scaling exponents
between 3 and 4 [3, 4, 18, 8, 31, 10, 2, 32, 33, 34, 35, 40, 38, 15, 21] and
even observed a temperature dependence [37] despite claims of universality. A
similar state of affairs was previously observed in the context of fracture
where a range of exponents could be extracted from naive power-law fitting of
data that was known to be exactly described by a rational function, that is, a
ratio of two polynomials [94, 95]. Here too, we suggest that the range of
various anomalous exponents is merely an artefact of the effect of the large
$n$ linear behaviour on the small $n$ leading nonlinear character.

Figure 2: The time-dependent structure factor
$\left\langle\bar{\xi}_{\vec{n}}(0)\bar{\xi}_{-\vec{n}}(\bar{t})\right\rangle$
as a function of $\bar{t}$ for $n=1,2,3,4,5$ when $g=0.5$. The solid lines
denote fits of exponential decays to each data set. Since the data appears to
fall on straight-lines on the logarithmic y-axis, we see that the time-
dependent two-point function indeed decays exponentially, as predicted by the
theory.  Figure 3: The exponential decay rate
$\Gamma_{|\vec{n}|}^{\mathrm{dyn}}$ of the time-dependent structure factor as
a function of $n$ for various values of $g$. The data points denote the
simulation data while the solid lines denote fits of equation (118) to the
data. The dashed and dotted lines are guidelines proportional to $n^{3}$ and
$n^{4}$ respectively. For large values of $g$, the $n^{4}$ behaviour
dominates, confirming the prediction of the linear theory while for small
values of $g$, the $n^{3}$ behaviour dominates, confirming the prediction of
the nonlinear theory. For intermediate values, a crossover between $n^{3}$
behaviour for small $n$ and $n^{4}$ behaviour for large $n$ is observed.

To confirm our predictions for the dynamical properties of the system, the
time-dependent structure factor can be plotted as a function of $\bar{t}$ and
fitted by an exponential decay from which the decay rate
$\Gamma_{|\vec{n}|}^{\mathrm{dyn}}$ can be extracted as suggested by equation
(108). Figure 2 shows this decay for several values of $n$ for $g=0.5$ (safely
in the nonlinear regime) with the solid lines denoting fits of exponential
decays to each data set. Since the y-scale is logarithmic, the decays appear
as straight lines confirming our theory that indeed, the time-dependent
structure factor decays exponentially with time. Similar results can be
obtained for other values of $n$ and $g$ and the results are presented in
figure 3 (symbols) together with fits (solid lines) of the equation

$\Gamma_{\left|\vec{n}\right|}^{\mathrm{dyn}}=Cn^{3}+An^{4}\,.$ (118)

The dashed and dotted lines are guidelines proportional to $n^{3}$ and $n^{4}$
respectively and it can be seen that as $g$ increases, the behaviour of the
data shifts from being parallel to the $n^{3}$ guideline, predicted for the
nonlinear theory, to being parallel to the $n^{4}$ guideline, predicted by the
linear theory. Indeed, the fact that the fits of equation (118) are excellent
over all values of $n$ for all values of $g$ confirms that equation (118)
correctly predicts the functional form of the decay rate of the time-dependent
structure factor. The fit parameters $C(g)$ and $A(g)$ are plotted in figure 5
and will be discussed shortly.

Figure 4: The ratio of the decay rate of the time-dependent structure factor,
$\Gamma_{|\vec{n}|}^{\mathrm{dyn}}$, to the effective linear coupling,
$\Gamma_{|\vec{n}|}^{\mathrm{stat}}$, given by equation (119) as a function
$n$ for various $g$. The ratio can be seen to be roughly equal to 1 for all
values of $g$ except when $n$ becomes small at which point the ratio increases
slightly.

In addition to predicting the functional forms of the static structure factor
(equation (117)) and the decay rate of the time-dependent structure factor
(equation (118)), our theory predicts that the system is quasilinear at all
scales, that is to say, that the effective linear coupling
$\Gamma_{\left|\vec{n}\right|}^{\mathrm{stat}}$, which can be defined using

$\left\langle\bar{\xi}_{\vec{n}}\bar{\xi}_{-\vec{n}}\right\rangle=\frac{1}{2\Gamma_{\left|\vec{n}\right|}^{\mathrm{stat}}}\,,$
(119)

should equal the decay rate of the time-dependent structure factor
$\Gamma_{\left|\vec{n}\right|}^{\mathrm{dyn}}$, i.e. we should have

$\frac{\Gamma_{\left|\vec{n}\right|}^{\mathrm{dyn}}}{\Gamma_{\left|\vec{n}\right|}^{\mathrm{stat}}}=1$
(120)

for all $\vec{n}$ and $g$. Figure 4 shows this ratio for various values of $g$
as a function of $n$ and indeed, for all $g$, this ratio is close to 1 except
when $\vec{n}$ is small at which point the ratio becomes slightly yet
noticeably larger than 1. The fact this ratio is close to 1 for most $n$ shows
that our prediction that the system is quasi-linear is indeed mostly correct.
For large $g$, we expect the system to become increasingly linear and indeed,
the deviation for small $n$ appears to shrink with increasing $g$. The
deviation itself is surprising considering the fact that figure 1 and figure 3
each individually show excellent agreement across all values of $n$. We are
thus forced to conclude that though equations (117) and (118) are each
functionally correct, the fit parameters $C(g)$ and $A(g)$ are not necessarily
identical.

Figure 5: The fit parameters, (a) $C(g)$, (b) $A(g)$, obtained by fitting
equations (117) and (118) to the static data presented in figure 1 (blue +)
and the dynamic data presented in figure 3 (red $\circ$) respectively. In (a),
the solid and dashed lines are the average value of the static and dynamic
$C(g)$ respectively. In (b), the solid and dashed lines are linear fits to the
static and dynamic $A(g)$ respectively.

Figure 5 shows the fit parameters $C(g)$ and $A(g)$ as functions of $g$ with
the static values (blue +) coming from fits of equation (117) to the data
displayed in figure 1 and the dynamic values (red $\circ$) coming from fits of
equation (118) to the data displayed in figure 3. As predicted by the theory,
the values of $C(g)$ remain relatively constant while the values of $A(g)$
grow with increasing $g$.

Since the theory predicts that $C(g)$ is relatively constant, the solid and
dashed lines in figure 5(a) denote the average value of $C(g)$ obtained from
the static and dynamic data respectively and, as can be seen, the deviations
around this average value are indeed small over the entire range of $g$ except
when $g$ becomes large. Our theory predicts that $C(g)\sim 1$ and this too is
indeed the case for most values of $g$ for both sets of $C(g)$. Since our
theory predicts that for large $g$, $\Gamma_{|\vec{n}|}$ will grow mostly like
$n^{4}$, the sudden growth in $C(g)$ when $g$ is large is merely an artefact
of the result that $C(g)$ isn’t needed as a fit parameter in this regime.

According to equation (99), $A(g)$ is expected to grow linearly. The solid and
dashed lines in figure 5(b) denote linear fits to the $A(g)$ obtained from the
static and dynamic data respectively and can be seen to be excellent over all
values of $g$. It is worth noting that according to equation (99), $A(g)$ is
expected to become negative for sufficiently small values of $g$ and this
behaviour is clearly visible on the log-log plot for the dynamic $A(g)$ where
the dashed line suddenly plummets to $-\infty$. Indeed, the static $A(g)$ and
dynamic $A(g)$ are almost identical for large values of $g$ and only really
begin to differ for small values of $g$ with the dynamic $A(g)$ becoming
negative while the static $A(g)$ remains positive.

Accordingly, it seems that the small $n$ deviation visible in figure 4 is
primarily due to the small discrepancy between the values of $A(g)$ obtained
from the static and dynamic data. This discrepancy in the coefficient of the
quartic behaviour, indicating that the quasi-linearity is not perfect, is
completely unpredicted by our theory though it is worth emphasising that it
only appears to become significant for small values of $g$. Conceptually, we
seem to have found that the quartic behaviour of the decay rate of the time-
dependent two-point function is a renormalisation of the quartic behaviour of
the effective linear coupling. Properly accounting for this renormalisation is
left as an open question.

##  6 Discussion

We have used the formalism of statistical mechanics to develop a Langevin
equation describing an elastic thin sheet being driven by white or thermal
noise. By investigating this so-called dynamic Föppl-von Kármán equation, we
have found that the behaviour of such a sheet can be characterised as
belonging to one of three “classes” depending on the value of the single
dimensionless parameter $g$, given by equation (20). For $g\gg 1$, the system
resides in the mostly linear regime where the sheet’s curvature plays little
role and thus the equal-time structure factor in Fourier space
$\langle\tilde{\xi}_{\vec{n}}\tilde{\xi}_{-\vec{n}}\rangle$ is found to decay
like $\sim n^{-4}$ with corresponding roughness exponent $\zeta=1$. For $g\ll
2/n_{\textrm{max}}$, the system resides in the highly nonlinear regime where
the effects of curvature dominate at all scales and thus, consistent with
previously analyses [7, 19], the equal-time structure factor is renormalised
to decay like $\sim n^{-3}$ with corresponding roughness exponent $\zeta=1/2$.
In between these two extremes, an intermediate behaviour is obtained where the
equal-time structure factor is found to decay like $\sim n^{-3}$ for small
frequencies (large scales) and $\sim n^{-4}$ for large frequencies (small
scales). This intermediate scale is of particular interest as it seems likely
that many ultra-thin, novel two-dimensional materials such as graphene lie
within it.

Indeed, there is much research indicating that for thin and ultra-thin
materials, the equal-time structure factor decays with an exponent between 3
and 4 [3, 4, 18, 8, 31, 10, 2, 32, 33, 34, 35, 40, 38, 15, 21], with
corresponding roughness exponent $\zeta$ between 1/2 and 1, or is even
temperature dependent [37]. Here, we argue that many of these results can be
better explained as insufficiently accounting for the crossover to the linear
tail. The transition from an $n^{-3}$ behaviour to an $n^{-4}$ behaviour can
be quite subtle and misidentifying the location of the crossover by even a
small amount can cause naive power-law fits to produce anomalous exponents. A
similar phenomenon is reported in the context of fracture in references [94,
95]. The confirmation of our results provided by our simulations supports this
position.

In addition to recovering and extending the predicted theory for the static
structure factor, we have used our model to provide the first concrete
prediction for the dynamic structure factor
$\langle\tilde{\xi}_{\vec{n}}\left(0\right)\tilde{\xi}_{-\vec{n}}\left(t\right)\rangle$,
finding analytic predictions in all regimes for the exponential decay towards
equilibrium and predicting a dynamic exponent $z$ given by $z=2\zeta+2$.
Importantly, this relationship between the dynamic exponent $z$ and the
roughness exponent $\zeta$ holds in all regimes and thus even in the most
extreme nonlinear regime, the model we have developed can be thought of as
quasi-linear.

This quasi-linearity can to some extent be interpreted as a consequence of a
pair of known exponent inequalities [48, 47]. In particular, it is known that
for systems whose dynamics can be derived from a Hamiltonian, that the dynamic
exponent $z$ must obey the inequality $z\geq 2\zeta+d$ where $d$ denotes the
dimension of the system and is thus equal to 2 in our case. Simultaneously, it
is known that Galilean invariant systems, i.e. systems for which the dynamics
of higher modes do not depend on the behaviour of the zeroth mode
$\tilde{\xi}_{0}$, must obey the reverse inequality $z\leq 2\zeta+d$. Since
our system can indeed be derived from a Hamiltonian [61, 2, 57, 19] and has
been observed to be manifestly Galilean invariant [see the discussion
immediately following equation (16)], these two inequalities collapse into the
single equality $z=2\zeta+d$, independent of the size of the nonlinear
coupling. From this perspective, it is an important and nontrivial example of
a physical system where both inequalities apply. It is worth noting however
that the application of these inequalities is limited in scope to only the
largest scales while we have found in our particular system that the quasi-
linearity applies at all scales. That is, the exponent inequalities mentioned
in [48, 47] only guarantee that the leading order powers $\zeta$ and $z$ will
be related by $z=2\zeta+d$ while we have found that the prefactors and
subleading orders also coincide. Indeed, our simulations show that the only
significant deviation from quasi-linear behaviour is in the large frequency,
small scale coefficient $A(g)$ and thus the exponent inequalities are not
contradicted. Accordingly, this quasi-linearity across all scales regardless
of the size of the nonlinear coupling is both unexpected and highly
significant. To the best of our knowledge, a theory which could a priori
predict this broad degree of quasi-linearity directly from the structure of
the equations is lacking despite its obvious great value.

It is also worth mentioning that in fact, not only do the higher modes not
depend on the zeroth mode (Galilean invariance) but the zeroth mode also does
not depend on the higher modes. Accordingly, the sheet’s center of mass
performs a random walk in space without being affected by the structure of the
sheet, while also not affecting its structure. This prediction is amenable to
experimental test and could help decipher the dynamics of the sheet. In
particular, if no such decoupling is observed, or if the stochastic motion of
the center of mass exhibits a different pattern, this could be valuable
information.

The results discussed here could also have an interesting application to the
optical response of thin sheets, or actually any waves interacting with the
structure, including acoustic waves. In particular, crumpled papers are known
to exhibit speckle patterns [96] and to emit a well-distinguishable crackling
noise [97, 98]. Here too, it would be interesting to see if optical effects
could be observed in graphene, which might be of practical use.

Finally, while the focus of this work has been on understanding the structure
of thin sheets by studying the impact of thermal fluctuations on their
curvature, there are many other curvature mediated means of effecting the
structure of thin matter. For instance, it is known that the surface
composition of thin films also modifies their bending moduli in a manner not
dissimilar to that of thermal fluctuations [99]. The equilibrium composition
of polymer-blend films, however, remains only partially understood (see for
instance [100]). Accordingly, it would be interesting to ask whether the
theory and methods developed here could provide insight into other types of
curvature-structure interactions.

The authors wish to thank Arezki Boudaoud for useful discussions. This work
was supported by the Israel Science Foundation Grant No. 1682/18.

## References

  * [1] Kostya S Novoselov et al.  “Two-dimensional atomic crystals”  In _Proceedings of the National Academy of Sciences_ 102.30 National Acad Sciences, 2005, pp. 10451–10453
  * [2] David Nelson, Tsvi Piran and Steven Weinberg  “Statistical mechanics of membranes and surfaces”  Singapore: World Scientific, 2004
  * [3] Yacov Kantor, Mehran Kardar and David R. Nelson  “Statistical Mechanics of Tethered Surfaces”  In _Phys. Rev. Lett._ 57 American Physical Society, 1986, pp. 791–794  DOI: 10.1103/PhysRevLett.57.791
  * [4] Yacov Kantor, Mehran Kardar and David R. Nelson  “Tethered surfaces: Statics and dynamics”  In _Phys. Rev. A_ 35 American Physical Society, 1987, pp. 3056–3071  DOI: 10.1103/PhysRevA.35.3056
  * [5] Yacov Kantor and David R. Nelson  “Crumpling transition in polymerized membranes”  In _Phys. Rev. Lett._ 58 American Physical Society, 1987, pp. 2774–2777  DOI: 10.1103/PhysRevLett.58.2774
  * [6] Maya Paczuski, Mehran Kardar and David R. Nelson  “Landau Theory of the Crumpling Transition”  In _Phys. Rev. Lett._ 60 American Physical Society, 1988, pp. 2638–2640  DOI: 10.1103/PhysRevLett.60.2638
  * [7] DR Nelson and L Peliti  “Fluctuations in membranes with crystalline and hexatic order”  In _J. Phys. France_ 48.7 Société française de physique, 1987, pp. 1085–1092
  * [8] Joseph Aronovitz, Leonardo Golubovic and TC Lubensky  “Fluctuations and lower critical dimensions of crystalline membranes”  In _Journal de Physique_ 50.6 Société française de physique, 1989, pp. 609–631
  * [9] Bertrand Duplantier, Terence Hwa and Mehran Kardar  “Self-avoiding crumpled manifolds: Perturbative analysis and renormalizability”  In _Phys. Rev. Lett._ 64 American Physical Society, 1990, pp. 2022–2025  DOI: 10.1103/PhysRevLett.64.2022
  * [10] Pierre Le Doussal and Leo Radzihovsky  “Self-consistent theory of polymerized membranes”  In _Physical review letters_ 69.8 APS, 1992, pp. 1209
  * [11] Andrej Košmrlj and David R. Nelson  “Mechanical properties of warped membranes”  In _Phys. Rev. E_ 88 American Physical Society, 2013, pp. 012136  DOI: 10.1103/PhysRevE.88.012136
  * [12] Andrej Košmrlj and David R. Nelson  “Thermal excitations of warped membranes”  In _Phys. Rev. E_ 89 American Physical Society, 2014, pp. 022126  DOI: 10.1103/PhysRevE.89.022126
  * [13] Andrej Košmrlj and David R. Nelson  “Response of thermalized ribbons to pulling and bending”  In _Phys. Rev. B_ 93 American Physical Society, 2016, pp. 125431  DOI: 10.1103/PhysRevB.93.125431
  * [14] Andrej Košmrlj and David R. Nelson  “Statistical Mechanics of Thin Spherical Shells”  In _Phys. Rev. X_ 7 American Physical Society, 2017, pp. 011002  DOI: 10.1103/PhysRevX.7.011002
  * [15] Pierre Le Doussal and Leo Radzihovsky  “Anomalous elasticity, fluctuations and disorder in elastic membranes”  In _Annals of Physics_ 392 Elsevier, 2018, pp. 340–410
  * [16] Suraj Shankar and David R. Nelson  “Thermalized buckling of isotropically compressed thin sheets”  In _Phys. Rev. E_ 104 American Physical Society, 2021, pp. 054141  DOI: 10.1103/PhysRevE.104.054141
  * [17] Ali Morshedifard, Miguel Ruiz-García, Mohammad Javad Abdolhosseini Qomi and Andrej Košmrlj  “Buckling of thermalized elastic sheets”  In _Journal of the Mechanics and Physics of Solids_ 149 Elsevier, 2021, pp. 104296
  * [18] Mehran Kardar and David R. Nelson  “$\epsilon$ expansions for crumpled manifolds”  In _Phys. Rev. Lett._ 58 American Physical Society, 1987, pp. 1289–1292  DOI: 10.1103/PhysRevLett.58.1289
  * [19] Fatemeh Ahmadpoor, Peng Wang, Rui Huang and Pradeep Sharma  “Thermal fluctuations and effective bending stiffness of elastic thin sheets and graphene: A nonlinear analysis”  In _Journal of the Mechanics and Physics of Solids_ 107 Elsevier, 2017, pp. 294–319
  * [20] KR Elder et al.  “Modeling buckling and topological defects in stacked two-dimensional layers of graphene and hexagonal boron nitride”  In _Physical Review Materials_ 5.3 APS, 2021, pp. 034004
  * [21] Enzo Granato et al.  “Dynamic scaling of out-of-plane fluctuations in freestanding graphene”  In _Physical Review B_ 105.20 APS, 2022, pp. L201409
  * [22] Jannik C Meyer et al.  “The structure of suspended graphene sheets”  In _Nature_ 446.7131 Nature Publishing Group, 2007, pp. 60–63
  * [23] Jannik C Meyer et al.  “On the roughness of single-and bi-layer graphene membranes”  In _Solid State Communications_ 143.1-2 Elsevier, 2007, pp. 101–109
  * [24] Franck Plouraboué and Stéphane Roux  “Experimental study of the roughness of crumpled surfaces”  In _Physica A_ 227, 1996, pp. 173–182
  * [25] Daniel L. Blair and Arshad Kudrolli  “Geometry of Crumpled Paper”  In _Phys. Rev. Lett._ 94 American Physical Society, 2005, pp. 166107  DOI: 10.1103/PhysRevLett.94.166107
  * [26] Alexander S. Balankin et al.  “Intrinsically anomalous roughness of randomly crumpled thin sheets”  In _Phys. Rev. E_ 74 American Physical Society, 2006, pp. 061602  DOI: 10.1103/PhysRevE.74.061602
  * [27] Alexander S. Balankin and Orlando Susarrey Huerta  “Entropic rigidity of a crumpling network in a randomly folded thin sheet”  In _Phys. Rev. E_ 77 American Physical Society, 2008, pp. 051124  DOI: 10.1103/PhysRevE.77.051124
  * [28] Alexander S. Balankin et al.  “Fractal features of a crumpling network in randomly folded thin matter and mechanics of sheet crushing”  In _Phys. Rev. E_ 87 American Physical Society, 2013, pp. 052806  DOI: 10.1103/PhysRevE.87.052806
  * [29] Stefan Bilbao  “A family of conservative finite difference schemes for the dynamical von Karman plate equations”  In _Numerical Methods for Partial Differential Equations: An International Journal_ 24.1 Wiley Online Library, 2008, pp. 193–216
  * [30] Stefan Bilbao, Olivier Thomas, Cyril Touzé and Michele Ducceschi  “Conservative numerical methods for the full von Kármán plate equations”  In _Numerical Methods for Partial Differential Equations_ 31.6 Wiley Online Library, 2015, pp. 1948–1970
  * [31] Mark J Bowick et al.  “The flat phase of crystalline membranes”  In _Journal de Physique I_ 6.10 EDP Sciences, 1996, pp. 1321–1345
  * [32] Annalisa Fasolino, JH Los and Mikhail I Katsnelson  “Intrinsic ripples in graphene”  In _Nature materials_ 6.11 Nature Publishing Group, 2007, pp. 858–861
  * [33] JH Los et al.  “Scaling properties of flexible membranes from atomistic simulations: application to graphene”  In _Physical Review B_ 80.12 APS, 2009, pp. 121405
  * [34] Andreas Tröster  “High-precision Fourier Monte Carlo simulation of crystalline membranes”  In _Physical Review B_ 87.10 APS, 2013, pp. 104112
  * [35] Andreas Tröster  “Fourier Monte Carlo simulation of crystalline membranes in the flat phase”  In _Journal of Physics: Conference Series_ 454.1, 2013, pp. 012032  IOP Publishing
  * [36] Rebecca C Thompson-Flagg, Maria JB Moura and M Marder  “Rippling of graphene”  In _EPL (Europhysics Letters)_ 85.4 IOP Publishing, 2009, pp. 46002
  * [37] Siby Thomas, KM Ajith, Sharat Chandra and MC Valsakumar  “Temperature dependent structural properties and bending rigidity of pristine and defective hexagonal boron nitride”  In _Journal of Physics: Condensed Matter_ 27.31 IOP Publishing, 2015, pp. 315302
  * [38] JH Los, A Fasolino and MI Katsnelson  “Scaling behavior and strain dependence of in-plane elastic properties of graphene”  In _Physical review letters_ 116.1 APS, 2016, pp. 015901
  * [39] Erwin Frey and David L. Nelson  “Dynamics of flat membranes and flickering in red blood cells”  In _Journal de Physique I_ 1.12 EDP Sciences, 1991, pp. 1715–1757  DOI: 10.1051/jp1:1991238
  * [40] Ken-ichi Mizuochi, Hiizu Nakanishi and Takahiro Sakaue  “Dynamical scaling of polymerized membranes”  In _EPL (Europhysics Letters)_ 107.3 IOP Publishing, 2014, pp. 38003
  * [41] Chanania Steinbock, Eytan Katzav and Arezki Boudaoud  “Structure of fluctuating thin sheets under random forcing”  In _Phys. Rev. Research_ 4 American Physical Society, 2022, pp. 033096  DOI: 10.1103/PhysRevResearch.4.033096
  * [42] Chanania Steinbock and Eytan Katzav  “Dynamics of fluctuating thin sheets under random forcing”  In _Phys. Rev. E_ 107 American Physical Society, 2023, pp. 025002  DOI: 10.1103/PhysRevE.107.025002
  * [43] P.. Hohenberg and B.. Halperin  “Theory of dynamic critical phenomena”  In _Rev. Mod. Phys._ 49 American Physical Society, 1977, pp. 435–479  DOI: 10.1103/RevModPhys.49.435
  * [44] Tamás Vicsek and Fereydoon Family  “Dynamic scaling for aggregation of clusters”  In _Physical Review Letters_ 52.19 APS, 1984, pp. 1669
  * [45] Fereydoon Family and Tamas Vicsek  “Scaling of the active zone in the Eden process on percolation networks and the ballistic deposition model”  In _Journal of Physics A: Mathematical and General_ 18.2 IOP Publishing, 1985, pp. L75
  * [46] A-L Barabási and Harry Eugene Stanley  “Fractal concepts in surface growth”  Cambridge, England: Cambridge university press, 1995
  * [47] Eytan Katzav and Moshe Schwartz  “Exponent Inequalities in Dynamical Systems”  In _Phys. Rev. Lett._ 107 American Physical Society, 2011, pp. 125701  DOI: 10.1103/PhysRevLett.107.125701
  * [48] Eytan Katzav and Moshe Schwartz  “Dynamical inequality in growth models”  In _EPL (Europhysics Letters)_ 95.6 IOP Publishing, 2011, pp. 66003  DOI: 10.1209/0295-5075/95/66003
  * [49] Gustavo Düring, Christophe Josserand and Sergio Rica  “Weak Turbulence for a Vibrating Plate: Can One Hear a Kolmogorov Spectrum?”  In _Phys. Rev. Lett._ 97 American Physical Society, 2006, pp. 025503  DOI: 10.1103/PhysRevLett.97.025503
  * [50] Arezki Boudaoud, Olivier Cadot, Benoît Odille and Cyril Touzé  “Observation of Wave Turbulence in Vibrating Plates”  In _Phys. Rev. Lett._ 100 American Physical Society, 2008, pp. 234504  DOI: 10.1103/PhysRevLett.100.234504
  * [51] Nicolas Mordant  “Are There Waves in Elastic Wave Turbulence?”  In _Phys. Rev. Lett._ 100 American Physical Society, 2008, pp. 234505  DOI: 10.1103/PhysRevLett.100.234505
  * [52] Olivier Cadot, Arezki Boudaoud and Cyril Touzé  “Statistics of power injection in a plate set into chaotic vibration”  In _The European Physical Journal B_ 66.3 Springer, 2008, pp. 399–407
  * [53] Pablo Cobelli et al.  “Space-Time Resolved Wave Turbulence in a Vibrating Plate”  In _Phys. Rev. Lett._ 103 American Physical Society, 2009, pp. 204301  DOI: 10.1103/PhysRevLett.103.204301
  * [54] Thomas Humbert et al.  “Wave turbulence in vibrating plates: the effect of damping”  In _EPL (Europhysics Letters)_ 102.3 IOP Publishing, 2013, pp. 30002
  * [55] Benjamin Miquel, Alexandros Alexakis, Christophe Josserand and Nicolas Mordant  “Transition from Wave Turbulence to Dynamical Crumpling in Vibrated Elastic Plates”  In _Phys. Rev. Lett._ 111 American Physical Society, 2013, pp. 054302  DOI: 10.1103/PhysRevLett.111.054302
  * [56] Gustavo Düring, Christophe Josserand and Sergio Rica  “Self-similar formation of an inverse cascade in vibrating elastic plates”  In _Phys. Rev. E_ 91 American Physical Society, 2015, pp. 052916  DOI: 10.1103/PhysRevE.91.052916
  * [57] Gustavo Düring, Christophe Josserand and Sergio Rica  “Wave turbulence theory of elastic plates”  In _Physica D: Nonlinear Phenomena_ 347 Elsevier, 2017, pp. 42–73
  * [58] Gustavo Düring, Christophe Josserand, Giorgio Krstulovic and Sergio Rica  “Strong turbulence for vibrating plates: Emergence of a Kolmogorov spectrum”  In _Physical Review Fluids_ 4.6 APS, 2019, pp. 064804
  * [59] Roumaissa Hassaini et al.  “Elastic weak turbulence: From the vibrating plate to the drum”  In _Phys. Rev. E_ 99 American Physical Society, 2019, pp. 033002  DOI: 10.1103/PhysRevE.99.033002
  * [60] Ignacio Pavez and Gustavo Düring  “Path integral description and direct interaction approximation for elastic plate turbulence”  In _Chaos, Solitons & Fractals_ 166 Elsevier, 2023, pp. 112911
  * [61] Lev Davidovich Landau, Evgenii Mikhailovich Lifshitz, Arnol Markovich Kosevich and Lev Petrovich Pitaevskii  “Theory of elasticity”  New York: Elsevier, 1986
  * [62] Alex Lobkovsky et al.  “Scaling properties of stretching ridges in a crumpled elastic sheet”  In _Science_ 270.5241 American Association for the Advancement of Science, 1995, pp. 1482–1485
  * [63] Alexander E. Lobkovsky and T.. Witten  “Properties of ridges in elastic membranes”  In _Phys. Rev. E_ 55 American Physical Society, 1997, pp. 1577–1589  DOI: 10.1103/PhysRevE.55.1577
  * [64] M Ben Amar and Y Pomeau  “Crumpled paper”  In _Proceedings of the Royal Society of London. Series A: Mathematical, Physical and Engineering Sciences_ 453.1959 The Royal Society, 1997, pp. 729–755
  * [65] Thomas A Witten  “Stress focusing in elastic sheets”  In _Reviews of Modern Physics_ 79.2 APS, 2007, pp. 643
  * [66] Francesco Cagnetta, Viktor Škultéty, Martin R Evans and Davide Marenduzzo  “Universal properties of active membranes”  In _Physical Review E_ 105.1 APS, 2022, pp. L012604
  * [67] Hagen Kleinert and Verena Schulte-Frohlinde  “Critical properties of $\phi$4-theories”  Singapore: World Scientific, 2001
  * [68] Hannes Risken  “Fokker-planck equation”  In _The Fokker-Planck Equation_ New York: Springer, 1996, pp. 63–95
  * [69] L Isserlis  “On certain probable errors and correlation coefficients of multiple frequency distributions with skew regression”  In _Biometrika_ 11.3 JSTOR, 1916, pp. 185–190
  * [70] Leon Isserlis  “On a formula for the product-moment coefficient of any order of a normal frequency distribution in any number of variables”  In _Biometrika_ 12.1/2 JSTOR, 1918, pp. 134–139
  * [71] Mehran Kardar  “Statistical physics of particles”  Cambridge, England: Cambridge University Press, 2007
  * [72] George A Baker and Peter Graves-Morris  “Pade Approximants: Encyclopedia of Mathematics and It’s Applications, Vol. 59”  Cambridge, England: Cambridge University Press, 1996
  * [73] Eytan Katzav and Moshe Schwartz  “Kardar-Parisi-Zhang equation with temporally correlated noise: A self-consistent approach”  In _Phys. Rev. E_ 70 American Physical Society, 2004, pp. 011601  DOI: 10.1103/PhysRevE.70.011601
  * [74] Henry Poincaré  “Les Méthodes Nouvelles de la Mécanique Célèste, vol. II”  Dover, 1957 [1893], pp. $\S$123–$\S$128
  * [75] A. Lindstedt  In _Abh. K. Akad. Wiss. St. Petersburg_ 31, 1882
  * [76] P.G. Drazin  “Nonlinear systems”  Cambridge, England: Cambridge University Press, 1992, pp. 181–186
  * [77] C… Rao, A.. Sood, K.. Subrahmanyam and Achutharao Govindaraj  “Graphene: the new two-dimensional nanomaterial”  In _Angewandte Chemie International Edition_ 48.42 Wiley Online Library, 2009, pp. 7752–7777
  * [78] Jae-Ung Lee, Duhee Yoon and Hyeonsik Cheong  “Estimation of Young’s modulus of graphene by Raman spectroscopy”  In _Nano letters_ 12.9 ACS Publications, 2012, pp. 4444–4448
  * [79] Yujie Wei et al.  “Bending rigidity and Gaussian bending stiffness of single-layered graphene”  In _Nano letters_ 13.1 ACS Publications, 2013, pp. 26–30
  * [80] Laura Picas, Felix Rico and Simon Scheuring  “Direct measurement of the mechanical properties of lipid phases in supported bilayers”  In _Biophysical journal_ 102.1 Elsevier, 2012, pp. L01–L03
  * [81] Moshe Schwartz and SF Edwards  “Nonlinear deposition: A new approach”  In _EPL (Europhysics Letters)_ 20.4 IOP Publishing, 1992, pp. 301
  * [82] Moshe Schwartz and S.. Edwards  “Peierls-Boltzmann equation for ballistic deposition”  In _Phys. Rev. E_ 57 American Physical Society, 1998, pp. 5730–5739  DOI: 10.1103/PhysRevE.57.5730
  * [83] Eytan Katzav and Moshe Schwartz  “Self-consistent expansion for the Kardar-Parisi-Zhang equation with correlated noise”  In _Phys. Rev. E_ 60 American Physical Society, 1999, pp. 5677–5680  DOI: 10.1103/PhysRevE.60.5677
  * [84] Eytan Katzav and Moshe Schwartz  “Existence of the upper critical dimension of the Kardar–Parisi–Zhang equation”  In _Physica A: Statistical Mechanics and its Applications_ 309.1-2 Elsevier, 2002, pp. 69–78
  * [85] Moshe Schwartz and SF Edwards  “Stretched exponential in non-linear stochastic field theories”  In _Physica A: Statistical Mechanics and its Applications_ 312.3-4 Elsevier, 2002, pp. 363–368
  * [86] Eytan Katzav  “Self-consistent expansion for the molecular beam epitaxy equation”  In _Phys. Rev. E_ 65 American Physical Society, 2002, pp. 032103  DOI: 10.1103/PhysRevE.65.032103
  * [87] Eytan Katzav  “Self-consistent expansion results for the nonlocal Kardar-Parisi-Zhang equation”  In _Phys. Rev. E_ 68 American Physical Society, 2003, pp. 046113  DOI: 10.1103/PhysRevE.68.046113
  * [88] Eytan Katzav and Moshe Schwartz  “Numerical evidence for stretched exponential relaxations in the Kardar-Parisi-Zhang equation”  In _Phys. Rev. E_ 69 American Physical Society, 2004, pp. 052603  DOI: 10.1103/PhysRevE.69.052603
  * [89] Sam F Edwards and Moshe Schwartz  “Lagrangian statistical mechanics applied to non-linear stochastic field equations”  In _Physica A: Statistical Mechanics and its Applications_ 303.3-4 Elsevier, 2002, pp. 357–386
  * [90] E Katzav and M Adda-Bedia  “Roughness of tensile crack fronts in heterogenous materials”  In _EPL (Europhysics Letters)_ 76.3 IOP Publishing, 2006, pp. 450
  * [91] E. Katzav, M. Adda-Bedia, M. Ben Amar and A. Boudaoud  “Roughness of moving elastic lines: Crack and wetting fronts”  In _Phys. Rev. E_ 76 American Physical Society, 2007, pp. 051601  DOI: 10.1103/PhysRevE.76.051601
  * [92] Moshe Schwartz and Eytan Katzav  “The ideas behind self-consistent expansion”  In _Journal of Statistical Mechanics: Theory and Experiment_ IOP Publishing, 2008, pp. P04023
  * [93] Benjamin Remez and Moshe Goldstein  “From divergent perturbation theory to an exponentially convergent self-consistent expansion”  In _Phys. Rev. D_ 98 American Physical Society, 2018, pp. 056017  DOI: 10.1103/PhysRevD.98.056017
  * [94] E Katzav, M Adda-Bedia and B Derrida  “Fracture surfaces of heterogeneous materials: a 2D solvable model”  In _EPL (Europhysics Letters)_ 78.4 IOP Publishing, 2007, pp. 46006
  * [95] E Katzav and M Adda-Bedia  “Stability and roughness of tensile cracks in disordered materials”  In _Physical Review E_ 88.5 APS, 2013, pp. 052402
  * [96] Vahideh Farzam Rad et al.  “Speckle pattern analysis of crumpled papers”  In _Applied optics_ 58.24 Optical Society of America, 2019, pp. 6549–6554
  * [97] Eric M. Kramer and Alexander E. Lobkovsky  “Universal power law in the noise from a crumpled elastic sheet”  In _Phys. Rev. E_ 53 American Physical Society, 1996, pp. 1465–1469  DOI: 10.1103/PhysRevE.53.1465
  * [98] Paul A. Houle and James P. Sethna  “Acoustic emission from crumpling paper”  In _Phys. Rev. E_ 54 American Physical Society, 1996, pp. 278–283  DOI: 10.1103/PhysRevE.54.278
  * [99] SA Safran  “Curvature elasticity of thin films”  In _Advances in physics_ 48.4 Taylor & Francis, 1999, pp. 395–448
  * [100] Giuseppe Pellicane, Mireille Megnidio-Tchoukouegno, Genene T Mola and Mesfin Tsige  “Surface enrichment driven by polymer topology”  In _Physical Review E_ 93.5 APS, 2016, pp. 050501