# Transitions in Stochastic Non-equilibrium Systems:
Efficient Reduction and Analysis

Mickaël D. Chekroun  Department of Atmospheric and Oceanic Sciences,
University of California, Los Angeles, CA 90095-1565, USA, and
Department of Earth and Planetary Sciences, Weizmann Institute, Rehovot 76100,
Israel<EMAIL_ADDRESS>,  Honghu Liu  Department of Mathematics,
Virginia Tech, Blacksburg, VA 24061, USA<EMAIL_ADDRESS>,  James C. McWilliams
Department of Atmospheric and Oceanic Sciences, University of California, Los
Angeles, CA 90095-1565, USA<EMAIL_ADDRESS>and  Shouhong Wang  Department
of Mathematics, Indiana University, Bloomington, IN 47405<EMAIL_ADDRESS>
###### Abstract.

A central challenge in physics is to describe non-equilibrium systems driven
by randomness, such as a randomly growing interface, or fluids subject to
random fluctuations that account e.g. for local stresses and heat fluxes in
the fluid which are not related to the velocity and temperature gradients. For
deterministic systems with infinitely many degrees of freedom, normal form and
center manifold theory have shown a prodigious efficiency to often completely
characterize how the onset of linear instability translates into the emergence
of nonlinear patterns, associated with genuine physical regimes. However, in
presence of random fluctuations, the underlying reduction principle to the
center manifold is seriously challenged due to large excursions caused by the
noise, and the approach needs to be revisited.

In this study, we present an alternative framework to cope with these
difficulties exploiting the approximation theory of stochastic invariant
manifolds, on one hand, and energy estimates measuring the defect of
parameterization of the high-modes, on the other. To operate for fluid
problems subject to stochastic stirring forces, these error estimates are
derived under assumptions regarding dissipation effects brought by the high-
modes in order to suitably counterbalance the loss of regularity due to the
nonlinear terms. As a result, the approach enables us to predict, from reduced
equations of the stochastic fluid problem, the occurrence in large probability
of a stochastic analogue to the pitchfork bifurcation, as long as the noise’s
intensity and the eigenvalue’s magnitude of the mildly unstable mode scale
accordingly.

In the case of SPDEs forced by a multiplicative noise in the orthogonal
subspace of e.g. its mildly unstable mode, our parameterization formulas show
that the noise gets transmitted to this mode via non-Markovian coefficients,
and that the reduced equation is only stochastically driven by the latter.
These coefficients depend explicitly on the noise path’s history, and their
memory content is self-consistently determined by the intensity of the random
force and its interaction through the SPDE’s nonlinear terms. Applications to
a stochastic Rayleigh-Bénard problem are detailed, for which conditions for a
stochastic pitchfork bifurcation (in large probability) to occur, are
clarified.

###### Contents

  1. 1 Introduction
  2. 2 Reduced-order dynamics of stochastic PDEs near the onset of instability
    1. 2.1 SPDEs driven by multiplicative noise
    2. 2.2 Leading-order approximation of stochastic invariant manifolds, and non-Markovian terms
    3. 2.3 Tracking the large-excursions through the non-Markovian $M$-terms
  3. 3 Modeling error estimates from reduced equations
    1. 3.1 Assumptions and pitchfork scenario
    2. 3.2 Preparatory estimates for the pitchfork scenario
    3. 3.3 Residual error estimates in large probability
  4. 4 Low- and high-mode error estimates
    1. 4.1 A priori bounds on SPDE solutions
    2. 4.2 Fluid problems subject to fluctuations: Low- and high-mode error estimates
  5. 5 Applications to Rayleigh-Bénard convection
    1. 5.1 The stochastic Rayleigh-Bénard model and its mathematical formulation
    2. 5.2 Verification of Conditions 1-3 and 1
    3. 5.3 Verification of Conditions 1 and 2
    4. 5.4 Stochastic pitchfork bifurcation in large probability
  6. 6 Concluding remarks
  7. A $M$-term: First moments, non-Gaussian statistics, and large probability estimates
    1. A.1 First moments
    2. A.2 Non-Gaussian statistics and log-normal approximations
    3. A.3 Probability estimates in the large
  8. B A priori estimates in large probability: Surrogate system (3.19)
  9. C High-mode estimates in large probability
    1. C.1 High-mode estimate (4.7)
    2. C.2 High-mode estimate (4.9)
  10. D The non-Markovian normal form of a pitchfork bifurcation
  11. E Proof of Lemma 5.1

##  1\. Introduction

The development of instability in a physical system may often be described in
terms of the temporal evolution of the amplitudes of certain normal modes,
namely those that are mildly unstable and those that are only slightly damped
in linear theory. When the number of these nearly marginal modes is finite,
their amplitudes are governed by ordinary differential equations (ODEs) in
which the growth rates of the linear theory have been renormalized by
nonlinear terms [Cra91, CK91]. Intuitively, the reason for this reduction is a
simple separation of time scales. Modes that have just crossed the imaginary
axis have a small real part and are evolving on long time scales, all the
other fast modes rapidly adapting themselves to the slow modes.

For deterministic dynamical physical systems with infinitely many degrees of
freedom (d.o.f.), normal form and center manifold theory [Hen81, MW05] along
with their recent extension [MW14] have shown a prodigious efficiency to often
completely characterize how the onset of linear instability translates into
the emergence of nonlinear patterns, associated with genuine physical regimes.
To do so, the theory along with its recent advances identify, for a given
problem, a nonlinear reduction mapping that permits to reduce the many d.o.f.
to a few essential variables (governed by reduced ODEs), able to predict in
turn the dynamical transitions [MW14].

However, in presence of random fluctuations, the reduced ODEs produced by
these deterministic theories are no longer valid objects to account for
phenomena pertaining to the stochastic realm such as noise-induced
transitions. In the stochastic context, the reduction mapping needs indeed to
account for these random fluctuations. This mapping consists of e.g. a
parameterization of the stable modes which becomes then stochastic due to
noise. The existence of such reduction mappings has been abundantly studied
under various settings for stochastic partial differential equations (SPDEs)
[BF95, DPD96, DLS03, DLS04], but the question of their efficient calculation
or approximation has been much less addressed, mainly in some special cases
[BHP07] or by means of deterministic (and thus incomplete) formulas [BW10].

Only in the recent years, the derivation of general, stochastic, approximation
formulas to the leading-order (in the nonlinear terms) of stochastic invariant
manifolds (IMs) appeared in [CLW15a] for a broad class of SPDEs subject to a
multiplicative noise. Yet, in spite of their relevance demonstrated through
numerical examples [CLW15b], the usage of such formulas to derive reduced
stochastic differential equations (SDEs) aimed at predicting dynamical
transitions in SPDEs, calls for more understanding.111In fact, the bifurcation
analysis for SPDEs is available only for a few particular cases [CLR01,
Wan14]. Even for stochastic ODEs, the question of how to describe a stochastic
bifurcation is not completely settled; see e.g. [Arn98, Chap. 9] and [CF98].

Indeed, what makes the success of center/unstable manifold theory to describe
bifurcations for PDEs—the so-called reduction principle [Cra91, SY02]—from
low-dimensional reduced systems fails to operate in the stochastic realm,
calling for a new approach to apprehend bifurcation analysis in presence of
noise. The reason is that a breakdown of the reduction principle to the
center/unstable manifold now takes place in the stochastic setting, and an
exponential slaving of the stable modes onto the unstable ones is no longer
valid, due to large excursions caused by the (white) noise, even when its
intensity is small. At a technical level, the exponential attraction by
invariant manifolds that holds in presence of a sufficiently large spectral
gap and that extends to SPDEs with globally Lipschitz nonlinearity [CLW15a,
Corollary 4.3], is violated in the stochastic setting for SPDEs with
nonlinearities that are only locally Lipschitz, such as encountered in
stochastic fluid problems. In these cases, the solution is not guaranteed to
stay within the neighborhood over which the local stochastic IMs exist, due to
noise-induced deviations.

In this article we bring justification regarding the usage of reduced SDEs
built upon stochastic IMs and their approximations formulas based on [CLW15a]
to predict dynamical transitions in SPDEs, via estimates measuring the error
between the solutions to the reduced SDEs and those of the original SPDE. More
generally, the goal of this article is twofold: (i) To advance the fundamental
understanding for the efficient derivation of reduced models able to describe
spatiotemporal pattern formations and their transitions subject to random
fluctuations in non-equilibrium systems described by SPDEs, and (ii) to derive
error estimates relevant for the analysis of the spatio-temporal transitions,
in a random environment.

Other approaches have shown their usefulness to derive low-dimensional reduced
systems able to capture the essential macroscopic dynamics of SPDEs. Examples
include the averaging method [CF09, Kif04, PS08, WR13], the amplitude equation
approach [BH04, Blö07, BF20, BHP07], and approaches rooted in the singular
perturbation theory of Markov processes [MTVE01, Kur73].

However, it remains challenging to derive reduced systems and rigorous error
estimates to predict dynamical transitions in e.g. stochastic fluid problems
involving typically quadratic nonlinearities. For such problems, the main
issue lies in the fact that the nonlinearity does not act only on the dominant
modes but also influences them through non-dominant modes [Blö05, BHP07], and
thus makes central the way the noise interacts with the nonlinear dynamics,
either forcing the dominant or the neglected modes, or both. This is where, as
explained below, the stochastic reduction mappings provided by the
approximation theory of stochastic IMs [CLW15a] show their usefulness. As we
will see, these mappings parameterize the neglected variables beyond Ornstein-
Uhlenbeck approximations (cf. [BH05, MTVE01, BM09]) by producing essential
variables able to predict dynamical transitions (with large probability) for a
broad class of SPDEs with quadratic nonlinearities (such as arising in fluid
problems), regardless of whether the noise acts only on the resolved or
unresolved scales, or both (cf. [BHP07]).

The main difficulty to perform an efficient stochastic reduction for SPDEs is
thus to track properly the noise-path dependence which in many situations
requires to be properly parameterized to render account for noise-induced
phenomena and other interactions between the noise and nonlinear dynamics.
This is the purpose of Theorem 2.1 recalled below that summarizes from
[CLW15a] general formulas of stochastic reduction mappings based on the
approximation results of [CLW15a, Thm. 6.1 & Cor. 6.1]. These results entail
that the corresponding stochastic reduction mappings not only provide the
leading-order approximation of the underlying stochastic IM but also the
tracking of its time-stochastic dependence via explicit random coefficients;
see (2.17) below. As already pointed out in [CLW15b, Chap. 5] and more
substantiated here, these coefficients are non-Markovian as dependent
explicitly on the history of the noise path, conveying thus exogeneous memory
effects [HO07, Hai09] rooted in the (very) small-scale fluctuations
(noise).222These exogeneous memory effects are to be contrasted with the
endogeneous ones encountered in the reduction of nonlinear autonomous systems
as predicted by the Mori-Zwanzig (MZ) theory; see e.g. [CHK02, GKS04, ZDV18,
FQS+19, SGLCG21]. The former are functionals of the “past” of the noise and
emerge in the reduction of stochastic systems by means of stochastic IMs
[CLW15b] whereas the latter are functionals of the past of the resolved
variables that arise in the reduction of autonomous systems when the validity
of (deterministic) IMs or, more generally, of the conditional expectation
alone in the MZ expansion, breaks down; see [Sti06, CLM17, CLM20, CLM21].

Mathematically, these coefficients are shown to correspond to stationary
solutions of auxiliary SDEs whose drift part depends on the distance to
certain resonances between the critical eigenvalues (losing stability) and the
stable part of the spectrum; see (2.19) below. Such a distance to resonance
controls in turns the decay of temporal correlations of these coefficients, as
well as their statistics, allowing for non-Gaussian ones—unlike Ornstein-
Uhlenbeck processes encountered in other reduction methods [BH05, MTVE01,
BM09]—with tails all the more pronounced than their ratio with the noise
amplitude approaches unity; see Appendix A.

As already shown through numerical examples in [CLW15b], these non-Markovian
coefficients constitute in fact the key ingredients for the reduced SDEs to
achieve good performance, and can be seen as part of the essential variables
produced by the underlying stochastic reduction mapping. For instance, in the
case when the noise acts only on the (unresolved) stable modes, such a
reduction mapping shows that the noise gets transmitted to the (resolved)
unstable modes via these non-Markovian coefficients, and the reduced equation
is only stochastically driven by the latter. In contrast, any deterministic
parameterization of the unresolved modes will lead in this case to a
deterministic ODE system, which is insufficient to capture any noise-induced
phenomena in the original system. In Sec. 2.3 below, we illustrate on a simple
two-dimensional system, the ability to capture its stochastic dynamics and its
large excursions from such non-Markovian reduced equations.

In this article, we go beyond these numerical examples and show that these
non-Markovian coefficients play furthermore a key role to conduct a rigorous
error analysis thanks to the auxiliary SDEs they satisfy. Our focus is on
SPDEs involving energy-preserving quadratic nonlinearities motivated by
physical problems such as the prediction and characterization of pattern
formation arising in fluid problems after the onset of instability and subject
to noise disturbance. In order to circumscribe the difficulties while
conveying the main ideas, we place ourselves in the case of a supercritical
pitchfork bifurcation, assuming that exactly one eigenmode loses its stability
at the critical parameter value. An attentive reader will notice that the
framework laid out here is actually not limited to the pitchfork scenario and
adaptable to many other bifurcation/transition scenarios, including the Hopf
bifurcation and other attractor bifurcations.

At the basis of our error estimates is a “Markovianisation” of the one-
dimensional non-Markovian reduced SDE governing the local dynamics in the case
of a pitchfork bifurcation scenario; cf. (3.8) and (3.19) below. In this
reformulation, the dynamical equation satisfied by the stochastic
parameterization $\Phi(t)$ is obtained by Itô’s formula, which is made
explicit thanks to the aforementioned auxiliary SDEs determining the random
coefficients involved in $\Phi$; see again (2.19). By recasting the original
SPDE into a system of two coupled equations satisfied by the unstable and
stable components of the SPDE solution (see (3.23)), the estimates become
facilitated, exploiting in particular the dynamical equation satisfied by
$\Phi(t)$ (in (3.19)) aimed at parameterizing the stable component of the SPDE
solution.

The error estimates assessing the quality of our non-Markovian reduced SDEs to
predict dynamical transitions in SPDEs, are then organized through Secns. 3
and 4 as follows. After establishing in Sec. 3.2 a suitable a priori bound for
the solution $(X(t),\Phi(t))$ of the Markovianized reduced system (3.19), we
derive in Sec. 3.3 a residual error estimate for the non-Markovian reduced SDE
(3.8); see Lemma 3.1 and Theorem 3.1 below. This residual error estimate
gives—with large probability and over long time intervals—the error made by
$(X(t),\Phi(t))$ in satisfying the SPDE model (3.23). Our result shows that
this error is all the smaller as the eigenvalue $\epsilon$ associated with the
unstable mode is small and the noise intensity scales as $\sqrt{\epsilon}$.

We then derive in Sec. 4 rigorous error estimates between the solutions of the
reduced systems and those of the original SPDE emanating from small initial
data. The result presented in Theorem 4.1 shows that the same bounds
controlling the residual error in Theorem 3.1 control also in large
probability and over long time intervals, the error between the solutions of
the reduced systems and those of the SPDE.

The derivation of this error estimate requires suitable a priori bounds of the
SPDE solution itself. As explained in Sec. 4.1 below, the main difficulty in
deriving these a priori bounds compared to those for $(X(t),\Phi(t))$ lies in
the infinite-dimensional character of the projection onto the stable subspace
of the SPDE solution to which applies the loss of regularity via the
nonlinearity $B$, requiring de facto new estimates. In contrast, $\Phi(t)$ is
of finite-dimensional range due to our working assumptions; see condition
(3.4) below. This latter condition is typically encountered for fluid problems
defined e.g. on rectangular spatial domains and is made here to avoid
unnecessary technicalities. We show that with the addition of mild dissipative
assumptions on the linear operator $L_{\lambda}$, one can deal with more
general situations and counterbalance the loss of regularity induced by the
nonlinear term $B$ to derive thereof the required a priori bounds on the SPDE
solution; see Conditions 1-2 and Remark 4.1 below.

Equipped with these error estimates, the non-Markovian reduced SDE (3.8)
provides thus a natural normal form to describe in large probability the
notion of a stochastic pitchfork bifurcation for SPDEs. Indeed, thanks to the
error estimates of Theorem 4.1, the non-Markovian reduced SDEs derived from
the approximation formulas of stochastic IMs given in Theorem 2.1, provide an
effective way to predict stochastic pitchfork bifurcations (in large
probability) for a broad class of SPDEs issued from fluid dynamics, as long as
the noise intensity and the eigenvalue’s magnitude of the mildly unstable
mode, scale accordingly.

As an application, we illustrate in Sec. 5 below that our framework is able to
predict a stochastic pitchfork bifurcation in a stochastic Rayleigh-Bénard
model, in terms of natural conditions involving the distance of the Rayleigh
number to its first critical value; see Theorem 5.1. The article is then
concluded by Sec. 6 with a discussion about future developments of the
framework introduced here.

##  2\. Reduced-order dynamics of stochastic PDEs near the onset of
instability

###  2.1. SPDEs driven by multiplicative noise

In [CLW15a, CLW15b] was undertaken the endeavor of deriving general, analytic
formulas for the approximation of (local) stochastic invariant manifolds of
the following class of nonlinear stochastic evolution equations posed on a
Hilbert space $H$ and driven by a multiplicative noise according to

$\mathrm{d}u=\big{(}L_{\lambda}u+F(u)\big{)}\mathrm{d}t+\sigma
u\circ\mathrm{d}W_{t},\;u\in H,$ (2.1)

where $\sigma u\circ\mathrm{d}W_{t}$ means that this stochastic term (of
magnitude $\sigma>0$) is considered in the Stratonovich sense. Here $W_{t}$
denotes a standard Brownian motion whose probability space is denoted by
$\Omega$, and whose probability law is denoted by $\mathbb{P}$. A noise’s
realization, also called a noise path will be thus labeled by an element
$\omega$ in $\Omega$. This multiplicative term is such that the perturbation
“scales” with the solution, such that the total energy is e.g. (formally)
$\mathbb{P}$-almost surely preserved in the case of a nonlinearity $F$ that is
energy preserving.

Such equations results from recasting an SPDE with its boundary conditions, in
which the unknown $u$ evolves in a functional space, taken to be here a
Hilbert space $H$. The operator $L_{\lambda}$ represents a linear differential
operator (parameterized by a scalar $\lambda$) while $F$ is a nonlinear
operator that accounts for the nonlinear terms. Both of these operators imply
loss of regularity when applied to a function $u$ in $H$. To have a consistent
existence theory of solutions and their stochastic invariant manifolds for
SPDEs requires to take into account these loss of regularity effects.

To do so, we place ourselves in the functional setting of sectorial operators
and analytic semigroups [Hen81, Paz83]. We assume that

$L_{\lambda}=-A+P_{\lambda},$

where $A$ is sectorial with domain $D(A)\subset H$ which is compactly and
densely embedded in $H$. We assume also that $-A$ is stable, while
$P_{\lambda}$ is a low-order perturbation of $A$, i.e.
$P_{\lambda}:D(A^{\alpha})\rightarrow H$ is a family of bounded linear
operators depending continuously on $\lambda$, for some $\alpha$ in $[0,1)$;
see [Hen81, Sec. 1.4] for definitions of fractional power of an operator and
the function space $D(A^{\alpha})$. In practice, the choice of $\alpha$ should
match the loss of regularity effects caused by the nonlinear terms so that

$F\colon D(A^{\alpha})\rightarrow H,$

is a well-defined $C^{p}$-smooth mapping that satisfies $F(0)=0$, $DF(0)=0$
(tangency condition), and

$F(u)=F_{k}(u,\cdots,u)+O(\|u\|^{k+1}_{\alpha}),$ (2.2)

with $p>k\geq 2$, and $F_{k}$ denoting the leading-order operator (on
$D(A^{\alpha})$) in the Taylor expansion of $F$, near the origin. A broad
class of stochastic equations that arise in many branches of physics can be
recasted into this abstract functional framework. Such equations driven by
(possibly spatially inhomogeneous) multiplicative noise are met in various
contexts such as in turbulence theory (intermittency phenomena) [Bir13], ocean
dynamics [SFL01], climate dynamics [KJSK08, MFC09, CSG11, AN19, BAB+17], non-
equilibrium phase transitions [CH93, Hin00, Muñ04, SSGO07], statistical
physics [SB79] or population dynamics [HSZ02, CPT16].

For later usage, we denote by $V$, the function space $D(A^{\alpha})$. In
particular $D(A)\subset V\subset H$. The norm of $H$ is denoted by
$\|\cdot\|$, while the norm of $V$ is denoted by $\|\cdot\|_{V}.$ Under mild
conditions often met in applications about the spectrum of $L_{\lambda}$ (in
particular near the instability onset, [CLW15a, Secns. 3.3 and 6.2]), one can
then prove that for any initial condition $u_{0}$ in $V$ (sufficiently small)
333The existence and uniqueness result proved in [CLW15a] for globally
Lipschitz nonlinearity, applies to the case of a standard cutoff argument of
the nonlinearity ensures existence and uniqueness for sufficiently small
initial data., emanates a (unique) solution to Eq. (2.1) in the sense that
$v=e^{-Z(t,\omega)}u$ is the unique classical solution of

$\frac{\mathrm{d}v}{\mathrm{d}t}=L_{\lambda}v+Z(t,\omega)v+e^{-Z(t,\omega)}F(e^{Z(t,\omega)}v),\;v(0)=e^{-Z(0,\omega)}u_{0},$
(2.3)

that lies, $\mathbb{P}$-a.s., in

$\displaystyle C((0,\tau);D(A))\cap C([0,\tau);V)\cap C^{1}((0,\tau);H),$
(2.4)

where $\tau$ is a stopping time444Defined, loosely speaking, as the (random)
minimum time after which $v(t,\omega)$ exits the neighborhood of the origin
over which the required estimates apply. and $Z(t,\omega)$ denotes the
Ornstein-Uhlenbeck (OU) process, stationary solution of the scalar Langevin
equation

$\,\mathrm{d}r=-r\,\mathrm{d}t+\sigma\,\mathrm{d}W_{t}(\omega).$ (2.5)

Furthermore, $u$ depends continuously on $\lambda$ and on $u_{0}$, and is a
stochastic $H$-valued adapted process.

Note that to derive (2.3) from (2.1), we need to note that, by using the Itô’s
formula,

$\displaystyle\,\mathrm{d}e^{-Z(t,\omega)}$
$\displaystyle=(Z(t,\omega)e^{-Z(t,\omega)}+\frac{\sigma^{2}}{2}e^{-Z(t,\omega)})\,\mathrm{d}t-\sigma
e^{-Z(t,\omega)}\,\mathrm{d}W_{t}$
$\displaystyle=Z(t,\omega)e^{-Z(t,\omega)}\,\mathrm{d}t-\sigma
e^{-Z(t,\omega)}\circ\,\mathrm{d}W_{t},$

where the last equality follows from the conversion between Itô and
Stratonovich integrals. On the other hand,

$\,\mathrm{d}v=\,\mathrm{d}(e^{-Z(t,\omega)}u)=u\circ\,\mathrm{d}e^{-Z(t,\omega)}+e^{-Z(t,\omega)}\circ\,\mathrm{d}u,$

and (2.3) results from (2.1), after simplification.

Note that the usage of the OU process $Z(t,\omega)$ through the transformation
$v=e^{-Z(t,\omega)}u$ is here motivated by the fact that the PDE with random
coefficients (2.3) enjoys a better temporal regularity than the SPDE (2.1) due
to the Hölder regularity of the OU process [CLW15a, Lemma 3.1], which allows
in turn for applying standard existence and uniqueness results from the theory
of non-autonomous PDEs; see [Hen81, Theorem 3.3.3 and Corollary 3.5.3]. We
refer to [CLW15a, Proposition 3.1] and [CLW15a, Appendix A], for more details.

In this article, we are concerned with the study of bifurcations for the class
of problems that can be recast into the abstract formulation of Eq. (2.1). In
other words, we are interested in describing how linear instability translates
to nonlinear dynamics, in presence of multiplicative noise. We are in
particular concerned with the dynamical reduction problem near the onset of
instability, i.e. to describe by a low-dimensional object the SPDE dynamics in
case of e.g. the existence of a critical value $\lambda_{c}$ at which $m$
eigenvalues $\beta_{j}(\lambda)$ of $L_{\lambda}$ (counting algebraic
multiplicity) cross the imaginary axis. In this case, it is known, under the
assumptions recalled above, that a finite-dimensional $C^{p}$-smooth
stochastic (local) invariant manifold, $\mathfrak{M}_{\lambda}(\omega)$,
exists near the origin [CLW15a, Corollary 5.1 and Prop. 6.1]. This manifold is
obtained as the graph of a random mapping, $h^{\lambda}$ over the space
$H_{\mathfrak{c}}$ of eigenmodes $\bm{e}_{j}$ losing stability. Inspired by
the deterministic theory, such manifolds are natural objects to consider for
reduction of SPDEs near the instability onset. In [CLW15a, CLW15b], the
approximation problem of such manifolds was thus investigated. Theorem 2.1
below gives a summary of the main approximation results derived in [CLW15a,
Theorem 6.1] (see also [CLW15a, Corollary 6.1]) which allows, in practice, for
deriving explicit reduced systems.

However, due to noise, intrinsic difficulties arise and the prediction of
transitions for Stochastic PDEs, based on such reduced systems, requires a
particular care. Indeed, what makes the success of center/unstable manifold
theory to describe bifurcations for PDEs from low-dimensional reduced systems
fails to operate in the stochastic realm, calling for a new approach to
apprehend bifurcation analysis in presence of noise. The reason is that a
breakdown of the reduction principle to the center/unstable manifold now takes
place, and exponential slaving of the stable modes onto the unstable ones is
no longer valid, due to large excursions caused by the (white) noise, even
when its intensity is small. We come back to this issue in Sec. 2.3 below and
the later sections (Secns. 3 and 4) about the error estimates. Next we recall
the approximation formulas that play a key role in the derivation of these
error estimates.

###  2.2. Leading-order approximation of stochastic invariant manifolds, and
non-Markovian terms

Note that $L_{\lambda}$ has a compact resolvent by recalling that $D(A)$ is
compactly and densely embedded in $H$ [EN00, Prop. II.4.25]. As a consequence,
since $L_{\lambda}\colon D(A)\rightarrow H$ is a closed operator due to the
sectorial property of $-L_{\lambda}$, we have that for each $\lambda$, the
spectrum of $L_{\lambda}$, $\sigma(L_{\lambda})$, consists only of isolated
eigenvalues with finite algebraic multiplicities; see [Kat95, Thm. III-6.29]
(see also [EN00, Corollary IV.1.19]). This spectral property combined with the
sectorial property of $-L_{\lambda}$ implies that there are at most finitely
many eigenvalues with a given real part. The sectorial property of
$-L_{\lambda}$ also implies that the real part of the spectrum, ${\rm
Re\,}\sigma(L_{\lambda})$, is bounded above; see also [EN00, Thm. II.4.18].

These two properties about ${\rm Re\,}\sigma(L_{\lambda})$ allow us in turn to
label elements in $\sigma(L_{\lambda})$ according to the lexicographical
order. According to this rearrangement, we can label the eigenvalues by a
single positive integer $n$, so that

$\displaystyle\sigma(L_{\lambda})=\\{\beta_{n}(\lambda)\>|\>n\in\mathbb{N}^{\ast}\\},$
(2.6)

with, for any $1\leq n<n^{\prime}$, either

$\displaystyle{\rm Re\,}\beta_{n}(\lambda)>{\rm
Re\,}\beta_{n^{\prime}}(\lambda),$ (2.7)

or

$\displaystyle{\rm Re\,}\beta_{n}(\lambda)={\rm
Re\,}\beta_{n^{\prime}}(\lambda),\;\text{ and }\;{\rm
Im\,}\beta_{n}(\lambda)\geq{\rm Im\,}\beta_{n^{\prime}}(\lambda).$ (2.8)

In this convention, an eigenvalue of algebraic multiplicity $m$, is repeated
$m$ times. Hereafter, this rearrangement is mainly used for simplifying the
notations of some theoretical developments, whereas the labeling with
wavenumbers is often restored when dealing with applications; see Sec. 5
below.

As already mentioned, we are concerned with describing how linear
instabilities translate to the nonlinear dynamics, in presence of
multiplicative noise. To do so, the onset of instability is described in terms
of the principle of exchange of stabilities (PES) [MW05], concerned with the
loss of stability of the basic steady state. More precisely, the PES describes
situations for which the spectrum of $L_{\lambda}$ experiences the following
change at a critical parameter $\lambda_{c}$:

$\displaystyle{\rm Re\,}\beta_{j}(\lambda)\begin{cases}<0&\mbox{if
}\lambda<\lambda_{c},\\\ =0&\mbox{if }\lambda=\lambda_{c},\\\ >0&\mbox{if
}\lambda>\lambda_{c},\end{cases}$ $\displaystyle 1\leq j\leq m,$ (2.9)
$\displaystyle{\rm Re\,}\beta_{j}(\lambda_{c})<0,$ $\displaystyle j\geq m+1,$

for some $m>0$, and for $\lambda$ in some neighborhood $\Lambda$ of
$\lambda_{c}$. To this PES condition (2.9) is associated the following
decomposition of $\sigma(L_{\lambda})$:

$\displaystyle\sigma_{\mathfrak{c}}(L_{\lambda})=\\{\beta_{j}(\lambda)\>|\>j=1,\
2,\ \cdots,m\\},$ (2.10)
$\displaystyle\sigma_{\mathfrak{s}}(L_{\lambda})=\\{\beta_{j}(\lambda)\>|\>j=m+1,\
m+2,\ \cdots\\}.$

The PES condition prevents eigenvalues from
$\sigma_{\mathfrak{s}}(L_{\lambda})$ to cross the imaginary axis as $\lambda$
varies in $\Lambda$. Hence, no eigenvalues other than those of
$\sigma_{\mathfrak{c}}(L_{\lambda})$ change sign in $\Lambda$. Furthermore,
the PES condition implies the following uniform spectral gap by reducing
$\Lambda$ accordingly [CLW15a, Lemma 6.1],

$0>2k\eta_{\mathfrak{c}}>\eta_{\mathfrak{s}},$ (2.11)

where $k$ is the leading order of $F$ (see (2.2)), and

$\displaystyle\eta_{\mathfrak{c}}=\inf_{\lambda\in\Lambda}\inf_{j=1,\cdots,m}\\{\mathrm{Re}(\beta_{j}(\lambda))\\},$
$\displaystyle\eta_{\mathfrak{s}}=\sup_{\lambda\in\Lambda}\sup_{j\geq
m+1}\\{\mathrm{Re}(\beta_{j}(\lambda))\\}.$

Such a uniform spectral gap condition is required for the approximation
formulas of stochastic invariant manifolds recalled in Theorem 2.1 below; see
[CLW15a, Theorem 6.1] and [CLW15a, Corollary 6.1]. It also implies
(exponential) dichotomy estimates (see [CLW15a, Eqns. (3.24a,b,c)] ) to be
satisfied by the semigroup generated by $L_{\lambda}$ and that are key to
ensure the existence of stochastic invariant manifolds (see [CLW15a, Appendix
B]) associated with the $m$ modes losing their stability according to (2.9).

To these modes, we associate the reduced state space, $H_{\mathfrak{c}}$,
given by

$H_{\mathfrak{c}}=\mathrm{span}\\{\bm{e}_{1},\cdots,\bm{e}_{m}\\},$ (2.12)

while a mode $\bm{e}_{n}$ with $n\geq m+1$ denotes a stable mode. In what
follows, the indices $j_{1},\cdots,j_{k}$ correspond to the $m$ critical
wavenumbers losing stability (allowing repetition) at $\lambda=\lambda_{c}$.
We finally introduce the notation

$(\bm{k},\bm{\beta}_{c}(\lambda))=\sum_{\ell=1}^{k}\beta_{j_{\ell}}(\lambda),$
(2.13)

with $\bm{k}$ denoting the $k$-tuple, $(1,\cdots,1)$. The projector onto the
subspace $H_{\mathfrak{s}}$ (resp. $H_{\mathfrak{c}}$) spanned by the stable
modes (resp. given by (2.12)) is denoted by $\Pi_{\mathfrak{s}}$ (resp.
$\Pi_{\mathfrak{c}}$). The inner product in $H$ is denoted by
$\langle\cdot,\cdot\rangle$. We have then the following stochastic invariant
manifold approximation theorem which lies at the core of our reduction
approach near the onset of instability for non-equilibrium systems subject to
random fluctuations.

######  Theorem 2.1.

Assume that $F$ and $L_{\lambda}$ satisfy the assumptions recalled above, and
that the PES condition (2.9) is satisfied. Then for each $\lambda$ in a
neighborhood $\Lambda$ of $\lambda_{c}$, Eq. (2.1) admits a stochastic (local)
invariant manifold,
$\mathfrak{M}_{\omega}^{\lambda}=\mbox{graph}(h_{\omega}^{\lambda})$, with
$h_{\omega}^{\lambda}$ that maps $H_{\mathfrak{c}}$ into the regular space
$D(A^{\alpha})$, for every $\omega$ in $\Omega$.

The mapping $h_{\omega}^{\lambda}$ characterizing
$\mathfrak{M}_{\omega}^{\lambda}$ is a time-dependent stochastic mapping which
is approximated by another such mapping $\Phi_{\omega}^{\lambda}$ satisfying
for any $\epsilon>0$,

$\|h_{\omega}^{\lambda}(X,t)-\Phi_{\omega}^{\lambda}(X,t)\|_{V}\leq\epsilon\|X\|^{k},\;\omega\in\Omega,\;t\in\mathbb{R},$
(2.14)

as long as $X$ lies in a neighborhood $\mathcal{N}_{\epsilon}$ of the origin
in $H_{\mathfrak{c}}$ spanned by the $m$ eigenmodes $\bm{e}_{j}$ losing their
stability (see (2.12)), as lambda crosses $\lambda_{c}.$

The mapping $\Phi_{\omega}^{\lambda}$ in (2.14) is given explicitly by the
following random Lyapunov-Perron integral:

$\displaystyle\Phi_{\omega}^{\lambda}(X,t)=\int_{-\infty}^{0}e^{\sigma(k-1)W_{t,s}(\omega){\rm
Id}}e^{-sL_{\lambda}}P_{\mathfrak{s}}F_{k}(e^{sL_{\lambda}}X)\,\mathrm{d}s,$
(2.15)

with $W_{t,s}(\omega)=W_{t+s}(\omega)-W_{t}(\omega).$

This mapping is well defined if
$(\bm{k},\bm{\beta}_{c}(\lambda))>\beta_{n}(\lambda)$ for $n\geq m+1$ and
possesses the following expansion:

$\Phi^{\lambda}_{\omega}(X,t)=\sum_{n\geq
m+1}\Phi^{n,\lambda}_{\omega}(X,t)\bm{e}_{n},\;X\in
H_{\mathfrak{c}},\;\omega\in\Omega.$ (2.16)

Here

$\Phi^{n,\lambda}_{\omega}(X,t)=\sum_{(j_{1},\cdots,j_{k})}F_{j_{1}\cdots
j_{k}}^{n}M_{j_{1}\cdots j_{k}}^{n,\lambda}(t,\omega)X_{j_{1}}\cdots
X_{j_{k}},$ (2.17)

where $1\leq j_{1},\cdots,j_{k}\leq m$, $X_{j}=\langle X,e_{j}\rangle$, and
the $F_{j_{1}\cdots j_{k}}^{n}$ are coefficients accounting for the $n$-th
component of the nonlinear interactions (through the leading-order term
$F_{k}$) between the low modes $\bm{e}_{j_{1}}$, $\cdots$, $\bm{e}_{j_{k}}$
(in $H_{c}$), namely:

$F_{j_{1}\cdots j_{k}}^{n}=\langle
F_{k}(\bm{e}_{j_{1}},\cdots,\bm{e}_{j_{k}}),\bm{e}_{n}^{\ast}\rangle,\;1\leq
j_{1},\cdots,j_{k}\leq m,$ (2.18)

where the $\bm{e}_{n}^{\ast}$ denote the eigenmodes of the adjoint operator of
$L_{\lambda}$.

The $M_{j_{1}\cdots j_{k}}^{n,\lambda}$-terms are path-dependent coefficients
making explicit the $(t,\omega)$-dependence of $\Phi_{\omega}^{\lambda}$,
which are obtained as stationary solutions of the auxiliary scalar SDE:

$\mathrm{d}M=\left(1-\Big{(}(\bm{k},\bm{\beta}_{c}(\lambda))-\beta_{n}(\lambda)\Big{)}M\right)\mathrm{d}t-\sigma(k-1)M\circ\mathrm{d}W_{t}.$
(2.19)

This theorem is a reformulation of [CLW15a, Proposition 6.1 and Theorem 6.1]
whose proofs can be found in [CLW15a, Section 6.4]. Note that although proved
for $L_{\lambda}$ self-adjoint, the conclusions of these results extend to the
case $L_{\lambda}$ diagonalizable in $\mathbb{C}$. We emphasize that such
approximation results have been also extended in [CLW15a, Corollary 7.1] to
situations in which stable modes are included in the reduced state space which
has important consequences for certain applications requiring higher-
dimensional reduced state spaces to ensure better approximation of the SPDE
dynamics as previously illustrated in [CLW15b, Chapter 6] in the context of
capturing by reduced equations noise-induced large excursions.

The above theorem provides via the coefficients of
$\Phi_{\omega}^{\lambda}(X,t)$ an explicit tracking of the time-stochastic
dependence of the (unknown) random mapping $h_{\omega}^{\lambda}(X,t)$. Prior
works proposed deterministic approximation formulas for
$h_{\omega}^{\lambda}(X,t)$ in the case $F(u)=B(u,u)$ [BW10], or $F(u)=u^{p}$
(both with error $O(\|X\|)$) [CDZ11]; see also [CLW15b, Section 4.1]. However,
such deterministic approximation formulas are of limited virtues for
applications; see e.g. example (LABEL:AB_syst).

Instead, the path-dependent coefficients $M_{j_{1}\cdots j_{k}}^{n,\lambda}$
(simply referred to as $M$-terms hereafter), depend on the past of the noise,
when, for $\lambda$ in a neighborhood of $\lambda_{c}$, the following non-
resonance conditions are satisfied

$\Bigl{(}F_{j_{1}\cdots j_{k}}^{n}\neq 0\Bigr{)}\Longrightarrow\biggl{(}{\rm
Re\,}(\bm{k},\bm{\beta}_{c}(\lambda))>{\rm
Re\,}\beta_{n}(\lambda)\biggr{)},\;1\leq j_{1},\cdots,j_{k}\leq m,\;n\geq
m+1.$ (NR)

Under such an (NR)-condition, the $M$-terms are indeed given by

$M_{j_{1}\cdots
j_{k}}^{n,\lambda}(t,\omega)=\int_{-\infty}^{0}e^{\big{(}(\bm{k},\bm{\beta}_{c}(\lambda))-\beta_{n}(\lambda)\big{)}s+\sigma(k-1)W_{t,s}(\omega)}\mathrm{d}s,$
(2.20)

with $W_{t,s}(\omega)=W_{t+s}(\omega)-W_{t}(\omega)$. Note that due to the NR-
condition, the integral in (2.20) is well defined almost surely, since the
Brownian motion $W_{t}$ satisfies the sublinear growth condition
$\lim_{t\rightarrow\pm\infty}W_{t}/t=0,$ almost surely (see [CLW15a, Lemma
3.1]). Finally, note that the PES condition implies that the NR-condition is
satisfied for $\lambda$ in neighborhood of $\lambda_{c}$, as
$\sum_{\ell=1}^{k}{\rm
Re\,}\beta_{j_{\ell}}(\lambda)>k\eta_{\mathfrak{c}}>2k\eta_{\mathfrak{c}}>\eta_{\mathfrak{s}}\geq{\rm
Re\,}\beta_{n}(\lambda)$ due to (2.11).

The exponential decaying integrand in (2.20) depends on the distance,
$(\bm{k},\bm{\beta}_{c}(\lambda))-\beta_{n}(\lambda)$, between the
$\overbrace{(1:\cdots:1)}^{(k)}$-resonance made up from critical eigenvalues,
and the stable eigenvalue $\beta_{n}(\lambda)$. Such a distance to resonance
controls in turns the decay of temporal correlations of the $M$-terms,
provided that the noise intensity, $\sigma$, lies in some admissible range;
see Lemma A.1 in Appendix A. As a result, these terms exhibit decay of
correlations of a “reddish” nature, although allowing for non-Gaussian
statistics with heavy tails all the more pronounced than
$\sigma/\sqrt{(\bm{k},\bm{\beta}_{c}(\lambda))-\beta_{n}(\lambda)}$ is close
to $1$ from below (in the case of real eigenvalues); see Appendix A.2.

The $M$-coefficients are thus non-Markovian and convey exogenous memory
effects in the sense of [HO07, Hai09]. It is worthwhile noting that the “noise
bath” is essential for these coefficients to exhibit decay of correlations,
the $M$-terms being reduced to simply the constant term
$((\bm{k},\bm{\beta}_{c}(\lambda))-\beta_{n}(\lambda))^{-1}$, when $\sigma=0$,
recovering standard approximations formulas of (deterministic) invariant
manifolds; see [CLM20, Theorem 2].

These non-Markovian terms are new ingredients produced by the interactions
between the noise with the nonlinear effects. The reduced systems built from
the corresponding non-Markovian parameterization,
$\Phi^{\lambda}_{\omega}(X,t)$ given by Theorem 2.1, are thus, in general,
non-Markovian SDEs. In the case where $F(u)$ is a quadratic nonlinearity,
$B(u,u)$, such as arising in fluid problems (see Secns. 3 and 4 below), such a
non-Markovian reduced SDE takes the form

$\displaystyle\mathrm{d}X=\bigg{(}\Pi_{\mathfrak{c}}L_{\lambda}X+\Pi_{\mathfrak{c}}\Big{(}B(X,X)+$
$\displaystyle B(\Phi(X,t),X)+B(X,\Phi(X,t))$ (2.21)
$\displaystyle+B(\Phi(X,t),\Phi(X,t))\Big{)}\bigg{)}\mathrm{d}t+\sigma
X\circ\,\mathrm{d}W_{t},\;X\in\mathbb{R}^{m},$

which in coordinate forms reads

$\displaystyle\mathrm{d}X_{\ell}$
$\displaystyle=\bigg{(}\beta_{\ell}(\lambda)X_{\ell}+\overbrace{\sum_{i,j=1}^{m}B^{\ell}_{ij}X_{i}X_{j}}^{(a)}+\overbrace{\sum_{p,i,j=1}^{m}\sum_{n=m+1}^{\infty}(B^{\ell}_{pn}+B^{\ell}_{np})B^{n}_{ij}M_{ij}^{n,\lambda}(t,\omega)X_{i}X_{j}X_{p}}^{(b)}$
(2.22)
$\displaystyle+\underbrace{\sum_{n,n^{\prime}=m+1}^{\infty}\sum_{i,j=1}^{m}\sum_{p,q=1}^{m}B^{\ell}_{nn^{\prime}}B_{ij}^{n}B_{pq}^{n^{\prime}}M_{ij}^{n,\lambda}(t,\omega)M_{pq}^{n^{\prime},\lambda}(t,\omega)X_{i}X_{j}X_{p}X_{q}}_{(c)}\bigg{)}\,\mathrm{d}t$
$\displaystyle\hskip 85.35826pt+\sigma X_{\ell}\circ\mathrm{d}W_{t},\quad\
1\leq\ell\leq m,\;\;\text{with $B^{\ell}_{\nu\mu}=\langle
B(\bm{e}_{\nu},\bm{e}_{\mu}),\bm{e}_{\ell}^{\ast}\rangle$}.$

Here the $X_{\ell}$ are aimed at approximating the amplitudes of the SPDE
solution onto the $m$ modes $\bm{e}_{\ell}$ ($1\leq\ell\leq m$) losing
stability, as $\lambda$ crosses $\lambda_{c}$.

In the expansion (2.22), the $(a)$-terms account for self-interactions between
the unstable modes, the $(b)$-terms account for cross-interactions between
unstable and stable modes, and the $(c)$-terms for self-interactions between
the stable modes. The non-Markovian terms appear only out of the two latter
types of interactions.

As we will see, the $M$-terms are key, under certain circumstances, to
appropriately track the fluctuations and large-excursions caused by the
(white) noise, although the original SPDE is Markovian. This is particularly
true for SPDEs forced by a multiplicative noise in the orthogonal subspace of
e.g. its mildly unstable mode. Our parameterization formulas show that the
noise gets transmitted to this mode via non-Markovian coefficients, and that
the reduced equation is only stochastically driven by the latter. This
situation is a special case of more general spatially inhomogeneous
multiplicative noises, for which the remark below points out the required
modifications of Theorem 2.1.

######  Remark 2.1.

Theorem 2.1 above can be extended to the case of SPDEs driven by spatially
inhomogeneous multiplicative noise, of the form

$\mathrm{d}u=\big{(}L_{\lambda}u+F(u)\big{)}\mathrm{d}t+\sigma(u)\circ\mathrm{d}W_{t},$
(2.23)

where $L_{\lambda}$ and $F$ satisfy the assumptions of Theorem 2.1 and

$\sigma(u)=\sum_{j=1}^{N_{f}}\sigma_{j}\langle
u,\bm{e}_{j}^{\ast}\rangle\bm{e}_{j},$ (2.24)

with $\sigma_{j}\geq 0$ and $N_{f}\geq 1$. While the proof of [CLW15a, Theorem
6.1] requires a certain care to handle such multiplicative noises, the same
approach applies and leads to approximation formulas of stochastic invariant
manifolds for Eq. (2.23) which are still given by (2.16)–(2.17), with the only
difference that the $M$-terms are now replaced by

$\displaystyle M_{j_{1}\cdots
j_{k}}^{n,\lambda}(t,\omega)=\int_{-\infty}^{0}e^{\big{(}(\bm{k},\bm{\beta}_{c}(\lambda))-\beta_{n}(\lambda)\big{)}s+(\bm{k},\bm{\sigma})W_{t,s}(\omega)}\mathrm{d}s,$
(2.25)

where

$(\bm{k},\bm{\sigma})=\sum_{\ell=1}^{k}\sigma_{j_{\ell}}-\sigma_{n}.$ (2.26)

These $M$-terms are obtained as stationary solutions of the auxiliary scalar
SDE:

$\displaystyle\mathrm{d}M=\left(1-\Big{(}(\bm{k},\bm{\beta}_{c}(\lambda))-\beta_{n}(\lambda)\Big{)}M\right)\mathrm{d}t-(\bm{k},\bm{\sigma})M\circ\mathrm{d}W_{t}.$
(2.27)

######  Remark 2.2.

It is also interesting to point out that when Theorem 2.1 is applied to
(finite-dimensional) SDE systems of the form
$\mathrm{d}y=\big{(}Ay+F(y)\big{)}\mathrm{d}t+\sigma y\circ\mathrm{d}W_{t}$
with $A$ being a $d\times d$ matrix and
$F:\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}$ a sufficiently smooth function
satisfying $F(0)=0$ and $DF(0)=0$, our approximation formula (2.15) coincides
with the one determined by the cohomological equation given by [Arn98, Eq.
(8.4.9)] for the lowest-order non-trivial $H^{s}_{p0}$ therein. Indeed, for
such SDEs, after performing the change of variables $x=e^{-Z(t,\omega)}y$, we
obtain the following RDE system

$\frac{\,\mathrm{d}x}{\mathrm{d}t}=Ax+Z(t,\omega)x+e^{-Z(t,\omega)}F(e^{Z(t,\omega)}x);$
(2.28)

see again (2.3). Denoting the analogue of $\Phi$ for this transformed RDE by
$\Psi$ and adopting the invariance equation approach such as reviewed in
[CLM20, Sec. 2.2], the homological equation satisfied by $\Psi(t,\omega)$
reads

$\frac{\,\mathrm{d}\Psi}{\,\mathrm{d}t}+D_{x_{\mathfrak{c}}}\Psi(x_{\mathfrak{c}},t)(A_{\mathfrak{c}}x_{\mathfrak{c}}+Z(t,\omega)x_{\mathfrak{c}})-(A_{\mathfrak{s}}\Psi+Z(t,\omega)\Psi)=\Pi_{\mathfrak{c}}e^{-Z(t,\omega)}F_{k}(e^{Z(t,\omega)}x_{\mathfrak{c}}),$
(2.29)

where $k\geq 2$ denotes the lowest-order term in the Taylor expansion of $F$.
Note that compared with [CLM20, Eq. (2.27)] for the deterministic case, the
above equation involves also a time-derivative term
$\frac{\,\mathrm{d}\Psi}{\,\mathrm{d}t}$ arising from the time-dependence (or
noise-path dependence) nature of $\Psi$.

Note that (2.29) is the same as [Arn98, Eq. (8.4.9)] when the RDE system
(2.28) is considered and when $p$ therein equals the lowest order $k$ in the
Taylor expansion of $F$ here, since at this lowest order the term
$R_{p0}^{c}(x_{c})$ in [Arn98, Eq. (8.4.9)] simply equals
$\Pi_{\mathfrak{c}}e^{-Z(t,\omega)}F_{k}(e^{Z(t,\omega)}x_{\mathfrak{c}})$.

One one hand, the approach adopted in [Arn98] for finite-dimensional SDEs aims
to handle more general SDE settings when e.g. the basic reference state is
more complicated than a steady state, by relying on the multiplicative ergodic
theorem and the associated Lyapunov spectra and Oseledets subspaces to perform
the decomposition of the stochastic flow. Extension to SPDEs using the
Oseledets ergodic theorem has been pursued in [MZZ08].

On the other hand, the practical aspects in the general setting are still too
intricate to derive reduced systems with easily computable coefficients even
for high-dimensional SDEs (let alone SPDEs). Indeed, as pointed out in [Arn98,
Sec. 8.4.3] already for the application to a two-dimensional SDE system
associated with the Duffing-van der Pol oscillator in the Pitchfork
bifurcation scenario, “the computational effort for these results is enormous
and could only be accomplished by using the computer algebra program MAPLE ….”

Instead, by adopting a standard normal modes framework, although we lose the
characterization of perturbations near complicated states, we gain insights
about (i) the derivation of explicit reduced equations based on rigorous
approximation formulas (as reviewed in Theorem 2.1), and (ii) error estimates
(see Theorem 4.1 below) for SPDEs that cover the challenging case of fluid
problems.

Figure 1. Schematic of the stochastic invariant manifold approximation
summarized by Theorem 2.1. The time-dependent stochastic manifold,
$\Phi^{\lambda}_{\omega}(X,t)$ (red curve) approximates the stochastic
invariant manifold $h^{\lambda}_{\omega}(X,t)$ (blue curve). For a given
noise’s realization $\omega$, the shape of the latter changes in the course of
time. The approximation $\Phi^{\lambda}_{\omega}(X,t)$ is able to track these
changes due to the path-dependent, non-Markovian, coefficients,
$M_{j_{1}\cdots j_{k}}^{n,\lambda}(t,\omega)$, solving the SDE (2.19). In this
schematic, the black curve shows a solution path to the SPDE that evolves on
the stochastic invariant manifold.

###  2.3. Tracking the large-excursions through the non-Markovian $M$-terms

As already mentioned, although Theorem 2.1 provides an analytic approximation
of the (unknown) slaving function $h_{\omega}^{\lambda}$, and in particular
regarding its dependence with respect to time and the noise path, it suffers
from a technical restriction, namely that the approximation (2.14) is valid
over a deterministic neighborhood $\mathcal{N}_{\epsilon}$ of the origin.

From an analysis viewpoint, this theorem concedes thus some taste of
dissatisfaction since as soon as $X$ becomes a time-dependent random
variable555Such as when $X$ is a solution to non-Markovian reduced equation
(2.21)., large excursions caused by the (white) noise are expected to take
place even for small noise intensity, pushing $X$ outside of
$\mathcal{N}_{\epsilon}$ in the course of time and questioning de facto the
validity of the approximation (2.14). This restriction is technical (from the
proof), but even if the neighborhood would be allowed to fluctuate with time,
still noise could drive the solutions outside of it.

Similarly, exponential attraction results of invariant manifolds that occur in
presence of a sufficiently large spectral gap and that extend to SPDEs with
globally Lipschitz nonlinearity ([CLW15a, Corollary 4.3]), encounter however
similar technical restrictions for SPDEs with nonlinearities that are only
locally Lipschitz, such as for fluid problems. The current proofs do not allow
indeed for establishing a pathwise exponential attraction of local stochastic
invariant manifolds, preventing to rely on a reduction principle as in the
deterministic theory [Cra91], to analyze bifurcations. This important but
subtle point, misled some authors in their conclusions regarding the analogue
of a reduction principle for SPDEs; see [CLW15a, Remark 4.3]. We provide in
Secns. 3 and 4 below error estimates accounting for large-excursions and that
are furthermore relevant for stochastic fluid problems, allowing in turn to
make precise the use of stochastic invariant manifolds in bifurcation analysis
of such SPDEs.

Before embarking in these error estimates, we provide below a simple example
to illustrate these issues, and to show that the non-Markovian $M$-terms allow
for tracking the large excursions caused by the noise. For this purpose, we
consider the stochastic system:

$\displaystyle\,\mathrm{d}A=(\lambda A+AB)\,\mathrm{d}t,\;\lambda>0$ (2.30)
$\displaystyle\,\mathrm{d}B=(\kappa B-A^{2})\,\mathrm{d}t+\sigma
B\circ\,\mathrm{d}W_{t},\;\kappa<0,\sigma>0,$

in which the noise forces only the equation of the stable mode, here the
$B$-equation. This system arises in diverse fields, with or without stochastic
forces. For instance it arises in the study of nonlinear crystals [SB79], and
more recently in the study of the finite time blow-up problem for the 3-D
Navier Stokes equations [Tao16], where such systems are used as elemental
“pump gates” to execute transition of energy from one mode to another.

Here, we are looking for a closed equation describing the evolution of $A$
without resolving $B$. If one adopts ideas of deterministic unstable/center-
manifold reduction to seek for such a closure [KW83], we arrive using this
theory at the parameterization $B=-(2\lambda-\kappa)^{-1}A^{2}$, which leads
to the closure

$\dot{X}=\lambda X-(2\lambda-\kappa)^{-1}X^{3}.$ (2.31)

This closure is here only deterministic and thus cannot capture the stochastic
nature of the dynamics of (LABEL:AB_syst).

We would like instead to use an explicit stochastic parameterization of $B$.
To be of practical interest, this stochastic parameterization should be able
to track the large excursions caused by the white noise. Application of our
parameterization formulas to this particular system (see Remark 2.1), provides
such a stochastic parameterization. The latter is given here by the non-
Markovian parameterization $B=-M_{11}^{2,\lambda}(t,\omega)A^{2}$, with
$M_{11}^{2,\lambda}$ being the stationary solution of (see (2.27))

$\,\mathrm{d}M=(1-(2\lambda-\kappa)M)\,\mathrm{d}t+\sigma
M\circ\,\mathrm{d}W_{t},$ (2.32)

namely

$M_{11}^{2,\lambda}(t,\omega)=\int_{-\infty}^{0}e^{(2\lambda-\kappa)s-\sigma(W_{t+s}(\omega)-W_{s}(\omega))}\,\mathrm{d}s.$
(2.33)

The closure becomes then the following ODE with a path-dependent coefficient:

$\dot{X}=\lambda X-M_{11}^{2,\lambda}(t,\omega)X^{3}.$ (2.34)

Recall that the AB-system is forced by a noise in the orthogonal subspace of
its mildly unstable mode. In this case, our parameterization formulas show
that the noise gets transmitted to this mode via the non-Markovian coefficient
$M_{11}^{2,\lambda}$, and that the reduced equation is only stochastically
driven by the latter.

The practical efficiency of such parameterizations to track the large
excursions induced by the white noise is illustrated in Fig. 2. The left panel
shows a solution path of (LABEL:AB_syst) displayed in the (A,B)-plane (black
dash curve), and two time instances of the non-Markovian unstable manifold
approximation, $B=-M_{11}^{2,\lambda}(t,\omega)A^{2}$ (red curves). In
contrast this figures shows that the deterministic unstable manifold (blue
curve) gives only a poor approximation of the solution path’s average motion.
The consequence is that, given a noise realization $\omega$, the $A$\- and
$B$-trajectories as simulated from (LABEL:AB_syst) (black curves), are
approximated to a good precision by the non-Markovian closure (2.34) (red
curves), whereas the closure (2.31) based on the deterministic unstable
manifold is only proposing a deterministic constant approximation.

Figure 2. Left panel: A solution path of (LABEL:AB_syst) shown in the
(A,B)-plane (black dash curve), and two time instances of the non-Markovian
unstable manifold approximation, $B=-M_{11}^{2,\lambda}(t,\omega)A^{2}$ (red
curves). The latter is able to track the large excursions of the solution
path, while the deterministic unstable manifold gives only a poor
approximation to its average motion. Right, top-to-bottom panels: $A$\- and
$B$-time series as simulated from (LABEL:AB_syst) (black curves), the
deterministic closure (2.31) (blue curves), and non-Markovian closure (2.34)
(red curves).

The reasons for such a successful closure by Eq. (2.34), can be worked it out
on this example as system (LABEL:AB_syst) possesses a natural closure which
allows for providing a more intuitive “raison d’être” of the non-Markovian
parameterization $B=-M_{11}^{2,\lambda}(t,\omega)A^{2}$. Indeed, by
integrating the $B$-equation of system (LABEL:AB_syst), a closed form equation
can be naturally derived for describing the evolution of $B$, given by the
following integro-differential equation with path-dependent kernel

$\frac{\,\mathrm{d}A}{\,\mathrm{d}t}=\lambda
A-\bigg{(}\overbrace{\int_{-\infty}^{t}e^{\kappa(t-s)+\sigma(W_{t}(\omega)-W_{s}(\omega))}A^{2}(s)\,\mathrm{d}s}^{I(t)}\bigg{)}A.$
(2.35)

Thus the quality of approximation of the integral term in (2.35) by the non-
Markovian parameterization $B=-M_{11}^{2,\lambda}(t,\omega)A^{2}$ is
conditioning the efficiency of the closure (2.34).

It is worthwhile noting that this non-Markovian parameterization can be
obtained by integration of a backward-forward system [CLW15b, Chap. 4], which
for this example, consists of approximating $A(s)$ (for $s\leq t$) — in the
integral term of (2.35) and only for this term — by the solution of the
following backward linear equation

$\displaystyle\frac{\,\mathrm{d}A}{\,\mathrm{d}s}^{(1)}=\lambda
A^{(1)}(s),\;s<t$ (2.36) $\displaystyle A^{(1)}(t)=A(t).$

Thus, by approximating $A(s)$ in the integral term of (2.35), with

$A(s)=e^{\lambda(s-t)}A(t)+o(A(t)),$ (2.37)

we obtain the following equation (up to order 3),

$\frac{\,\mathrm{d}A}{\,\mathrm{d}t}=\lambda
A-\bigg{(}\overbrace{\int_{-\infty}^{t}e^{(\kappa-2\lambda)(t-s)+\sigma(W_{t}(\omega)-W_{s}(\omega))}\,\mathrm{d}s}^{J(t)}\bigg{)}A^{3},$

from which we recover the closure (2.34) (see also (2.33)), by using the
change of variable $s^{\prime}=s-t$ in the integral term. The reader may
wonder why approximating the exact closure (2.35) by (2.34), on this
particular example. The reasons are multiple, especially in terms of insights
gleaned by this operation. First, it shows the relevance of the non-Markovian
parameterization resulting from application of Theorem 2.1, for certain
parameter regimes like that corresponding of Fig. 1 (for which $J(t)\approx
I(t)$). Second, it shows the relevance — for stochastic systems — of
approximating a fully non-Markovian closure like (2.35) by a non-Markovian
closure (2.34) that involves only the past of the noise, the latter allowing
for an efficient simulation of the path-dependent coefficient
$M_{11}^{2,\lambda}(t,\omega)$ in (2.33) by an SDE decoupled from the
$A$-equation (i.e. by simulating (2.32)). Instead, the closure (2.35) would
require to resolve the corresponding time integral $I(t)$ at each time step to
compute $A(t)$, an operation that is costly numerically. This decoupling from
the $A$-equation to simulate the $M$-term and thus $J(t)$ approximating $I(t)$
may be seen as a small advantage here due to the small size of system
(LABEL:AB_syst), but such a feature is highly beneficial numerically when the
size of the original system gets much larger, as for fluid models; see Sec. 5
below.

As much as striking might be the ansatz (LABEL:Eq_ansatz), the non-Markovian
closure (2.34) depending only on the past of the noise, provides thus, under
certain circumstances, a good approximation of $A(t)$ obtained from the closed
form equation (2.35) depending furthermore on the past-history of $A(t)$. The
general error estimates of Secns. 3 and 4 below allow us to clarify these
circumstances, for this particular example and beyond. In particular, they
show that when $\lambda$ is sufficiently small and the noise’s intensity
$\sigma$, scales as $\sqrt{\lambda}$, the non-Markovian closure (2.34)
(depending on the past of the noise) provides a good approximation of the
exact non-Markovian closure (2.35) (depending not only on the past of the
noise but also the past of $A$).

##  3\. Modeling error estimates from reduced equations

In this section we apply the approximation formulas of Theorem 2.1 to derive
general error estimates, establishing further legitimacy of the reduced
equations that can be derived by application of Theorem 2.1, near the onset of
instability.

To simplify the presentation of the ideas, we focus on the following class of
SPDE

$\mathrm{d}u=\big{(}L_{\lambda}u+B(u,u)\big{)}\mathrm{d}t+\sigma
u\circ\,\mathrm{d}W_{t},$ (3.1)

and place ourselves in a pitchfork bifurcation scenario when $\sigma=0$.

###  3.1. Assumptions and pitchfork scenario

Recall that, $L_{\lambda}$ denotes a linear operator on the Hilbert space $H$
such that $L_{\lambda}=-A+P_{\lambda}$, where $A$ is the generator of a
contraction semigroup in $H$, while $P_{\lambda}$ is a low-order perturbation
of $A$, i.e. $P_{\lambda}:{\color[rgb]{0,0,0}V}\rightarrow H$ is a family of
bounded linear operators depending continuously on the parameter $\lambda$,
with $D(A)\subset V\subset H$. Recall also that $D(A)$ is compactly and
densely embedded in $H$.

Denote the eigenelements of $L_{\lambda}$ by
$\\{(\beta_{k}(\lambda),\bm{e}_{k})\;:\;k\in\mathbb{N}\\}$. To derive our
error estimates, we assume

######  Condition 1.

The linear operator $L_{\lambda}$ is self-adjoint.

######  Condition 2.

The principle of exchange of stabilities (PES) condition is satisfied with
$m=1$ in (LABEL:splitting), namely there exists $\lambda_{c}$ such that:

$\displaystyle\beta_{1}(\lambda)\begin{cases}<0&\mbox{if
}\lambda<\lambda_{c},\\\ =0&\mbox{if }\lambda=\lambda_{c},\\\ >0&\mbox{if
}\lambda>\lambda_{c},\end{cases}$ (3.2)
$\displaystyle\beta_{j}(\lambda_{c})<0,\quad\textrm{for all }j\geq 2,$

and all the eigenvalues $\beta_{j}(\lambda)$, for $j\geq 2$, remain negative
for $\lambda$ in some interval $[\lambda_{c},\lambda^{*}]$, with
$\lambda^{*}>\lambda_{c}$.

######  Condition 3.

The nonlinear term $B$ is a continuous bilinear map from $V\times V$ to $H$.
Furthermore the following conditions hold:

$\langle B(u,u),u\rangle=0,\;\textrm{ for any }u\textrm{ in }V,$ (3.3)

and there exists $N\geq 1$ for which

$\langle B(\bm{e}_{1},\bm{e}_{1}),\bm{e}_{n}\rangle=0,\;n>N,$ (3.4)

while for at least one $2\leq n\leq N$, $\langle
B(\bm{e}_{1},\bm{e}_{1}),\bm{e}_{n}\rangle\neq 0.$

Note that if $\langle B(\bm{e}_{1},\bm{e}_{1}),\bm{e}_{n}\rangle=0$ for all
$n\geq 2$, then the manifold $\Phi_{\lambda}$ will be identically zero and one
needs then to consider higher-order approximation formulas; see [CLW15b].

###  3.2. Preparatory estimates for the pitchfork scenario

The parameterization of the modes $\bm{e}_{2}$ to $\bm{e}_{N}$ is given by

$\Phi(X,t,\omega)=\sum_{n=2}^{N}\Phi_{n}(X,t,\omega)\bm{e}_{n},$ (3.5)

with

$\Phi_{n}(X,t,\omega)=B_{11}^{n}M_{11}^{n,\lambda}(t,\omega)X^{2},$ (3.6)

where

$B_{ij}^{n}=\langle B(\bm{e}_{i},\bm{e}_{j}),\bm{e}_{n}\rangle.$ (3.7)

Here $X$ denotes the amplitude of $\bm{e}_{1}$. In what follows we will omit
the dependence on $\omega$, reserving this explicit dependence only when
necessary.

Due to (3.3), we have $\langle B(\bm{e}_{1},\bm{e}_{1}),\bm{e}_{1}\rangle=0$
and the invariant manifold (IM) reduced equation, up to order $O(X^{3})$,
takes the following form

$\mathrm{d}X=\Big{(}\epsilon
X+X\mathcal{F}_{\mathfrak{c}}(\Phi(X,t))\Big{)}\mathrm{d}t+\sigma
X\circ\,\mathrm{d}W_{t},$ (3.8)

with

$\mathcal{F}_{\mathfrak{c}}(v)=\left\langle\big{(}B(\bm{e}_{1},v)+B(v,\bm{e}_{1})\big{)},\bm{e}_{1}\right\rangle,\quad\forall
v\in V_{\mathfrak{s}}.$ (3.9)

Note that

$\displaystyle X\mathcal{F}_{\mathfrak{c}}(\Phi(X,t))$
$\displaystyle=X\bigg{(}\sum_{n=2}^{N}\Phi_{n}(X,t)\big{(}B_{1n}^{1}+B_{n1}^{1}\big{)}\bigg{)}$
(3.10)
$\displaystyle=\bigg{(}\sum_{n=2}^{N}B_{11}^{n}M_{11}^{n,\lambda}(t)\big{(}B_{1n}^{1}+B_{n1}^{1}\big{)}\bigg{)}X^{3}.$

We assume then

######  Condition 1.

The coefficients in (3.10), satisfy for $\lambda$ in
$[\lambda_{c},\lambda^{*}]$, the dissipation condition

$\sum_{n=2}^{N}B_{11}^{n}M_{11}^{n,\lambda}(t)\big{(}B_{1n}^{1}+B_{n1}^{1}\big{)}<0,\;\text{
for all }t.$ (3.11)

Condition 1 is a sufficient condition for the reduced equation (3.8) to
experience a stochastic supercritical pitchfork bifurcation as $\epsilon$
changes sign; see Sec. 5.4 below and Appendix D.

Since $M_{11}^{n,\lambda}$ is always positive, (3.11) can be ensured if we
have for all $2\leq n\leq N$,

$B_{11}^{n}\big{(}B_{1n}^{1}+B_{n1}^{1}\big{)}\leq 0,$ (3.12)

and at least one of the inequalities in (3.12) is a strict one. In case
$\sigma=0$, such a condition is a sufficient condition for the deterministic
reduced equation Eq. (3.8) to experience a pitchfork bifurcation [Cra91,
MW05].

We turn now to the derivation of an equation satisfied by $\Phi$ which is key
in our estimates presented below. From (3.6), we obtain by using the
Stratonovich form of Itô’s formula [Kun90]:

$\,\mathrm{d}\Phi_{n}=B_{11}^{n}X^{2}\circ\,\mathrm{d}M_{11}^{n,\lambda}+2B_{11}^{n}M_{11}^{n,\lambda}X\circ\,\mathrm{d}X.$
(3.13)

Since here $(\bm{k},\beta_{c}(\lambda))=2\epsilon$ and $k=2$ in (2.19), we
have that $M_{11}^{n,\lambda}$ satisfies

$\mathrm{d}M=\left(1-\big{(}2\epsilon-\beta_{n}(\lambda)\big{)}M\right)\mathrm{d}t-\sigma
M\circ\mathrm{d}W_{t}.$ (3.14)

Using (3.14) and (3.8) in (3.13), we get after simplification that

$\displaystyle\,\mathrm{d}\Phi_{n}=\big{(}\beta_{n}(\lambda)\Phi_{n}+B_{11}^{n}X^{2}$
$\displaystyle+2\Phi_{n}\mathcal{F}_{c}(\Phi(X,t))\big{)}\,\mathrm{d}t+\sigma\Phi_{n}\circ\,\mathrm{d}W_{t}.$
(3.15)

Then, by (3.5) and (3.15), the equation satisfied by $\Phi$ is:

$\displaystyle\,\mathrm{d}\Phi=\Big{(}L^{\mathfrak{s}}_{\lambda}\Phi+X^{2}B_{11}^{\mathfrak{s}}+2\Phi\mathcal{F}_{\mathfrak{c}}($
$\displaystyle\Phi(X,t))\Big{)}\,\mathrm{d}t+\sigma\Phi\circ\,\mathrm{d}W_{t},$
(3.16)

where $L^{\mathfrak{s}}_{\lambda}=\Pi_{\mathfrak{s}}L_{\lambda}$ and we have
used

$\sum_{n=2}^{N}B_{11}^{n}X^{2}\bm{e}_{n}=X^{2}\Pi_{\mathfrak{s}}B(\bm{e}_{1},\bm{e}_{1}),$
(3.17)

which holds thanks to (3.4), and we have introduced

$B_{11}^{\mathfrak{s}}=\Pi_{\mathfrak{s}}B(\bm{e}_{1},\bm{e}_{1}).$ (3.18)

To summarize, we get from (3.8) and (3.16) that $(X,\Phi)$ satisfy the
following system:

$\displaystyle\mathrm{d}X$ $\displaystyle=\Big{(}\epsilon
X+X\mathcal{F}_{\mathfrak{c}}(\Phi(X,t))\Big{)}\mathrm{d}t+\sigma
X\circ\,\mathrm{d}W_{t},$ (3.19) $\displaystyle\,\mathrm{d}\Phi$
$\displaystyle=\Big{(}L^{\mathfrak{s}}_{\lambda}\Phi+X^{2}B_{11}^{\mathfrak{s}}+2\Phi\mathcal{F}_{\mathfrak{c}}(\Phi(X,t))\Big{)}\,\mathrm{d}t+\sigma\Phi\circ\,\mathrm{d}W_{t}.$

We have then the following lemma, whose proof is provided in Appendix B.

######  Lemma 3.1 (A priori estimates for $X$ and $\Phi$).

Assume that Condition 1 holds. Consider the system (3.19) for which the
initial condition $(X(0),\Phi(0,\omega))$ satisfies, for any $\omega$ in
$\Omega$,

$\displaystyle X(0)$ $\displaystyle=X_{0},$ (3.20)
$\displaystyle\Phi(0,\omega)$
$\displaystyle=\sum_{n=2}^{N}\Big{(}B_{11}^{n}M_{11}^{n,\lambda}(0,\omega)X_{0}^{2}\Big{)}\bm{e}_{n},$

where $N$ is the integer that appears in (3.4).

Assume furthermore that Conditions 1, 2, and 3 hold, with
$\beta_{1}(\lambda)=\epsilon$ in 2. Assume finally that
$\sigma=\sqrt{\epsilon}$ and that

$|X_{0}|\sim\sqrt{\epsilon}.$ (3.21)

Then, for any given $T>0$ and any $\chi$ in $(0,1)$, there exists a constant
$C>0$ depending on $\chi$ but independent of $\epsilon$, for which the
following estimate holds:

$\mathbb{P}\left(\sup_{[0,T/\epsilon]}|X(t)|\leq
C\sqrt{\epsilon}\;\;\textrm{and}\sup_{[0,T/\epsilon]}\|\Phi(t)\|\leq
C\,\epsilon\right)\geq 1-\chi.$ (3.22)

######  Remark 3.1.

  * (i)

The proof of Lemma 3.1 given in Appendix B relies on Lemmas A.2 and A.3,
themselves proved in their respective appendices. These latter lemmas show
that the constant $C$ in (3.22) grows as $\chi$ decreases, i.e. the larger the
probability with which (3.22) holds, the closer one needs to push $\epsilon$
towards $0$ in order for the corresponding a priori bounds about $|X(t)|$ and
$\|\Phi(t)\|$ to be below a target value.

  * (ii)

Note that a solution $(X(t),\Phi(t))$ to the system (3.19) emanating from an
initial condition satisfying (3.20), implies that $\Phi(X,t)$ is given by the
parameterization (3.5), for all $t>0$. As it can be observed in the proof of
Lemma 3.1, it is the dynamic interpretation of $\Phi(X,t)$ (i.e. the equation
it solves in system (3.19)) that allows for deriving estimates about
$\Phi(X,t)$, while $X$ is a solution to the reduced equation (3.8).

###  3.3. Residual error estimates in large probability

We first rewrite the original SPDE (3.1) into a coupled system for
$u_{\mathfrak{c}}(t)=\Pi_{\mathfrak{c}}u(t)=x(t)\bm{e}_{1}$ and
$u_{\mathfrak{s}}=\Pi_{\mathfrak{s}}u$:

$\displaystyle\mathrm{d}x=\big{(}\epsilon
x+x\mathcal{F}_{\mathfrak{c}}(u_{\mathfrak{s}})+\Pi_{1}B(u_{\mathfrak{s}},u_{\mathfrak{s}})\big{)}\mathrm{d}t+\sigma
x\\!\circ\\!\,\mathrm{d}W_{t},$ (3.23a)
$\displaystyle\mathrm{d}u_{\mathfrak{s}}\\!=\\!\big{(}L_{\lambda}^{\mathfrak{s}}u_{\mathfrak{s}}+x^{2}B_{11}^{\mathfrak{s}}+x\mathcal{F}_{\mathfrak{s}}(u_{\mathfrak{s}})+\Pi_{\mathfrak{s}}B(u_{\mathfrak{s}},u_{\mathfrak{s}})\big{)}\mathrm{d}t+\sigma
u_{\mathfrak{s}}\circ\,\mathrm{d}W_{t},$ (3.23b)

where $F_{\mathfrak{c}}$ is defined in (3.9) and

$\displaystyle\Pi_{1}B(u_{\mathfrak{s}},u_{\mathfrak{s}})=\langle
B(u_{\mathfrak{s}},u_{\mathfrak{s}}),\bm{e}_{1}\rangle,$ (3.24)
$\displaystyle\mathcal{F}_{\mathfrak{s}}(v)=\Pi_{\mathfrak{s}}(B(\bm{e}_{1},v)+B(v,\bm{e}_{1})),\quad\forall
v\in V_{\mathfrak{s}}.$

Note that the term $\Pi_{1}\big{(}x^{2}B(\bm{e}_{1},\bm{e}_{1})\big{)}$ does
not appear in the equation for $x$ since it equals $x^{2}\langle
B(\bm{e}_{1},\bm{e}_{1}),\bm{e}_{1}\rangle$, which is zero thanks to (3.3).

For any arbitrary stochastic processes $X(t)\bm{e}_{1}$ and $Y(t)$ that are
adapted to the underlying stochastic basis and evolve respectively in
$H_{\mathfrak{c}}$ and in $V_{\mathfrak{s}}$, we define the residual
$\mathcal{R}(t)=(\mathcal{R}_{1}(t),\mathcal{R}_{2}(t))$, where

$\displaystyle\mathcal{R}_{1}(t)$ $\displaystyle=X(t)-e^{\epsilon t+\sigma
W_{t}}X(0)-\int_{0}^{t}e^{\epsilon(t-s)+\sigma(W_{t}-W_{s})}E_{1}(s)\,\mathrm{d}s,$
(3.25) $\displaystyle\mathcal{R}_{2}(t)$
$\displaystyle=Y(t)-e^{tL_{\lambda}^{\mathfrak{s}}+\sigma
W_{t}}Y(0)-\int_{0}^{t}e^{L_{\lambda}^{\mathfrak{s}}(t-s)+\sigma(W_{t}-W_{s})}E_{2}(s)\,\mathrm{d}s,\;\text{with}$
$\displaystyle E_{1}(s)$
$\displaystyle=X(s)\mathcal{F}_{\mathfrak{c}}(Y(s))+\Pi_{1}B(Y(s),Y(s)),$
$\displaystyle E_{2}(s)$
$\displaystyle=X^{2}(s)B_{11}^{\mathfrak{s}}+X(s)\mathcal{F}_{\mathfrak{s}}(Y(s))+\Pi_{\mathfrak{s}}B(Y(s),Y(s)).$

This residual measures the modeling error made by the process $(X(t),Y(t))$ in
satisfying the system (3.23), and thus the modeling error made by
$X(t)\bm{e}_{1}+Y(t)$ in satisfying the SPDE (3.1).

Now pick up any solution $(X(t),\Phi(t))$ of the surrogate system (3.19), for
which the initial condition satisfies the relation (3.6) at $t=0$. Note that
such a solution is well defined for all time $t\geq 0$, since under this
constraint on the initial data, the system (3.19) is equivalent to the reduced
equation (3.8).

We have the following bounds for the residual:

######  Theorem 3.1.

Consider the solution $(X(t),\Phi(t))$ to Eq. (3.19) and its residual error
defined in (3.25). Assume that the assumptions of Lemma 3.1 are satisfied.
Then, for any $T>0$ and any $\chi$ in $(0,1)$, there exists a constant $C>0$
independent of $\epsilon$ for which

$\mathbb{P}\left(\sup_{t\in[0,T/\epsilon]}|\mathcal{R}_{1}(t)|\leq
C\epsilon\right)\geq 1-\chi,$ (3.26)

and

$\mathbb{P}\left(\sup_{t\in[0,T/\epsilon]}\|\mathcal{R}_{2}(t)\|\leq
C\epsilon^{3/2}\right)\geq 1-\chi.$ (3.27)

###### Proof.

Let $(X(t),\Phi(t))$ be a solution to (3.19). After integration and
simplifications, we obtain that (3.25) becomes

$\displaystyle\mathcal{R}_{1}(t)=-\int_{0}^{t}e^{\epsilon(t-s)+\sigma(W_{t}-W_{s})}\Pi_{1}B(\Phi(s),\Phi(s))\,\mathrm{d}s,$
(3.28)
$\displaystyle\mathcal{R}_{2}(t)=\int_{0}^{t}e^{L_{\lambda}^{\mathfrak{s}}(t-s)+\sigma(W_{t}-W_{s})}E_{3}(s)\,\mathrm{d}s,$

with

$E_{3}=2\Phi\mathcal{F}_{\mathfrak{c}}(\Phi)-X\mathcal{F}_{\mathfrak{s}}(\Phi)-\Pi_{s}B(\Phi,\Phi).$
(3.29)

Then, the desired results follow directly from the a priori estimates (3.22)
together with the probabilistic estimate about the Brownian motion given in
Lemma A.3. Indeed, by letting $\Omega^{*}$ and $\gamma$ be the same as given
in the proof of Lemma 3.1 (see (B.19)), we have proved that
$\|\Phi(t,\omega)\|\leq C\epsilon$ for $t$ in $[0,T/\epsilon]$ and $\omega$ in
$\Omega^{*}$; cf. (B.23). Since $\Phi(t,\omega)$ as defined in (3.5) takes
value in a finite dimensional subspace of $V$ due to the assumption (3.4), the
two norms $\|\cdot\|_{V}$ and $\|\cdot\|$ are equivalent when restricted to
this subspace. We get then $\|\Phi(t,\omega)\|_{V}\leq C\epsilon$ by
redefining $C$. Thus,

$\displaystyle|\mathcal{R}_{1}(t,\omega)|$ $\displaystyle\leq
C_{B}\int_{0}^{t}e^{\epsilon(t-s)+2\gamma}\|\Phi(s,\omega)\|_{V}^{2}\,\mathrm{d}s$
(3.30) $\displaystyle\leq
C_{B}C^{2}e^{2\gamma}\int_{0}^{t}e^{\epsilon(t-s)}\,\mathrm{d}s\,\epsilon^{2}$
$\displaystyle=C_{B}C^{2}e^{2\gamma}(e^{\epsilon t}-1)\epsilon$
$\displaystyle\leq C_{B}C^{2}e^{2\gamma}(e^{T}-1)\epsilon,\quad\forall
t\in[0,T/\epsilon],\omega\in\Omega^{*}.$

The estimate (3.26) follows.

Finally, recall that for such $t$ and $\omega$, we have that
$|X(t,\omega)|\leq C\sqrt{\epsilon}$ for some $C>0$ independent of $\epsilon$;
see (B.20). Using the above estimates for $X$ and $\Phi$, we get then

$\displaystyle\|\mathcal{R}_{2}(t)\|$
$\displaystyle\leq\int_{0}^{t}e^{\beta_{2}(\lambda)(t-s)+2\gamma}E_{3}(s)\,\mathrm{d}s$
(3.31)
$\displaystyle\leq\widetilde{C}(\epsilon^{3/2}+\epsilon^{2}),\quad\forall
t\in[0,T/\epsilon],\omega\in\Omega^{*},$

where $\widetilde{C}>0$ is another constant independent of $\epsilon$. We have
thus derived the estimate (3.27). ∎

######  Remark 3.2.

When more general spatially inhomogeneous multiplicative noise such as given
in Remark 2.1 is considered, we can still derive residual estimates in the
style of (3.26) and (3.27) under the same conditions as required in Theorem
3.1 with $\sigma_{n}=c_{n}\sqrt{\epsilon}$ for some $c_{n}\geq 0$,
$n=1,\ldots,N$. The residual $\mathcal{R}_{2}$ needs to be adapted accordingly
by projecting onto each stable mode $n=2,\cdots,N$ as well.

The derivation of the associated residual estimates follows the same lines of
arguments presented in the proof of Theorem 3.1 once the a priori estimate for
$X(t)$ and $\Phi(t)$ given by (3.22) is established. The derivation of this
latter a priori estimate follows the same steps as presented in Appendix B,
but by relying on the following transformation rather than the one given in
(B.1):

$U=e^{-2\epsilon
t-2\sigma_{1}W_{t}}X^{2},\quad\Psi_{n}=e^{-2\beta_{n}t-2\sigma_{n}W_{t}}\Phi_{n}^{2}.$
(3.32)

##  4\. Low- and high-mode error estimates

In this section, we go beyond the residual error estimates presented above to
derive error estimates between the solutions of the finite-dimensional
surrogate system (3.19) and those of the original SPDE (3.1). To do so, it
requires a priori bounds for solutions to the SPDE (3.1). The main difficulty
compared to those conducted for solutions of (3.19) (cf. Lemma 3.1) lies in
the infinite-dimensional nature of $u_{\mathfrak{s}}$, whereas $\Phi$ is of
finite-dimensional range due to (3.4). As a consequence, the loss of
regularity via the nonlinearity $B$ is now effective which requires new
estimates.

To cope with these difficulties one needs extra assumptions on the linear part
$L_{\lambda}$. To the previous assumptions on this operator, we require in
particular an extra condition to be satisfied (Condition 2) for the linear
dissipative effects to suitably counterbalance the loss of regularity due to
the nonlinear terms in $B$. These conditions are described below.

######  Condition 1.

Let $L_{\lambda}=-A+P_{\lambda}$ satisfy the Conditions 1 and 2 of Sec. 3,
with $A$ given by

$A=-\Delta\mbox{diag}(\bm{\nu})$ (4.1)

where $\bm{\nu}$ is a $d$-dimensional vector with positive entries, denoted by
$\nu_{j}$, and $\mbox{diag}(\bm{\nu})$ denotes the $d\times d$ diagonal matrix
with $\bm{\nu}$ as the diagonal entries. We assume furthermore that the
boundary conditions and functional setting are such that for any $u$ in
$D(A)$, $\langle Au,u\rangle\geq 0$ and $\sqrt{\langle Au,u\rangle}$ defines
(after performing integration by parts) a norm $\|u\|_{\bm{\nu}}$ on $V$
equivalent to $\|u\|_{V}$.

######  Condition 2.

There exist $\eta$ in $(0,1)$ and $\delta>0$ such that for all $\lambda$ in
the interval $[\lambda_{c},\lambda^{*}]$ given in Condition 2, the following
damping relations hold for all $u_{\mathfrak{s}}$ in
$\mathcal{W}_{\mathfrak{s}}=D(A)\cap H_{\mathfrak{s}}$:

$\displaystyle\langle-\eta
Au_{\mathfrak{s}}+P_{\lambda}u_{\mathfrak{s}},u_{\mathfrak{s}}\rangle\leq-\delta\|u_{\mathfrak{s}}\|^{2},$
(4.2a) $\displaystyle\langle-\eta
Au_{\mathfrak{s}}+P_{\lambda}u_{\mathfrak{s}},Au_{\mathfrak{s}}\rangle\leq-\delta\|u_{\mathfrak{s}}\|^{2}_{\bm{\nu}}.$
(4.2b)

For this class of linear operators, the error estimates derived hereafter
allow us then, from the reduced equations based on Theorem 2.1, to conclude to
a pitchfork bifurcation scenario in large probability. Our results apply to a
broad class of SPDEs such as given by (3.1) (allowing even for inhomogeneous
noise), as long as the noise intensity $\sigma$ and the eigenvalue’s magnitude
of the mildly unstable mode, scale accordingly, while the linear and nonlinear
terms obey standard (energy-preserving) assumptions encountered in fluid
problems (Condition 3). Applications to Rayleigh-Bénard convection are
detailed in Sec. 5 below.

######  Remark 4.1.

The form of $A$ assumed in Condition 1 is chosen to simplify the derivation of
the error estimates in Theorem 4.1. More general differential operators $A$
than (4.1) could have been considered for which the conclusions of Theorem 4.1
would still hold. In that respect, an analogue of Conditions 1 and 2 would
still hold for a differential operator $A$ in divergence form, that (possibly)
includes higher-order derivatives and that satisfies a strong ellipticity
condition (see [Tem97] and [Paz83, Sec. 7.2]). Given a positive integer $m$,
such an operator writes for any $u$ in the Sobolev space
$(H^{2m}(\mathcal{D}))^{d}$ (over a smooth bounded domain $\mathcal{D}$)

$Au=\sum_{\begin{subarray}{c}\alpha,\beta\in\mathbb{N}^{p}\\\
[\alpha]=[\beta]=m\end{subarray}}(-1)^{m}D^{\beta}(a_{\alpha\beta}(x)D^{\alpha}u),$
(4.3)

where $\alpha=(\alpha_{1},\ldots,\alpha_{p})\in\mathbb{N}^{p}$,
$[\alpha]=\alpha_{1}+\cdots+\alpha_{p}$, the coefficients $a_{\alpha\beta}$
are matrix-valued mappings which are sufficiently smooth into the space of
$d\times d$ symmetric matrices, and $a_{\alpha\beta}=a_{\beta\alpha}$. If the
lower bound of $\langle Au,u\rangle$ is controlled as follows (for any $u$ in
$D(A)$)

$\langle
Au,u\rangle\geq\kappa\sum_{\alpha\in\mathbb{N}^{p},[\alpha]=m}\|D^{\alpha}u\|^{2},\;\text{
for some }\;\kappa>0,$ (4.4)

then with
$\|u\|_{V}=\sum_{\alpha\in\mathbb{N}^{p},[\alpha]=m}\|D^{\alpha}u\|^{2}$,
analogues of Conditions 1 and 2 could be formulated for which the conclusions
of Theorem 4.1 would still hold. In Sec. 5 below, it is shown that Conditions
1 and 2 as formulated here, i.e. for the operator $A$ defined by (4.1), and
with the linear damping estimates (4.2a)-(4.2b), are sufficient to apply
Theorem 4.1 to a standard Rayleigh-Bénard convection problem subject to
multiplicative noise.

###  4.1. A priori bounds on SPDE solutions

Roughly speaking, we consider the following a priori bounds on the SDPE
solution $(x(t),u_{\mathfrak{s}}(t))$

$\displaystyle\text{(I)}:\;\sup_{t\in[0,T/\epsilon]}|x(t)|\leq
C\sqrt{\epsilon},$ (4.5)
$\displaystyle\textrm{(II)}:\;\sup_{t\in[0,T/\epsilon]}\|u_{\mathfrak{s}}(t)\|\leq
C\epsilon,$
$\displaystyle\textrm{(III)}:\;\sup_{t\in[0,T/\epsilon]}\|u_{\mathfrak{s}}(t)\|_{V}\leq
C\epsilon,$

and are concerned whether these estimates hold in large probability. These a
priori bounds are key to derive our main Theorem about error estimates; see
Theorem 4.1 and its proof. We recall that the SPDE solutions are here
understood in the sense recalled in Sec. 2.1. In that respect, to derive a
priori bounds like (4.5), one first work with the transformed equation Eq.
(2.3), followed by the inverse transformation to infer back the desired bounds
about the SPDE solutions. We refer to Appendix C for more details and provide
here the main elements and ideas to derive (I), (II), and (III) in large
probability.

In that perspective, observe that standard energy estimates on the solutions
to (3.23) lead to

$x^{2}(t,\omega)+\|u_{\mathfrak{s}}(t,\omega)\|^{2}\leq e^{2\epsilon t+2\sigma
W_{t}(\omega)}\|u(0,\omega)\|^{2}.$ (4.6)

Thus, if one assumes $|u(0)|\sim\sqrt{\epsilon}$ one would obtain that, in
large probability, $|x(t)|\leq C\sqrt{\epsilon}$ and
$\|u_{\mathfrak{s}}(t)\|\leq C\sqrt{\epsilon}$ by conducting similar estimates
as in Step 1 of the proof of Lemma 3.1 (and using (3.3) for controlling the
nonlinear term); see Appendix B.

To derive error estimates between $u_{\mathfrak{s}}(t)$ and its
parameterization $\Phi(t)$, this is however insufficient, as the bound on
$\|u_{\mathfrak{s}}(t)\|$ should scale as $\epsilon$ and not $\sqrt{\epsilon}$
since $\|\Phi(t)\|$ scales like $\epsilon$ due (3.22). The main difficulty
compared to Step 2 of the proof of Lemma 3.1, lies in the infinite-dimensional
nature of $u_{\mathfrak{s}}$, whereas $\Phi$ is of finite-dimensional range
due to (3.4). As a consequence, the loss of regularity via the nonlinearity
$B$ is now effective which requires new estimates.

For instance, when estimating $\|u_{\mathfrak{s}}\|$, by taking the inner
product of the $u_{\mathfrak{s}}$-equation with $u_{\mathfrak{s}}$ in (3.23),
this loss of regularity causes for instance the appearance in the RHS of terms
depending on $\|u_{\mathfrak{s}}\|_{V}$ (and not only on
$\|u_{\mathfrak{s}}\|$) since e.g. $\langle
B(\bm{e}_{1},u_{\mathfrak{s}}),u_{\mathfrak{s}}\rangle\leq
C_{B}\|\bm{e}_{1}\|_{V}\|u_{\mathfrak{s}}\|_{V}\|u_{\mathfrak{s}}\|$.

Nevertheless, such difficulties can be handled through a careful exploitation
of the damping effects brought by the linear part allowing for bypassing to
estimate $\|u_{\mathfrak{s}}\|_{V}$ (at this stage), and leading us to prove
that, in large probability, and for $t$ in $[0,T/\epsilon]$,

$\|u_{\mathfrak{s}}(t)\|^{2}\leq
e^{2\gamma}\|u_{\mathfrak{s}}(0)\|^{2}+\frac{Ce^{4\gamma}}{\delta}\epsilon^{2},$
(4.7)

where $\gamma$ is as given in Lemma A.3, $\delta$ is given in (4.2), and $C>0$
is a generic constant independent of $\epsilon$. The proof of this key
estimate is given in Appendix C.1.

As a consequence, if $|x(0)|\sim\sqrt{\epsilon}$ and
$\|u_{\mathfrak{s}}(0)\|\sim\epsilon,$ one obtains that for any $T>0$ and any
$\chi$ in $(0,1)$,

$\mathbb{P}\Big{(}\textrm{(I) and (II)}\Big{)}\geq 1-\chi.$ (4.8)

Nevertheless, since $B:V\times V\rightarrow H$ (loss of regularity), to derive
a full set of a priori estimates about the SPDE solution, one cannot escape
for estimating $\|u_{\mathfrak{s}}\|_{V}$. To do so, we take the inner product
of the $u_{\mathfrak{s}}$-equation with $Au_{\mathfrak{s}}$ in (3.23), and
after conducting the relevant estimates (see Appendix C.2), we arrive at, in
large probability, and for $t$ in $[0,T/\epsilon]$,

$\|u_{\mathfrak{s}}(t)\|_{V}^{2}\leq
C\|u_{\mathfrak{s}}(0)\|_{V}^{2}+C\epsilon^{2}\int_{0}^{t}e^{-\delta(t-s)+2\sqrt{\epsilon}(W_{t}-W_{s})}\,\mathrm{d}s,$
(4.9)

namely $\|u_{\mathfrak{s}}(t)\|_{V}^{2}\leq C\epsilon^{2},$ provided that
$\|u_{\mathfrak{s}}(0,\omega)\|_{V}\sim\epsilon$ and
$0\leq\epsilon\leq\epsilon^{*}$ for some $\epsilon^{*}>0$ sufficiently small.

Thus, if

$|x(0)|\sim\sqrt{\epsilon},\;\;\|u_{\mathfrak{s}}(0)\|\sim\epsilon,\;\text{
and }\|u_{\mathfrak{s}}(0)\|_{V}\sim\epsilon.$ (4.10)

we have that for any $T>0$ and any $\chi$ in $(0,1)$,

$\mathbb{P}\Big{(}\textrm{(I), (II), and (III)}\Big{)}\geq 1-\chi.$ (4.11)

###  4.2. Fluid problems subject to fluctuations: Low- and high-mode error
estimates

From these probabilistic estimates we prove next our main error estimate
theorem. This theorem gives the error estimates made on the low- and high-mode
dynamics, when they are respectively approximated by the solution $X(t)$ to
the reduced equation Eq. (3.8), and the parameterization $\Phi(X(t),t)$ given
by (3.5)-(3.6).

######  Theorem 4.1.

Assume that Conditions 1, 2, and 3 are satisfied. Denote
$\epsilon=\beta_{1}(\lambda)>0$. Let $u(t)=x(t)\bm{e}_{1}+u_{\mathfrak{s}}(t)$
be a solution to the SPDE (3.1).

Assume furthermore that Condition 1 is satisfied, that
$\sigma=\sqrt{\epsilon}$ and that

$\displaystyle X(0)=x(0)\sim\sqrt{\epsilon},$ (4.12)
$\displaystyle\|u_{\mathfrak{s}}(0)\|\sim\epsilon,\;\|u_{\mathfrak{s}}(0)\|_{V}\sim\epsilon,$
$\displaystyle\;\|u_{\mathfrak{s}}(0)-\Phi(X,0)\|\sim\epsilon^{3/2}.$

Consider the errors

$\displaystyle a(t)=x(t)-X(t),\;\textrm{(low-mode error)}$ (4.13)
$\displaystyle b(t)=u_{\mathfrak{s}}(t)-\Phi(X,t)\;\textrm{(high-mode
error)},$

in which $X(t)$ solves the reduced equation Eq. (3.8), and the
parameterization $\Phi$ is given by (3.5)-(3.6).

Then, for any $T>0$ and any $\chi$ in $(0,1)$, there exist $\epsilon^{*}>0$
and $C>0$ such that for any $\epsilon$ in $[0,\epsilon^{*}]$, the following
error estimate holds:

$\mathbb{P}\Big{(}\sup_{[0,T/\epsilon]}|a(t)|\leq
C\epsilon\;\text{and}\;\sup_{[0,T/\epsilon]}\|b(t)\|\leq
C\epsilon^{3/2}\Big{)}\geq 1-\chi.$ (4.14)

###### Proof.

First, let us write down the system of equations satisfied by the low-mode and
high-mode errors, $a$ and $b$. By recalling that $(X,\Phi)$ and
$(x,u_{\mathfrak{s}})$ satisfy (3.19) and (3.23), respectively, we infer that
$(a,b)$ satisfies

$\displaystyle\mathrm{d}a$ $\displaystyle=\Big{(}\epsilon
a+x\mathcal{F}_{\mathfrak{c}}(u_{\mathfrak{s}})+\Pi_{1}B(u_{\mathfrak{s}},u_{\mathfrak{s}})-X\mathcal{F}_{\mathfrak{c}}(\Phi(X,t))\Big{)}\mathrm{d}t+\sigma
a\circ\,\mathrm{d}W_{t},$ (4.15) $\displaystyle\mathrm{d}b$
$\displaystyle=\Big{(}L_{\lambda}^{\mathfrak{s}}b+x^{2}B_{11}^{\mathfrak{s}}+x\mathcal{F}_{\mathfrak{s}}(u_{\mathfrak{s}})+\Pi_{\mathfrak{s}}B(u_{\mathfrak{s}},u_{\mathfrak{s}})-X^{2}B_{11}^{\mathfrak{s}}-2\Phi\mathcal{F}_{\mathfrak{c}}(\Phi(X,t))\Big{)}\mathrm{d}t+\sigma
b\circ\,\mathrm{d}W_{t},$

with $\mathcal{F}_{\mathfrak{c}}$ given by (3.9), $B_{11}^{\mathfrak{s}}$ by
(3.18), and $\Pi_{1}B$ and $\mathcal{F}_{\mathfrak{s}}$ defined in
(LABEL:Eq_Fs).

We detail below the main elements to derive (4.14) from (4.15). We start with
the estimates of the nonlinear terms controlling the amplitudes of $a(t)$ and
$b(t)$, in their respective dynamical equations. For the $a$-equation, this
term is given by $\Pi_{1}B(u_{s},u_{\mathfrak{s}})+I(t)$ with

$I(t)=x\mathcal{F}_{\mathfrak{c}}(u_{\mathfrak{s}})-X\mathcal{F}_{\mathfrak{c}}(\Phi(X,t)),$
(4.16)

while for the $b$-equation, it is given by
$J(t)+(x^{2}-X^{2})B_{11}^{\mathfrak{s}}$ with

$J(t)=x\mathcal{F}_{\mathfrak{s}}(u_{\mathfrak{s}})+\Pi_{\mathfrak{s}}B(u_{\mathfrak{s}},u_{\mathfrak{s}})-2\Phi\mathcal{F}_{\mathfrak{c}}(\Phi(X,t)).$
(4.17)

Observe that

$\displaystyle I(t)$
$\displaystyle=a\mathcal{F}_{\mathfrak{c}}(u_{\mathfrak{s}})+X(\mathcal{F}_{\mathfrak{c}}(u_{\mathfrak{s}})-\mathcal{F}_{\mathfrak{c}}(\Phi(X,t)))$
(4.18) $\displaystyle=a\mathcal{F}_{\mathfrak{c}}(u_{\mathfrak{s}})+X\langle
B(\bm{e}_{1},b)+B(b,\bm{e}_{1}),\bm{e}_{1}\rangle,$

and that

$\displaystyle aI(t)$
$\displaystyle\leq|\mathcal{F}_{\mathfrak{c}}(u_{\mathfrak{s}})|\,|a|^{2}+C|a|\,|X|\|b\|_{V}$
(4.19) $\displaystyle\leq
C\|u_{\mathfrak{s}}\|_{V}|a|^{2}+C\mathfrak{C}|a|\,|X|\|b\|_{\bm{\nu}},$

due to the definition of $\mathcal{F}_{\mathfrak{c}}$, and the norm
equivalence between $\|\cdot\|_{\bm{\nu}}$ and $\|b\|_{V}$ (Condition 1). In
the sequel, $C>0$ will denote a generic constant that is allowed to change in
the course of the estimates.

Now let us take $\Omega^{*}$ be the subset of $\Omega$ over which the a priori
bounds (4.11) and (3.22) about $(x,u_{\mathfrak{s}})$ and $(X,\Phi)$, hold.
Then over such a subset of events and for any $t$ in $[0,T/\epsilon]$, we
deduce from (4.19), that

$a(t)I(t)\leq C\epsilon a(t)^{2}+\frac{(1-\eta)}{2}\|b(t)\|_{\bm{\nu}}^{2},$
(4.20)

by using Young’s inequality, (4.11) to control $\|u_{\mathfrak{s}}\|_{V}$, and
(3.22) to control $X$. The factor $(1-\eta)/2$ in (4.20) is chosen through
application of the Young’s inequality, in order to be appropriately “absorbed”
by the linear damping effects on the high-modes, namely

$\displaystyle\langle L_{\lambda}^{\mathfrak{s}}b,b\rangle$
$\displaystyle=\langle L_{\lambda}b,b\rangle$ (4.21)
$\displaystyle=\langle(-\eta A+P_{\lambda})b,b\rangle-(1-\eta)\langle
Ab,b\rangle$ $\displaystyle\leq-\delta\|b\|^{2}-(1-\eta)\|b\|_{\bm{\nu}}^{2},$

where the latter inequality is a consequence of Conditions 1 and 2.

As a result, we have over $\Omega^{*}$ and for any $t$ in $[0,T/\epsilon]$,

$aI+\langle L_{\lambda}^{\mathfrak{s}}b,b\rangle\leq C\epsilon
a^{2}-\delta\|b\|^{2}-\frac{(1-\eta)}{2}\|b\|_{\bm{\nu}}^{2}\leq C\epsilon
a^{2}-\delta\|b\|^{2}.$ (4.22)

Similarly, we have over $\Omega^{*}$ and for any $t$ in $[0,T/\epsilon]$,

$a\Pi_{1}B(u_{s},u_{\mathfrak{s}})\leq C\epsilon^{2}|a|\leq
C\epsilon^{3}+\epsilon a^{2}.$ (4.23)

The estimates (4.22) and (4.23) are used to control $a(t)$. To control $b(t)$
one needs to estimate $\langle J(t),b\rangle$ and
$\langle(x^{2}-X^{2})B_{11}^{\mathfrak{s}},b\rangle$.

Here again by exploiting the a priori bounds (3.22), (4.11) and the Young’s
inequality, we infer that

$\langle(x^{2}-X^{2})B_{11}^{\mathfrak{s}},b\rangle\leq
C\sqrt{\epsilon}|a|\|b\|\leq C\epsilon a^{2}+\frac{\delta}{2}\|b\|^{2},$
(4.24)

and

$\langle J(t),b\rangle\leq C\epsilon^{3/2}\|b\|\leq
C\epsilon^{3}+\frac{\delta}{2}\|b\|^{2},$ (4.25)

still over over $\Omega^{*}$ and for any $t$ in $[0,T/\epsilon]$.

From these estimates, and denoting by
$M=aI+a\Pi_{1}B(u_{s},u_{\mathfrak{s}})$, we observe, using (4.22), that

$M+\langle(x^{2}-X^{2})B_{11}^{\mathfrak{s}},b\rangle+\langle
J,b\rangle+\langle L_{\lambda}^{\mathfrak{s}}b,b\rangle+\epsilon a^{2}\leq
C_{1}\epsilon^{3}+C_{2}\epsilon a^{2},$ (4.26)

where $C_{1}$ and $C_{2}$ are two positive constants.

Following similar steps as in Appendix C.1 (in particular after transformation
as in (C.6)), we arrive then at

$\displaystyle|a(t)|^{2}+\|b(t)\|^{2}$ $\displaystyle\leq e^{C_{2}\epsilon
t+2\sqrt{\epsilon}W_{t}}(|a(0)|^{2}+\|b(0)\|^{2})+C_{1}\epsilon^{3}\int_{0}^{t}e^{C_{2}\epsilon(t-s)+2\sqrt{\epsilon}(W_{t}-W_{s})}\,\mathrm{d}s$
(4.27) $\displaystyle\leq
e^{C_{2}T+2\gamma}(|a(0)|^{2}+\|b(0)\|^{2})+C_{1}\epsilon^{2},$

where $\gamma$ is given by Lemma A.3.

Note that $a(0)=0$. So if we assume $\|b(0)\|\sim\epsilon$, then,

$|a(t)|^{2}+\|b(t)\|^{2}\leq
C\epsilon^{2},\;\;\omega\in\Omega^{*},\;t\in[0,T/\epsilon].$ (4.28)

From this estimate, one can furthermore improve the estimate of $\|b\|$. To do
so, note that by using the Young’s inequality differently than for obtaining
(4.24), leads to

$\langle(x^{2}-X^{2})B_{11}^{\mathfrak{s}},b\rangle\leq
C\sqrt{\epsilon}|a|\|b\|\leq C\epsilon a^{2}+\frac{\delta}{4}\|b\|^{2}.$
(4.29)

Similarly we can arrange the constants through Young’s inequality to get
$\langle J(t),b\rangle\leq C\epsilon^{3/2}\|b\|\leq
C\epsilon^{3}+\delta\|b\|^{2}/4.$

Now since $|a(t)|^{2}\leq C\epsilon^{2}$ due to (4.28), we have

$\langle(x^{2}-X^{2})B_{11}^{\mathfrak{s}},b\rangle\leq
C\epsilon^{3}+\frac{\delta}{4}\|b\|^{2}.$ (4.30)

Then

$\langle(x^{2}-X^{2})B_{11}^{\mathfrak{s}},b\rangle+\langle J,b\rangle+\langle
L_{\lambda}^{\mathfrak{s}}b,b\rangle\leq
C\epsilon^{3}-\frac{\delta}{2}\|b\|^{2},$ (4.31)

which, by application of Gronwall’s inequality (still following similar steps
as in Appendix C.1), gives

$\|b(t)\|^{2}\leq e^{-\delta
t+2\gamma}\|b(0)\|^{2}+C\epsilon^{3},\quad\omega\in\Omega^{*},\;t\in[0,T/\epsilon].$
(4.32)

Thus, assuming $\|b(0)\|\sim\epsilon^{3/2}$, leads finally to

$\|b(t)\|^{2}\leq
C\epsilon^{3},\quad\omega\in\Omega^{*},\;t\in[0,T/\epsilon].$ (4.33)

The desired estimate (4.14) follows from (4.28) and (4.33). The proof is
complete.

∎

######  Remark 4.2.

  * i)

In the condition (LABEL:Eq_cond_init), if we drop the requirement
$\|u_{\mathfrak{s}}(0)-\Phi(X,0)\|\sim\epsilon^{3/2}$, the estimate (4.28)
still holds while (4.33) will be true after skipping a transient time of order
$|\ln(\epsilon)|/\delta$, with $\delta$ given by (4.2a).

Indeed, since it is assumed that $X(0)\sim\sqrt{\epsilon}$, then
$\Phi(X,0)\sim\epsilon$. This together with the condition
$\|u_{\mathfrak{s}}(0)\|\sim\epsilon$ leads to
$\|b(0)\|=\|u_{\mathfrak{s}}(0)-\Phi(X,0)\|\sim\epsilon$. As a result, (4.28)
still holds.

Note also that since $\|b(0)\|\sim\epsilon$, then $e^{-\delta
t+2\gamma}\|b(0)\|^{2}\sim\epsilon^{3}$ for all $t\geq|\ln(\epsilon)|/\delta$.
This together with (4.32) leads to

$\|b(t,\omega)\|^{2}\leq
C\epsilon^{3},\quad\forall\omega\in\Omega^{*},t\in[|\ln(\epsilon)|/\delta,T/\epsilon].$
(4.34)

  * ii)

Note also that if we drop the requirement
$\|u_{\mathfrak{s}}(0)-\Phi(X,0)\|\sim\epsilon^{3/2}$ in the condition
(LABEL:Eq_cond_init), then as a direct reinterpretation of (4.28), the
following relaxed version of the error estimate (4.14) holds without skipping
any transient dynamics:

$\mathbb{P}\Big{(}\sup_{[0,T/\epsilon]}|a(t)|\leq
C\epsilon\;\text{and}\;\sup_{[0,T/\epsilon]}\|b(t)\|\leq C\epsilon\Big{)}\geq
1-\chi.$ (4.35)

  * iii)

A close inspection of the proofs about the error estimates shows that the
requirement $\sigma=\sqrt{\epsilon}$ could be relaxed to the cases
$0<\sigma\leq\sqrt{\epsilon}$.

  * iv)

Note that $\delta<|\beta_{2}(\lambda)|$, as can be seen by replacing $b$ in
(4.21) with $\bm{e}_{2}$ and using $\langle
L_{\lambda}\bm{e}_{2},\bm{e}_{2}\rangle=\beta_{2}(\lambda)\|\bm{e}_{2}\|^{2}$.
We can actually reduce the transient time $|\ln(\epsilon)|/\delta$ in (4.34)
to $|\ln(\epsilon)|/|\beta_{2}(\lambda)|$. This is because the estimate (4.32)
still holds by replacing the exponent $-\delta t$ therein by
$-|\beta_{2}(\lambda)|t$. Indeed, we just need to use $|\beta_{2}(\lambda)|$
in place of $\delta$ when applying the Young’s inequality that leads to the
two estimates (4.29) and (4.30). These together with $\langle
L_{\lambda}^{\mathfrak{s}}b,b\rangle\leq\beta_{2}(\lambda)\|b\|^{2}$ lead to
(4.31) with $\delta$ therein replaced by $|\beta_{2}(\lambda)|$.

##  5\. Applications to Rayleigh-Bénard convection

###  5.1. The stochastic Rayleigh-Bénard model and its mathematical
formulation

We consider the following non-dimensionalized Boussinesq equations driven by a
linear multiplicative noise:

$\displaystyle\mathrm{d}\bm{u}=\big{(}\Delta\bm{u}-\nabla
p+\sqrt{R}\theta\vec{\bm{k}}-(\bm{u}\cdot\nabla)\bm{u}\big{)}\mathrm{d}t+\sigma\bm{u}\circ\mathrm{d}W_{t},$
(5.1)
$\displaystyle\mathrm{d}\theta=\big{(}\Delta\theta+\sqrt{R}w-(\bm{u}\cdot\nabla)\theta\big{)}\mathrm{d}t+\sigma\theta\circ\mathrm{d}W_{t},$
$\displaystyle\text{div\,}\bm{u}=0.$

Here, the unknown functions are the velocity field $\bm{u}=(u,w)$, the
pressure $p$, and the temperature fluctuation $\theta$. $\vec{\bm{k}}=(0,1)$
is the vertical unit vector, $\sigma>0$ is the noise amplitude. The
nondimensional parameter $R$ is the Rayleigh number, which serves as the
control parameter for the bifurcation. We consider the case that the fluid is
confined in a $2D$ nondimensional rectangular domain
$\mathcal{D}=(0,L)\times(0,1)$, where $L>0$ represents the aspect ratio
between the width and the height of the domain. To simplify the presentation,
the above equations (5.1) are presented for the case that the Prandtl number,
Pr, is taken to be one. In the following, we use $(x,z)$ to denote the
coordinates, $x$ for the horizontal direction and $z$ for the vertical
direction.

Various physically sound boundary conditions can be handled in the theoretic
setting; see [MW14, Sec. 4.1.3]. To fix ideas, we consider free-slip boundary
condition for $\bm{u}$, and Dirichlet boundary condition for $T$ on the top
and the bottom boundaries and Neumann boundary condition on the lateral
boundaries:

$\displaystyle u=0,\ \frac{\partial w}{\partial
x}=\frac{\partial\theta}{\partial x}=0$ $\displaystyle\text{at}\ x=0,L,$ (5.2)
$\displaystyle w=\theta=0,\ \frac{\partial u}{\partial z}=0$
$\displaystyle\text{at}\ z=0,1.$

For the problem (5.1)-(5.2), we set the spaces

$\displaystyle H=\\{(\bm{u},\theta)\in(L^{2}(\mathcal{D}))^{3}\ |\
\text{div\,}\bm{u}=0,\bm{u}\cdot\bm{n}|_{\partial\mathcal{D}}=0\\},$ (5.3)
$\displaystyle V=\big{\\{}(\bm{u},\theta)\in(H^{1}(\mathcal{D}))^{3}\cap H\
\big{|}\ \theta=0\text{ at }z=0,1\big{\\}},$ $\displaystyle
H_{1}=\big{\\{}(\bm{u},\theta)\in(H^{2}(\mathcal{D}))^{3}\ \big{|}\
\text{div\,}\bm{u}=0\text{ and \eqref{RBC_bdry condition} holds}\big{\\}},$

where $\bm{n}$ is the unit outward normal vector to $\partial\mathcal{D}$.

Let $L_{R}:H_{1}\rightarrow H$ be defined for any $\bm{\psi}=(\bm{u},\theta)$
in $H_{1}$ by

$\displaystyle L_{R}=-A+P_{R},\text{ with }$ (5.4) $\displaystyle
A\bm{\psi}=(-\mathbb{L}\Delta\bm{u},-\Delta\theta),$ $\displaystyle
P_{R}\bm{\psi}=(\sqrt{R}\,\mathbb{L}(\theta\vec{\bm{k}}),\sqrt{R}\,w),$

where $\mathbb{L}$ denotes the Leray projection defined on
$(L^{2}(\mathcal{D}))^{2}$, which projects each element into the divergence-
free subspace of $(L^{2}(\mathcal{D}))^{2}$. We define also $B:V\times
V\rightarrow H$ as

$B(\bm{\psi}_{1},\bm{\psi}_{2})=(-\mathbb{L}((\bm{u}_{1}\cdot\nabla)\bm{u}_{2}),-(\bm{u}_{1}\cdot\nabla)\theta_{2}),$
(5.5)

for all $\bm{\psi}_{1}=(\bm{u}_{1},\theta_{1})$ and
$\bm{\psi}_{2}=(\bm{u}_{2},\theta_{2})$ in $V$.

We denote by $\|\cdot\|$ the norm on $H$ defined for any
$\bm{\psi}=(u,w,\theta)$ by

$\|\bm{\psi}\|^{2}=|u|_{L^{2}(\mathcal{D})}^{2}+|w|_{L^{2}(\mathcal{D})}^{2}+|\theta|_{L^{2}(\mathcal{D})}^{2},$

and the norm on $V$, the norm $\|\cdot\|_{V}$ defined by

$\|\bm{\psi}\|_{V}^{2}=\|\nabla\bm{\psi}\|+\|\bm{\psi}\|.$

Then the problem (5.1)-(5.2) can be written as

$\displaystyle\mathrm{d}\bm{\psi}=\big{(}L_{R}\bm{\psi}+B(\bm{\psi},\bm{\psi})\big{)}\mathrm{d}t+\sigma\bm{\psi}\circ\mathrm{d}W_{t}.$
(5.6)

In the following, we assume for the non-dimensionalized spatial domain
$\mathcal{D}=[0,L]\times[0,1]$ that $L$ is chosen such that the deterministic
analogue of (5.1) admits a pitchfork bifurcation as $R$ crosses its critical
value $R_{c}$ from below. As pointed out in Sec. 5.2 when verifying Condition
2, there is exactly one critical wave vector for all values of $L$ except
those from a subset of $\mathbb{R}^{+}$ with measure zero; cf. (5.17) for the
latter exceptional set. So the pitchfork bifurcation scenario is a generic
situation for the deterministic case.

###  5.2. Verification of Conditions 1-3 and 1

In this section, we provide details on verifying the Conditions 1-3 and 1.

Verification of Condition 1. Since the linear operator $L_{R}$ defined by
(LABEL:RBC_lin_operators) involves the Leray projection $\mathbb{L}$, it is
not immediately clear that $L_{R}$ is self-adjoint. However, this property
will be verified once the eigenfunctions of $L_{R}$ are computed.

We consider thus the eigenvalue problem $L_{R}\bm{\psi}=\beta\bm{\psi}$ in the
space $H_{1}$. That is

$\displaystyle\mathbb{L}\Delta\bm{u}+\sqrt{R}\,\mathbb{L}(\theta\vec{\bm{k}})=\beta\bm{u},$
(5.7) $\displaystyle\Delta\theta+\sqrt{R}\,w=\beta\theta,$
$\displaystyle\text{div\,}\bm{u}=0.$

Recall that the eigenvalues and eigenfunctions for (5.7) are given as follows
(see e.g. [MW14, Section 4.1.6]):

  * Group one:

$\displaystyle\beta_{0k}=-k^{2}\pi^{2},\quad\bm{e}_{0k}=\sqrt{2/L}(0,0,\sin(k\pi
z)),\;\;k\in\mathbb{N}.$ (5.8)

  * Group two: The eigenvalues are given, for all $j,k$ in $\mathbb{N}$, by

$\displaystyle\beta_{jk}^{\pm}(R)=-\gamma_{jk}^{2}\pm\sqrt{\frac{R\alpha_{j}^{2}}{\gamma_{jk}^{2}}},$
(5.9)

with

$\alpha_{j}=j\pi/L,\;\gamma_{jk}=\sqrt{\alpha_{j}^{2}+k^{2}\pi^{2}};$ (5.10)

and the corresponding eigenfunctions are given by

$\displaystyle\bm{e}^{\pm}_{jk}=\Bigl{(}$ $\displaystyle
a_{jk}^{\pm}\sin(\alpha_{j}x)\cos(k\pi
z),b_{jk}^{\pm}\cos(\alpha_{j}x)\sin(k\pi
z),c_{jk}^{\pm}\cos(\alpha_{j}x)\sin(k\pi z)\Bigr{)},$ (5.11)

where

$a_{jk}^{\pm}=-k\pi N^{\pm}_{jk},\quad
b_{jk}^{\pm}=\alpha_{j}N^{\pm}_{jk},\quad
c_{jk}^{\pm}=\frac{\sqrt{R}\alpha_{j}}{\gamma_{jk}^{2}+\beta_{jk}^{\pm}(R)}N^{\pm}_{jk},$
(5.12)

with

$N^{\pm}_{jk}=\left(\frac{4(\gamma_{jk}^{2}+\beta_{jk}^{\pm}(R))^{2}}{L(\gamma_{jk}^{2}(\gamma_{jk}^{2}+\beta_{jk}^{\pm}(R))^{2}+R\alpha_{j}^{2})}\right)^{1/2}.$
(5.13)

Note that in Group two above, the coefficient
$\frac{1}{\gamma^{2}_{jk}+\beta^{\pm}_{jk}(R)}$ in the third component of
$\bm{e}^{\pm}_{jk}$ is well defined for $j,k$ in $\mathbb{N}$ and for all
$R>0$ because from (5.9), we know that
$\beta_{jk}^{\pm}(R)\neq-\gamma^{2}_{jk}$. The constant $N^{\pm}_{jk}$,
included in each of the three components of $\bm{e}^{\pm}_{jk}$, is a
normalization constant ensuring $\|\bm{e}^{\pm}_{jk}\|=1$.

Now, for any eigenfunctions $\bm{\psi}_{1}$ and $\bm{\psi}_{2}$ listed above,
it follows from a direct calculation that the following identity holds

$\langle
L_{R}\bm{\psi}_{1},\bm{\psi}_{2}\rangle=\langle\bm{\psi}_{1},L_{R}\bm{\psi}_{2}\rangle.$
(5.14)

Thus, $L_{R}$ is self-adjoint.

Verification of Condition 2. Note that only eigenvalues from the second group,
i.e., $\beta_{jk}^{\pm}(R)$ given by (5.9) can change signs when $R$ varies.
To identify the critical value of $R$, we equate $\beta_{jk}^{\pm}(R)$ to zero
to obtain

$R=\gamma_{jk}^{6}/\alpha_{j}^{2}.$ (5.15)

The critical Rayleigh number $R_{c}$ is thus given by:

$\displaystyle R_{c}$
$\displaystyle=\min_{j,k\in\mathbb{N}}\frac{\gamma_{jk}^{6}}{\alpha_{j}^{2}}=\min_{j,k\in\mathbb{N}}\frac{(k^{2}\pi^{2}+j^{2}\pi^{2}/L^{2})^{3}}{j^{2}\pi^{2}/L^{2}}$
(5.16)
$\displaystyle=\min_{j\in\mathbb{N}}\frac{\pi^{4}(1+j^{2}/L^{2})^{3}}{j^{2}/L^{2}}.$

One can readily check that the minimum is achieved at either $(j,k)=(\lfloor
L/\sqrt{2}\rfloor,1)$ or $(\lceil L/\sqrt{2}\rceil,1)$, or both, depending on
the value of $L$, where $\lfloor L/\sqrt{2}\rfloor$ denotes the largest
integer below $L/\sqrt{2}$, and $\lceil L/\sqrt{2}\,\rceil$ denotes the
smallest integer above $L/\sqrt{2}$. The cases with the minimum achieved at
both of the two indices $(j_{1},1)=(\lfloor L/\sqrt{2}\rfloor,1)$ and
$(j_{2},1)=(\lceil L/\sqrt{2}\,\rceil,1)$ occur when
$\gamma_{j_{1},1}^{6}/\alpha_{j_{1}}^{2}=\gamma_{j_{2},1}^{6}/\alpha_{j_{2}}^{2}$.
That is when $L$ is chosen from the following subset of $\mathbb{R}^{+}$:

$\mathcal{S}=\left\\{L>0\;\Big{|}\;\frac{(1+(j_{1})^{2}/L^{2})^{3}}{(j_{1})^{2}}=\frac{(1+(j_{2})^{2}/L^{2})^{3}}{(j_{2})^{2}}\text{
with }j_{1}=\lfloor L/\sqrt{2}\rfloor,j_{2}=\lceil
L/\sqrt{2}\,\rceil,j_{1}\neq j_{2}\right\\}.$ (5.17)

Note that $\mathcal{S}$ has measure zero. In the following, we focus on the
case of pitchfork bifurcation, and assume thus that $L$ takes values outside
of $\mathcal{S}$. We denote the unique index that achieves the minimum value
$R_{c}$ by $(j_{c},1)$.

Thus, Condition 2 is verified with $R$ playing the role of $\lambda$ and
$\beta_{1}(\lambda)$ in (2.9) taken to be $\beta_{j_{c}1}^{+}(R)$ here.
Although the eigenelements here are labeled using a double index associated
with the corresponding wave vectors, they can apparently also be labeled using
a single index to fit into the setting of (2.9).

Verification of Condition 3. Since the spatial domain is taken to be a (2D)
rectangle, then $B$ defined by (5.5) is a continuous bilinear map from
$V\times V$ to $H$ with $V$ being the the subspace of
$(H^{1}(\mathcal{D}))^{3}$ given in (LABEL:RBC_spaces).

Note also that $B$ satisfies that

$\langle B(\bm{\psi}_{1},\bm{\psi}_{2}),\bm{\psi}_{3}\rangle=-\langle
B(\bm{\psi}_{1},\bm{\psi}_{3}),\bm{\psi}_{2}\rangle,$ (5.18)

for all $\bm{\psi}_{1}$, $\bm{\psi}_{2}$, and $\bm{\psi}_{3}$ in $V$. The
condition (3.3) is a direct consequence of (5.18).

It follows also from a direct calculation that

$\displaystyle\langle
B(\bm{e}^{+}_{j_{c}1},\bm{e}^{+}_{j_{c}1}),\bm{e}^{\pm}_{jk})\rangle\\!=\\!0,\;\forall
j,k\in\mathbb{N},$ (5.19) $\displaystyle\langle
B(\bm{e}^{+}_{j_{c}1},\bm{e}^{+}_{j_{c}1}),\bm{e}_{0k})\rangle\\!=\\!\begin{cases}{\displaystyle-\frac{\pi\sqrt{2R}(\gamma_{j_{c}1}^{2}+\beta_{j_{c}1}^{+}(R))\alpha_{j_{c}}^{2}}{\sqrt{L}(\gamma^{2}_{j_{c}1}(\gamma_{j_{c}1}^{2}+\beta_{j_{c}1}^{+}(R))^{2}+R\alpha_{j_{c}}^{2})}}&\text{
if }k=2,\\\ 0&\text{ otherwise.}\end{cases}$

The condition (3.4) is thus verified with $N$ therein equals 2 here, after re-
arranging the double index into a single one.

Verification of Condition 1. Recall that $N$ in (3.11) is $2$ here (see again
(LABEL:Eq_interaction_table_RBC)). From (LABEL:Eq_interaction_table_RBC), we
also know

$B_{11}^{2}=\langle
B(\bm{e}^{+}_{j_{c}1},\bm{e}^{+}_{j_{c}1}),\bm{e}_{02})\rangle=-\frac{\pi\sqrt{2R}(\gamma_{j_{c}1}^{2}+\beta_{j_{c}1}^{+}(R))\alpha_{j_{c}}^{2}}{\sqrt{L}(\gamma^{2}_{j_{c}1}(\gamma_{j_{c}1}^{2}+\beta_{j_{c}1}^{+}(R))^{2}+R\alpha_{j_{c}}^{2})}<0.$
(5.20)

Note also that

$\displaystyle\big{(}B_{12}^{1}+B_{21}^{1}\big{)}$ $\displaystyle=\langle
B(\bm{e}^{+}_{j_{c}1},\bm{e}_{02}),\bm{e}^{+}_{j_{c}1}\rangle+\langle
B(\bm{e}_{02},\bm{e}^{+}_{j_{c}1}),\bm{e}^{+}_{j_{c}1}\rangle$ (5.21)
$\displaystyle=\frac{\pi\sqrt{2R}(\gamma_{j_{c}1}^{2}+\beta_{j_{c}1}^{+}(R))\alpha_{j_{c}}^{2}}{\sqrt{L}(\gamma^{2}_{j_{c}1}(\gamma_{j_{c}1}^{2}+\beta_{j_{c}1}^{+}(R))^{2}+R\alpha_{j_{c}}^{2})}.$

Since $M_{11}^{2,R}$ is always positive, we have then

$B_{11}^{2}M_{11}^{2,R}\big{(}B_{12}^{1}+B_{21}^{1}\big{)}<0.$ (5.22)

Condition 1 is thus verified.

###  5.3. Verification of Conditions 1 and 2

Since it has already been checked above that the linear operator $L_{R}$
defined by (LABEL:RBC_lin_operators) satisfies Conditions 1 and 2, to verify
Condition 1, it remains to show that $\sqrt{\langle
A\bm{\psi},\bm{\psi}\rangle}$ induces a norm on $V$ equivalent to
$\|\bm{\psi}\|_{V}$. For this purpose, we note that for $A$ defined in
(LABEL:RBC_lin_operators), one has, due to the boundary conditions, that for
any $\bm{\psi}=(u,w,\theta)$ in $D(A)=H_{1}$:

$\langle A\bm{\psi},\bm{\psi}\rangle=|\nabla
u|_{L^{2}(\mathcal{D})}^{2}+|\nabla
w|_{L^{2}(\mathcal{D})}^{2}+|\nabla\theta|_{L^{2}(\mathcal{D})}^{2}=\|\nabla\bm{\psi}\|^{2},$

Note that

$\|\nabla\bm{\psi}\|\leq\|\bm{\psi}\|_{V}\leq(1+C)\|\bm{\psi}\|,$

where $C>0$ is a generic constant in the Poincaré inequality. Thus,
$\|\nabla\bm{\psi}\|$ defines a norm equivalent to $\|\bm{\psi}\|_{V}$.

Condition 2 is a consequence of the following lemma.

######  Lemma 5.1.

There exist an $\eta$ in $(0,1)$, $\delta>0$ and $R^{*}>R_{c}$, for which the
condition (4.2) holds for the linear operator $L_{R}$ defined in
(LABEL:RBC_lin_operators) for all $R$ in $[R_{c},R^{*}]$.

See Appendix E for a proof.

###  5.4. Stochastic pitchfork bifurcation in large probability

Due to (LABEL:Eq_interaction_table_RBC), the parameterization $\Phi$ defined
in (3.5) reads here as

$\Phi(X,t,\omega)=\Phi_{2}(X,t,\omega)\bm{e}_{02},\quad\forall
X\in\mathbb{R},$ (5.23)

with

$\Phi_{2}(X,t,\omega)=B_{11}^{2}M_{11}^{2,R}(t,\omega)X^{2},$ (5.24)

and $B_{11}^{2}$ given by (5.20).

The abstract reduced equation (3.8) for the stochastic RBC problem takes then
the following explicit form

$\,\mathrm{d}X=(\beta_{j_{c}1}^{+}(R)X-\alpha(R)M_{11}^{2,R}(t,\omega)X^{3})\,\mathrm{d}t+\sigma
X\circ\,\mathrm{d}W_{t},$ (5.25)

with

$\alpha(R)=-B_{11}^{2}\big{(}B_{12}^{1}+B_{21}^{1}\big{)}=\frac{2\pi^{2}R(\gamma_{j_{c}1}^{2}+\beta_{j_{c}1}^{+}(R))^{2}\alpha_{j_{c}}^{4}}{L(\gamma^{2}_{j_{c}1}(\gamma_{j_{c}1}^{2}+\beta_{j_{c}1}^{+}(R))^{2}+R\alpha_{j_{c}}^{2})^{2}}.$
(5.26)

We call (5.25), the non-Markovian RBC reduced equation.

Note that the RBC reduced equation (5.25) fits into the non-Markovian normal
form of a supercritical pitchfork bifurcation as given by Appendix D. In
particular, the trivial steady state $X=0$ is globally stable when the
Rayleigh number $R$ is below the critical value $R_{c}$ given by (5.16). It
becomes unstable when $R>R_{c}$, and two locally stable random steady states
emerge. These bifurcated random equilibria are given by $\pm X_{R}$, where

$X_{R}(t,\omega)=\frac{1}{\sqrt{2\alpha(R)\int_{-\infty}^{t}M_{11}^{2,R}(s,\omega)\exp(f_{R}(t,s,\omega))\,\mathrm{d}s}},\quad
R\geq R_{c},$ (5.27)

with

$f_{R}(t,s,\omega)=-2\beta_{j_{c}1}^{+}(R)(t-s)-2\sigma(W_{t}(\omega)-W_{s}(\omega)).$
(5.28)

This stochastic pitchfork bifurcation is shown in Fig. 3 at the level of the
probability density function (PDF) of the (random) steady states for a
particular choice of the domain size $L=3$, with $\sigma=0.01$.

Figure 3. The stochastic supercritical pitchfork bifurcation for the RBC
reduced equation (5.25). When $\delta\\!R=R-R_{c}\leq 0$, the probability
density function (PDF) is a Dirac delta function since $X=0$ is globally
stable. When $\delta\\!R>0$, two locally stable random steady states $\pm
X_{R}$ emerge and $X=0$ loses its stability. The most probable values of $\pm
X_{R}$ as $\delta\\!R$ varies form the solid red curves in the
$(\delta\\!R,X)$-plane, and the unstable trivial steady state for
$\delta\\!R>0$ is marked by the dashed red line segment. The most probable
states scale like $\sqrt{\delta\\!R}$, which is consistent with the theoretic
results presented in Lemma D.1, by noting that
$\epsilon=\beta_{j_{c}1}^{+}(R)$ in Lemma D.1 scales linearly with
$\delta\\!R$; cf. (5.32) below. For the parameters, we have set $L=3$, leading
to $R_{c}\approx 660.52$ and the critical wave number $j_{c}=2$; and the noise
amplitude is set to $\sigma=0.01$.

In the following, we provide an interpretation of the general error estimate
results given by Theorem 4.1 for the stochastic RBC problem considered here.
By doing so, we establish a rigorous link between the pitchfork bifurcation
for the non-Markovian RBC reduced equation (5.25) and the dynamics of the full
stochastic RBC problem (5.6) when the noise amplitude parameter $\sigma$
scales like $\sqrt{\delta\\!R}$ with $\delta\\!R=R-R_{c}$. For this purpose,
we introduce

$\displaystyle A_{R}^{+}(t)=X_{R}(t)\bm{e}_{j_{c}1}^{+}+\Phi(X_{R}(t),t),$
(5.29a) $\displaystyle
A_{R}^{-}(t)=-X_{R}(t)\bm{e}^{+}_{j_{c}1}+\Phi(-X_{R}(t),t),$ (5.29b)

where $\bm{e}^{+}_{j_{c}1}$ is the eigenmode that becomes linearly unstable,

$\bm{e}^{+}_{j_{c}1}=\Bigl{(}-\pi\sin(\alpha_{j_{c}}x)\cos(\pi
z),\alpha_{j_{c}}\cos(\alpha_{j_{c}}x)\sin(\pi
z),\frac{\sqrt{R}\alpha_{j_{c}}}{\gamma_{j_{c}1}^{2}+\beta_{j_{c}1}^{+}(R)}\cos(\alpha_{j_{c}}x)\sin(\pi
z)\Bigr{)},$

and the manifold function $\Phi$ is given by (5.23). Then, the following
theorem holds.

######  Theorem 5.1.

Consider the stochastic RBC problem (5.6). Assume that the domain aspect ratio
parameter $L$ is chosen outside of the measure zero set $\mathcal{S}$ given by
(5.17), so that there is exactly one mode that becomes unstable as $R$ crosses
the first critical value $R_{c}$ defined by (5.16). Denote
$\delta\\!R=R-R_{c}$. Assume that $\sigma=(\beta_{j_{c}1}^{+}(R))^{1/2}$ and
that

$X(0)=\langle\bm{\psi}(0),\bm{e}^{+}_{j_{c}1}\rangle\sim\sqrt{\delta\\!R},\;\|\bm{\psi}_{\mathfrak{s}}(0)\|\sim\delta\\!R,\;\textrm{and}\;\|\bm{\psi}_{\mathfrak{s}}(0)\|_{V}\sim\delta\\!R.$
(5.30)

Then, if $X(0)>0$, the solution $\bm{\psi}(t)$ of (5.6) emanating from
$\bm{\psi}(0)$ satisfies that for any $\chi$ in $(0,1)$ and $T>0$, there exist
$R^{*}>R_{c}$ and $C>0$ such that the following estimate holds for all $R$ in
$[R_{c},R^{*}]$:

$\mathbb{P}\Big{(}\sup_{[0,\eta
T/\delta\\!R]}\|\bm{\psi}(t)-A_{R}^{+}(t)\|\leq C\delta\\!R\Big{)}\geq
1-\chi,$ (5.31)

where $A_{R}^{+}$ is defined in (5.29) and $\eta>0$ is a constant depending
only on the aspect ratio $L$. The same estimate holds with $A_{R}^{-}(t)$
replacing $A_{R}^{+}(t)$, if $X(0)<0$.

###### Proof.

The desired result is a direct consequence of Theorem 4.1 and Remark 4.2-(ii).
Indeed, the Conditions 1, 2, 3, and 1 required in Theorem 4.1 are verified for
the RBC problem considered here in Sections 5.2 and 5.3. To derive (5.31), we
are only left to show that the assumption (5.30) on the initial condition is
equivalent to the analogues given in (LABEL:Eq_cond_init) stated in terms of
the parameter $\epsilon=\beta_{j_{c}1}^{+}(R)$ instead.666Note that the
condition $\|u_{\mathfrak{s}}(0)-\Phi(X,0)\|\sim\epsilon^{3/2}$ required in
(LABEL:Eq_cond_init) is not needed here as pointed out in Remark 4.2-(ii).
Namely, we just need to show that $\delta\\!R$ is proportional to $\epsilon$.
The parameter $\eta$ in the estimate (5.31) is related to the associated
proportionality constant when converting from $T/\epsilon$. We are thus only
left with verifying the following scaling relation:

$\beta_{j_{c}1}^{+}(R)=\frac{(j_{c}\pi)^{2}/L^{2}}{(\pi^{2}+(j_{c}\pi)^{2}/L^{2})^{2}}(R-R_{c})+O(|R-R_{c}|^{2}),$
(5.32)

which can be derived through direct algebraic operations. Indeed, by
introducing $\delta_{j_{c}}=\alpha_{j_{c}}^{2}/\gamma_{j_{c}1}^{2}$, we have
for $R$ close to $R_{c}$ that

$\displaystyle\sqrt{\delta_{j_{c}}R}$
$\displaystyle=\sqrt{\delta_{j_{c}}R_{c}+\delta_{j_{c}}(R-R_{c})}=\sqrt{\delta_{j_{c}}R_{c}}\sqrt{1+(R-R_{c})/R_{c}}$
(5.33)
$\displaystyle=\sqrt{\delta_{j_{c}}R_{c}}\left(1+\frac{(R-R_{c})}{2R_{c}}+O(|R-R_{c}|^{2})\right).$

With $\delta_{j_{c}}$ defined above, we can rewrite $\beta_{j_{c}1}^{+}(R)$
given by (5.9) as
$\beta_{j_{c}1}^{+}(R)=-\gamma_{j_{c}1}^{2}+\sqrt{\delta_{j_{c}}R}$, leading
in turn to

$\beta_{j_{c}1}^{+}(R)=-\gamma_{j_{c}1}^{2}+\sqrt{\delta_{j_{c}}R_{c}}+\frac{\sqrt{\delta_{j_{c}}}(R-R_{c})}{2\sqrt{R_{c}}}+O(|R-R_{c}|^{2}).$
(5.34)

Since
$-\gamma_{j_{c}1}^{2}+\sqrt{\delta_{j_{c}}R_{c}}=\beta_{j_{c}1}^{+}(R_{c})=0$,
we get

$\beta_{j_{c}1}^{+}(R)=\frac{\sqrt{\delta_{j_{c}}}(R-R_{c})}{2\sqrt{R_{c}}}+O(|R-R_{c}|^{2}).$
(5.35)

Recall also from (5.16) that $R_{c}=\gamma_{j_{c}1}^{6}/\alpha_{j_{c}}^{2}$.
Using this identity in (5.35), we get
$\beta_{j_{c}1}^{+}(R)=\frac{\alpha_{j_{c}}^{2}}{2\gamma_{j_{c}1}^{4}}(R-R_{c})+O(|R-R_{c}|^{2})$,
and (5.32) follows. ∎

The above theorem shows that the SPDE solution $\bm{\psi}(t)$ of (5.6) is, in
large probability, within an $C\delta\\!R$-cones centered around either
$A_{R}^{+}$ or $A_{R}^{-}$ defined in (5.29) depending on the sign of the
projected initial data $X(0)=\langle\bm{\psi}(0),\bm{e}^{+}_{j_{c}1}\rangle$,
for $t$ in $[0,\eta T/\delta\\!R]$. Recall also from Lemma D.1 (and the
scaling relation (5.32)) that the magnitude of the bifurcated random steady
states $\pm X_{R}$ for the reduced equation (5.25) (and hence the
corresponding lifted states $A_{R}^{\pm}$ in $H$) are of order
$\sqrt{\delta\\!R}$ with large probability for all $t>0$.

As a result, there exists a $\delta\\!R^{\sharp}>0$ such that for any
$\delta\\!R$ in $(0,\delta\\!R^{\sharp})$, the two $C\delta\\!R$-cones, one
centered around $A_{R}^{+}$ and the other around $A_{R}^{-}$, do not intersect
with large probability; moreover, the SPDE solutions emanating from initial
data that satisfy the condition (5.30) stay within one of these two cones with
large probability and over a long time interval of order $1/\delta\\!R$. See
Fig. 4 for a schematic that illustrates this dynamical property in the
projected unstable subspace. This scenario is thus consistent with the idea of
a stochastic pitchfork transition for the SPDE dynamics that takes place when
$R$ crosses the threshold $R_{c}$ from below, albeit subject to large
probability and finite time intervals requirements.

Figure 4. A schematic to illustrate the stochastic pitchfork transition
scenario for the stochastic RBC problem (5.6). The red curves indicate the
most probable values, denoted by $\pm X^{*}_{R}$, of the random steady states
$\pm X_{R}$ for the reduced equation (5.25). Recall that $X_{R}$ scales like
$\sqrt{\delta\\!R}$ with large probability; see again Lemma D.1 and (5.32).
The gray zones indicate the $C\delta\\!R$-cones centered at $\pm X^{*}_{R}$
with the radius $C\delta\\!R$ being the upper bound appearing in the error
estimate (5.31). The vertical dashed line marks out the threshold
$\delta\\!R^{\sharp}$ at which $X^{*}_{R}$ intersects with the line segment
$y=C\delta\\!R$. For $\delta\\!R$ in $(0,\delta\\!R^{\sharp})$, the two
$C\delta\\!R$-cones are disjoint from each other. Theorem 5.1 ensures that for
$\delta\\!R\leq\min\\{\delta\\!R^{\sharp},\delta\\!R^{*}\\}$, the (projected)
SPDE solution from an initial condition satisfying (5.30) is confined in one
of the two cones with large probability and over a large time interval of
order $1/\delta\\!R$. Here, $\delta\\!R^{*}=R^{*}-R_{c}$ with $R^{*}$ being as
specified in Theorem 5.1.

##  6\. Concluding remarks

Based on a dynamical reformulation of the leading-order approximation formulas
for stochastic invariant manifolds, we derived rigorous energy estimates
between the solutions of the corresponding reduced system (3.19) and those of
the original SPDE (3.1). Such estimates as summarised in Theorem 4.1 hold with
large probability and are obtained for the case of a stochastic pitchfork
bifurcation scenario. The key assumption 2 that enables these estimates
exploits dissipation effects brought by the stable modes of the linear
operator in order to suitably counterbalance the loss of regularity due to the
nonlinear terms.

Although seeking a $\mathbb{P}$-almost sure description of the bifurcation for
a given SPDE is still out of reach for most problems arising from
applications777Recall that for stochastic ODEs already the question of how to
describe a stochastic bifurcation, is not completely settled; see e.g. [Arn98,
Chap. 9] and [CF98]., we showed in this article that the dynamical properties
of the SPDE solutions are actually captured by our reduced systems with large
probability and over long time intervals, as long as the noise’s intensity and
the eigenvalue’s magnitude of the mildly unstable mode, scale accordingly.
While this type of characterizations has been investigated in the literature
before, such as via the amplitude equation approach [Blö07], these approaches
do not seem to be directly applicable to problems considered here which are
with quadratic nonlinearities and subject to multiplicative noise. While the
focus here is on the pitchfork bifurcation, our approach can be suitably
adapted to handle stochastic disturbances of fluids experiencing more general
bifurcations such as $S^{1}$-attractor bifurcations [MW05], Hopf, and double-
Hopf bifurcations.

We have also shown that the underlying analytic formulas of our reduced
systems near the onset of instability for SPDE models involve, unlike for
PDEs, non-Markovian coefficients which depend explicitly on the noise path
(and its “past”), the model’s coefficients, the noise intensity, as well as
the model’s eigenelements. In that respect, the memory functions captured by
these coefficients are not set arbitrarily, but derived in a consistent
manner, so that in particular the characteristics of the memory function is
self-consistently determined by the intensity of the random force from the
model’s equation.

As aready pointed out in [CLW15b], higher-order approximations of the genuine
stochastic center manifold, $h^{\lambda}_{\omega}(X,t)$, involve also memory
coefficients albeit of more complicated functional dependence on the noise
path. The underlying high-order non-Markovian parameterizations have also
analytic expressions that can be derived by solving the appropriate backward-
forward (BF) systems [CLW15b, Chaps. 4 and 7.3]. The resulting random
coefficients solve then new auxiliary SDEs that can also be used for deriving
the error estimates as done for the leading-order approximation. Such high-
order non-Markovian parameterizations should lead to improved error estimates,
as numerically shown in [CLW15b, Chap. 7] in a stochastic Burgers example; see
also [CL15].

Finally and more generally, we mention that the usage of such backward-forward
systems allows also for moving away from criticality, by adopting and adapting
the variational approach and theory of optimal parameterizing manifolds (OPMs)
of [CLM17, CLM20, CLM21] to the stochastic context. The new parameterizations
obtained this way can be viewed as optimized homotopic deformations of those
valid near criticality recalled in Theorem 2.1 above. In the OPM approach, the
optimization step is performed using data of the fully resolved problem (in a
parsimonious way) by minimizing a least-square metric measuring typically the
parameterization defect, and where the parameters to optimize are reduced to
the backward integration times $\tau$ of the BF systems involved; see [CLM20,
Sec. 4.3]. In the stochastic context, it means that the $M_{n}$-coefficients
in (2.17) are replaced away from criticality by integrals of the form
$\int_{-\tau}^{0}e^{g_{n}(\lambda)s+\sigma W_{s}(\omega)}\,\mathrm{d}s$, that
are optimized (in $\tau$) per mode $\bm{e}_{n}$. Details about such an OPM
approach for the stochastic context will be communicated elsewhere.

## Acknowledgments

This work has been supported by the Office of Naval Research (ONR)
Multidisciplinary University Research Initiative (MURI) grant
N00014-20-1-2023, and by the National Science Foundation grant DMS-2108856.
This study was also supported by the European Research Council (ERC) under the
European Union’s Horizon 2020 research and innovation program [Grant Agreement
No. 810370], the Ben May Center grant for theoretical and/or computational
research and by the Israeli Council for Higher Education (CHE) via the
Weizmann Data Science Research Center.

##  Appendix A $M$-term: First moments, non-Gaussian statistics, and large
probability estimates

###  A.1. First moments

We consider the random variable of the form

$M(\omega)=\int_{-\infty}^{0}e^{g(\lambda)s+\sigma
W_{s}(\omega)}\,\mathrm{d}s,$ (A.1)

in which $\sigma$ is allowed to take negative values and $g(\lambda)$ is a
function of the control parameter $\lambda$. Such random variables (with
negative $\sigma$) arise for instance when the stable modes are forced
stochastically but not the unstable ones; see Eq. (2.25).

The following lemma provides exact formulas for the expectation, variance and
autocorrelation function, of such a general $M$-term.

######  Lemma A.1.

Let $g>0$ and

$\qquad\sigma_{*}=\sqrt{2g},\qquad\sigma_{\\#}=\sqrt{g}.$ (A.2)

Then $M$-term defined by (A.1) a wide-sense stationary random process provided
that $|\sigma|<\sigma_{\\#}$, and furthermore

  * (i)

The expectation of $M$ exists if and only if $|\sigma|<\sigma_{*}$, and is
given by

$\mathbb{E}(M)=\frac{2}{2g-\sigma^{2}},\quad|\sigma|<\sigma_{*}.$ (A.3)

  * (ii)

The variance of $M$ exists if and only if $|\sigma|<\sigma_{\\#}$, and is
given by

$\mathrm{Var}(M)=\frac{2\sigma^{2}}{(2g-\sigma^{2})^{2}(g-\sigma^{2})}.$ (A.4)

Finally, the autocorrelation $R(t)$ of the stochastic process

$M(t,\omega)=\int_{-\infty}^{0}e^{gs+\sigma(W_{s+t}(\omega)-W_{t}(\omega))}\,\,\mathrm{d}s,$
(A.5)

exists if and only if $|\sigma|<\sigma_{\\#}$, and is given by

$R(t)=\exp\Bigl{(}-\Bigl{(}g-\frac{\sigma^{2}}{2}\Bigr{)}|t|\Bigr{)},\quad
t\in\mathbb{R}.$ (A.6)

This lemma results from direct application of the Fubini Theorem, the
independent increment property of the Wiener process, and the fact that
$\mathbb{E}(e^{\sigma W_{t}(\cdot)})=e^{\sigma^{2}|t|/2}$ for any
$t\in\mathbb{R}$, as expectation of the geometric Brownian motion generated by
$\,\mathrm{d}S_{t}=\frac{\sigma^{2}}{2}S_{t}\,\mathrm{d}t+\sigma
S_{t}\,\mathrm{d}W_{t}$. See the Appendix of [CLW15b] for a proof.

###  A.2. Non-Gaussian statistics and log-normal approximations

We provide here a simple analytic understanding regarding the non-Gaussian
statistics followed by an arbitrary $M$-term, in the case of real eigenvalues.
Recall that an $M$-term is obtained as the stationary solution of the scalar
SDE,

$\displaystyle\mathrm{d}M=\left(1-gM\right)\mathrm{d}t-\sigma
M\circ\mathrm{d}W_{t},$ (A.7)

where $g$ denotes a distance to resonance, namely
$g=(\bm{k},\bm{\beta}_{c}(\lambda))-\beta_{n}(\lambda)$ for some $n$.

Rescale (A.7) by using

$W_{t}(\omega)=\frac{1}{\sqrt{g}}\widetilde{W}_{gt}(\omega),\quad\widetilde{t}=gt,\quad\widetilde{M}=gM.$
(A.8)

Then (A.7) becomes

$\,\mathrm{d}\widetilde{M}=(1-\widetilde{M})\,\mathrm{d}\widetilde{t}-\frac{\sigma}{\sqrt{g}}\widetilde{M}\circ\,\mathrm{d}\widetilde{W}_{\widetilde{t}}.$
(A.9)

To have an explicit formula for the probability distribution function (PDF) of
the random variable $\widetilde{M}$ is non-trivial. We illustrate nevertheless
that this PDF is well approximated by the PDF of an analytic, log-normal
distributed random variable $N$, provided that $\sigma<\sqrt{g}$.

Our analysis revealed that the parameters of this log-normal distribution are

$\displaystyle\mu_{\sigma,g}=\ln(\kappa_{\sigma,g})-\frac{1}{2}(\Sigma_{\sigma,g})^{2},$
(A.10)
$\displaystyle\Sigma_{\sigma,g}=\sqrt{\ln\left(1+\left(\frac{\eta_{\sigma,g}}{\kappa_{\sigma,g}}\right)^{2}\right)},$

where

$\displaystyle\kappa_{\sigma,g}=\frac{2}{2-\sigma^{2}/g},$ (A.11)
$\displaystyle\eta_{\sigma,g}=\sqrt{\frac{2\sigma^{2}/g}{(2-\sigma^{2}/g)^{2}(1-\sigma^{2}/g)}}.$

The PDF of $N$ is then given by ([JKB94, Chap. 14])

$p(N)=\frac{1}{N\Sigma_{\sigma,g}\sqrt{2\pi}}\exp\left(-\frac{(\ln(N)-\mu_{\sigma,g})^{2}}{2\Sigma_{\sigma,g}^{2}}\right),\quad
N>0,$ (A.12)

and turns out to approximate robustly the distribution of simulated $M$-term
for different values of $\sigma/\sqrt{g}$ as shown in Fig. 5.

Figure 5. PDF of $\widetilde{M}$ (black) vs. its lognormal approximation given
by (A.12) (red dash curve), for various values of $\sigma/\sqrt{g}$. A semi-
log scale is used to plot these PDFs. Note that $\sigma/\sqrt{g}<1$ is
required to ensure the mean and variance of $\widetilde{M}$ to exist; see
Lemma A.1.

###  A.3. Probability estimates in the large

######  Lemma A.2.

We assume that Conditions 1 and 2 as well as (3.4) of Sec. 3 hold. Consider
the $M_{11}^{n,\lambda}$-terms appearing in (3.6) for $n=2,\cdots,N$ and
$\lambda\in[\lambda_{c},\lambda^{*}]$, where $N$ is given by (3.4) and
$\lambda^{*}$ is given in Condition 2. Then, for any $\chi$ in $(0,1)$, there
exists $\kappa>0$ such that

$\mathbb{P}\left\\{0\leq
M_{11}^{n,\lambda}(\omega)<\kappa\;\big{|}\;n=2,\cdots,N,\lambda\in[\lambda_{c},\lambda^{*}]\right\\}\geq
1-\chi.$ (A.13)

###### Proof.

The result is a direct consequence of the Chebyshev’s inequality, which states
that if $Y$ is a random variable and $p\geq 1$, then for any $\kappa>0$, it
holds that

$\mathbb{P}\\{\omega:|Y|\geq\kappa\\}\leq\frac{1}{\kappa^{p}}\mathbb{E}(|Y|^{p}).$
(A.14)

Since
$M_{11}^{n,\lambda}(\omega)=\int_{-\infty}^{0}e^{(2\epsilon-\beta_{n}(\lambda))s+\sigma
W_{s}(\omega)}\mathrm{d}s$ and the eigenvalues are all real-valued and are
arranged in descending order, we have for all
$\omega\in\Omega,\lambda\in[\lambda_{c},\lambda^{*}]$,

$0\leq M_{11}^{N,\lambda}(0,\omega)\leq
M_{11}^{N-1,\lambda}(0,\omega)\leq\cdots\leq M_{11}^{3,\lambda}(0,\omega)\leq
M_{11}^{2,\lambda}(0,\omega).$ (A.15)

Denote

$\lambda^{\sharp}=\mathop{\mathrm{argmax}}_{\lambda\in[\lambda_{c},\lambda^{*}]}\beta_{2}(\lambda).$
(A.16)

We have then

$M_{11}^{2,\lambda}(\omega)\leq
M_{11}^{2,\lambda^{\sharp}}(\omega),\quad\forall\omega\in\Omega,\lambda\in[\lambda_{c},\lambda^{*}].$
(A.17)

Recall that due to Lemma A.1

$\mathbb{E}(M_{11}^{2,\lambda^{\sharp}})=\frac{2}{2(2\epsilon-\beta_{2}(\lambda^{\sharp}))-\sigma^{2}},$
(A.18)

which holds for all $\sigma<\sqrt{2(2\epsilon-\beta_{2}(\lambda^{\sharp}))}$.
Since $\sigma=\sqrt{\epsilon}$ and $\beta_{2}(\lambda^{\sharp})<0$ thanks to
Condition 2, this desired condition
$\sigma<\sqrt{2(2\epsilon-\beta_{2}(\lambda^{\sharp}))}$ is satisfied. Thus,
for any $\chi$ in $(0,1)$, thanks to (A.15) and (A.17), we just need to apply
(A.14) with $Y=M_{11}^{2,\lambda^{\sharp}}$, $p=1$, and

$\kappa=\frac{1}{|\beta_{2}(\lambda^{\sharp})|\chi},$ (A.19)

to ensure (A.13). ∎

######  Lemma A.3.

For any $\epsilon>0$, $T>0$ and $\chi$ in $(0,1]$, let
$\gamma=\sqrt{-2T\ln(1-\chi))}$, then

$\mathbb{P}\left\\{\sqrt{\epsilon}\sup_{0\leq t\leq
T/\epsilon}|W_{t}(\omega)|\leq\gamma\right\\}\geq 1-\chi.$ (A.20)

###### Proof.

Recall that for any $\gamma>0$ and $T^{*}>0$, it holds that

$\mathbb{P}\\{\sup_{0\leq t\leq T^{*}}|W_{t}(\omega)|\geq\gamma\\}\leq
2\mathbb{P}\\{|W_{T^{*}}(\omega)|\geq\gamma\\};$ (A.21)

see e.g. Billingsley [Bil86] page 529. So for any $\sigma>0$, we have

$\mathbb{P}\bigg{\\{}\sigma\sup_{0\leq t\leq
T^{*}}|W_{t}(\omega)|\geq\gamma\bigg{\\}}\leq
2\mathbb{P}\left\\{\sigma|W_{T^{*}}(\omega)|\geq\gamma\right\\}\leq\exp\left(-\frac{\gamma^{2}}{2\sigma^{2}T^{*}}\right),$
(A.22)

where in the last inequality above we used the facts that $W_{T^{*}}(\omega)$
follows a normal distribution with mean zero and variance equal to $T^{*}$ and
that

$\int_{a}^{\infty}e^{-x^{2}}\,\mathrm{d}x\leq\frac{\sqrt{\pi}}{2}e^{-a^{2}}.$
(A.23)

Hence, with $T^{*}=T/\epsilon$ and $\sigma=\sqrt{\epsilon}$, we have

$\mathbb{P}\left\\{\sqrt{\epsilon}\sup_{0\leq t\leq
T/\epsilon}|W_{t}(\omega)|\geq\gamma\right\\}\leq\exp\left(-\frac{\gamma^{2}}{2T}\right).$
(A.24)

Then, by choosing $\gamma=\sqrt{-2T\ln(1-\chi))}$, we get the desired estimate
(A.20). ∎

##  Appendix B A priori estimates in large probability: Surrogate system
(3.19)

###### Proof of Lemma 3.1.

The desired estimate (3.22) follows from (i) basic probabilistic estimates
about the $M$-terms and the Brownian motion (see Lemmas A.2 and A.3 above),
and (ii) elementary energy estimates based on the dissipation condition
(3.11).

To perform (and simplify) our estimates we first transform the SDE system
(3.19) into differential equations with random coefficients by using