Selection of quasi-stationary states in the stochastically forced Navier-Stokes equation on the torus

TL;DR

This paper develops a finite-dimensional model to understand how the parameter delta influences the selection of quasi-stationary states in the stochastically forced 2D Navier-Stokes vorticity equation on a torus, validated through simulations and analytical techniques.

Contribution

It introduces a center manifold projection capturing the state selection mechanism and applies averaging and homogenization to analyze long-term dynamics.

Findings

Model agrees with full vorticity equation in key long-time behaviors

Parameter delta critically influences state selection

Analytical techniques elucidate dominant dynamics for delta near 1

Abstract

The stochastically forced vorticity equation associated with the two dimensional incompressible Navier-Stokes equation on $D_\delta:=[0,2\pi\delta]\times [0,2\pi]$ is considered for $\delta\approx 1$, periodic boundary conditions, and viscocity $0<\nu\ll 1$. An explicit family of quasi-stationary states of the deterministic vorticity equation is known to play an important role in the long-time evolution of solutions both in the presence of and without noise. Recent results show the parameter $\delta$ plays a central role in selecting which of the quasi-stationary states is most important. In this paper, we aim to develop a finite dimensional model that captures this selection mechanism for the stochastic vorticity equation. This is done by projecting the vorticity equation in Fourier space onto a center manifold corresponding to the lowest eight Fourier modes. Through Monte Carlo simulation, the vorticity equation and the model are shown to be in agreement regarding key aspects of the long-time dynamics. Following this comparison, perturbation analysis is performed on the model via averaging and homogenization techniques to determine the leading order dynamics for statistics of interest for $\delta\approx1$.

Figures

Equations

$$\frac{\partial \mathbf{u}}{\partial t}=\nu \Delta \mathbf{u}$$

$$\nabla \cdot \mathbf{u}$$

$${\partial }_t\omega =\nu \Delta \omega -\mathbf{u}\cdot \nabla \omega ,\mathbf{u}=\left(\begin{array}{c}{\partial }_y(-{\Delta }^{-1})\\ -{\partial }_x(-{\Delta }^{-1})\end{array}\right)\omega .$$

$${\partial }_t\omega =\nu \Delta \omega -\mathbf{u}\cdot \nabla \omega +\frac{\partial \mathcal{W}}{\partial t},\mathbf{u}=\left(\begin{array}{c}{\partial }_y(-{\Delta }^{-1})\\ -{\partial }_x(-{\Delta }^{-1})\end{array}\right)\omega .$$

$$\mathcal{W}(t,x,y)=\sqrt{2\nu }\sum _{\vec{k}\in \mathcal{K}\subset {\mathbb{Z}}^2\backslash \{(0,0)\}}{\sigma }_{\vec{k}}{e}^{\mathrm{i}(k_{1}x/\delta +k_{2}y)}{\beta }_{\vec{k}}(t),$$

$$\omega (x,y,t)=e^{-\frac{\nu }{{\delta }^2}t}[a_1\cos (x/\delta )+a_2\sin (x/\delta )]+e^{-\nu t}[a_3\cos (y)+a_4\sin (y)].$$

$${\omega }_{xbar}(x,t)=e^{-\frac{\nu }{{\delta }^2}t}\sin (x/\delta ),{\omega }_{ybar}(y,t)=e^{-\nu t}\sin y,$$

$$\mathbf{u}_{xbar}(x,t)=-\delta e^{-\frac{\nu}{\delta^{2}}t}\left(\begin{array}[]{c}0\\ \cos(x/\delta)\\ \end{array}\right),\qquad\mathbf{u}_{ybar}(y,t)=e^{-\nu t}\left(\begin{array}[]{c}\cos y\\ 0\\ \end{array}\right),$$

$${\omega }_{dipole}(x,y,t)=e^{-\frac{\nu }{{\delta }^2}t}\sin (x/\delta )+e^{-\nu t}\sin y,$$

$$\mathbf{u}_{dipole}(x,y,t)=\left(\begin{array}[]{c}e^{-\nu t}\cos y\\ -\delta e^{-\frac{\nu}{\delta^{2}}t}\cos(x/\delta)\\ \end{array}\right).$$

$$\omega (x,y)=\sum _{\vec{k}\ne (0,0)}{\hat{\omega }}_{\vec{k}}{e}^{\mathrm{i}(k_{1}x/\delta +k_{2}y)},{\hat{\omega }}_{\vec{k}}=\frac{1}{4{\pi }^2\delta }{\int }_{D_{\delta }}\omega (x,y)e^{-\mathrm{i}(k_{1}x/\delta +k_{2}y)}\mathrm{d}x\mathrm{d}y,$$

$${\dot{\hat{\omega }}}_{\vec{k}}$$

$$=-\frac{\nu }{{\delta }^2}|\vec{k}{|}_{\delta }^2{\hat{\omega }}_{\vec{k}}-\frac{\delta }{2}\sum _{\vec{j}+\vec{l}=\vec{k}}⟨{\vec{j}}^{\perp },\vec{l}⟩\left(\frac{1}{|\vec{l}{|}_{\delta }^2}-\frac{1}{|\vec{j}{|}_{\delta }^2}\right){\hat{\omega }}_{\vec{j}}{\hat{\omega }}_{\vec{l}}+\sqrt{2\nu }{\sigma }_{\vec{k}}{\dot{\beta }}_{\vec{k}},$$

$$|\vec{k}{|}_{\delta }^2=k_1^2+{\delta }^2{k}_2^2,{\vec{k}}^{\perp }=(k_2,-k_1).$$

$$Z_{vort}(t):=\frac{|{\hat{\omega }}_{(1,0)}(t){|}^2}{|{\hat{\omega }}_{(1,0)}(t){|}^2+|{\hat{\omega }}_{(0,1)}(t){|}^2},$$

$${\omega }_1$$

$${\omega }_5$$

$${\omega }_1={\bar{\omega }}_2,{\omega }_3={\bar{\omega }}_4,{\omega }_5={\bar{\omega }}_8,{\omega }_7={\bar{\omega }}_8.$$

$$\mathcal{M}=\{\hat{\omega }:{\hat{\omega }}_{\vec{k}}=H({\omega }_1,\dots ,{\omega }_8;\vec{k}),\vec{k}\notin {\mathcal{K}}_0\}$$

$${\dot{\omega }}_1$$

$${\dot{\omega }}_3$$

$${\dot{\omega }}_5$$

$${\dot{\omega }}_7$$

$$Z_{red}(t):=\frac{|{\omega }_1(t){|}^2}{|{\omega }_1(t){|}^2+|{\omega }_3(t){|}^2},$$

$$|{\sigma }_{\vec{k}}|\le e^{-{\alpha }_0|\vec{k}{|}^2}.$$

$${\bar{Z}}_{vort}(t)=\frac{1}{N}\sum _{i=1}^N{Z}_{vort}^i(t),{\bar{Z}}_{red}(t)=\frac{1}{N}\sum _{i=1}^N{Z}_{red}^i(t).$$

$$V_{vort}(t)=\frac{1}{N-1}\sum _{i=1}^N(Z_{vort}^i(t)-{\bar{Z}}_{vort}(t){)}^2,V_{red}(t)=\frac{1}{N-1}\sum _{i=1}^N(Z_{red}^i(t)-{\bar{Z}}_{red}(t){)}^2.$$

$$A(f,t_{burn}):=\frac{1}{T-t_{burn}}{\int }_{t_{burn}}^{T}f(t)\text{mbox}\mathrm{d}t.$$

$$C{I}^{\pm }(t)={\bar{Z}}_{vort}(t)\pm 1.96\ast \sqrt{\frac{V_{vort}(t)}{N}}.$$

$${\dot{\omega }}_1$$

$${\dot{\omega }}_3$$

$${\dot{\omega }}_5$$

$$-\frac{{\delta }^6(3+{\delta }^2)}{2\nu (4+{\delta }^2)(1+{\delta }^2)}{\omega }_5{\omega }_1^2-\frac{1+3{\delta }^2}{2\nu {\delta }^2(1+4{\delta }^2)(1+{\delta }^2)}{\omega }_5{\omega }_3^2+\sqrt{2\nu }{\sigma }_5{\dot{W}}_5$$

$${\dot{\omega }}_7$$

$$-\frac{{\delta }^6(3+{\delta }^2)}{2\nu (4+{\delta }^2)(1+{\delta }^2)}{\omega }_7{\omega }_1^2-\frac{1+3{\delta }^2}{2\nu {\delta }^2(1+4{\delta }^2)(1+{\delta }^2)}{\omega }_7{\omega }_3^2+\sqrt{2\nu }{\sigma }_7{\dot{W}}_7.$$

$${\sigma }_{1,3}=e^{-{\alpha }_0}\text{mbox}and{\sigma }_{5,7}=e^{-2{\alpha }_0}.$$

Full text

1 Introduction

Consider the 2D incompressible Navier-Stokes equation,

[TABLE]

on the possibly asymmetric torus $(x,y)\in D_{\delta}:=[0,2\pi\delta]\times[0,2\pi]$ with $\delta\approx 1$, periodic boundary conditions, and viscosity $0<\nu\ll 1$. To obtain the equivalent vorticity formulation of the equation, take the curl of the vector field $\mathbf{u}$ and set $\omega=(0,0,1)\cdot(\nabla\times\mathbf{u})$ to find

[TABLE]

The relation between $\mathbf{u}$ and $\omega$ is known as the Biot-Savart law. The periodic boundary conditions insure $\int_{D_{\delta}}\omega=0$, and therefore $\Delta^{-1}\omega$ is well-defined.

Adding random forcing to the system allows to account for stochasticity/genericity in the system, see for example [17, 5, 16]. In particular, we add a stochastic forcing term to (1.2) to obtain the stochastic 2D vorticity equation,

[TABLE]

The noise is white in time, colored in space, and takes the form, for $\vec{k}=(k_{1},k_{2})\neq(0,0)$,

[TABLE]

with spatial correlation $\sigma_{\vec{k}}$ and $\mathcal{K}$ to be commented on below. Here $\beta(t)=\{\beta_{\vec{k}}(t)\}$ is a collection of i.i.d. Wiener processes.

Notice that with the noise (1.4), the equation (1.3) is now stochastic. To insure the random vorticity remains real valued for all times $t\geq 0$, the following complex conjugacy conditions are imposed, $\bar{\sigma}_{\vec{k}}=\sigma_{-\vec{k}}$ and $\bar{\beta}_{\vec{k}}=\beta_{-\vec{k}}$. Additional assumptions are often placed on the noise coefficients, $\sigma_{\vec{k}}$, to insure certain smoothness properties of solutions. In particular, we assume that there exist fixed positive constants $C_{0}$ and $\alpha_{0}$ such that $|\sigma_{\vec{k}}|\leq C_{0}e^{-\alpha_{0}|\vec{k}|^{2}}$ so that solutions will then be analytic in space [14]. Since the boundary conditions force solutions of the deterministic equation to satisfy $\int_{D_{\delta}}\omega=0$, we choose $\sigma_{(0,0)}=0$ so this property is preserved. Note that if $\sigma_{\vec{k}}=0$ for all $\vec{k}\in\mathbb{Z}^{2}$ then (1.3) reduces to the deterministic vorticity equation.

Although an $L^{2}$ energy estimate shows solutions of (1.2) have a time-asymptotic rest state of zero, certain quasi-stationary states, known as bars and dipoles, rapidly attract nearby solutions and correspond to transient structures that play a key role in the long-time evolution of solutions [1, 2, 6, 8, 9, 12, 19, 20]. These quasi-stationary states are members of an explicit family of functions given by,

[TABLE]

If $\delta=1$, then any member of this family is an exact solution to the deterministic vorticity equation. If $\delta\neq 1$, then (1.5) remains a solution if and only if $a_{1}=a_{2}=0$ or if $a_{3}=a_{4}=0$. These members, which only depend on one spatial variable, are called bar states, and they are also known as unidirectional or Kolmogorov flow. The $x$- and $y$-bar states are members of this family given by

[TABLE]

or similarly with sine replaced by cosine. The associated velocity fields are given by

[TABLE]

respectively. The dipoles are also members of the family (1.5) and are given by

[TABLE]

or similarly with sine replaced by cosine, with velocity field

[TABLE]

For illustration, contour plots for the bar and dipole states for fixed $t=0$ on the symmetric torus ($\delta=1$) are shown in Figure 1.

When $\nu=0$, equation (1.1) becomes the Euler equation. It is reasonable to expect that stationary solutions of the Euler equation could play an critical role in the evolution of the Navier-Stokes equation for $0<\nu\ll 1$. However, there are infinitely many stationary solutions, including the bars and dipoles, and so it is not immediately clear how to determine which would be most important. In [20], entropy arguments and extensive numerical studies were conducted in the case $\delta=1$ and suggested that the bars and dipoles should be the two most important stationary solutions of the Euler equations. Although both states were observed after initial transient periods in the evolution of the Navier-Stokes equation, interestingly the dipole seemed to emerge for a large class of initial data, whereas the bar states only emerged for a special class of initial data. Subsequent work, again for the deterministic system, showed that indeed the bar states attract nearby solutions at a rate much faster than the background global decay rate, confirming their importance as quasi-stationary states. Results in the case $\delta=1$ can be found in [2, 9] and results for more general values of $\delta$ are in [12, 19]. The stochastic system (1.3) was numerically analyzed in [6] where, after an initial transient period, metastable switching between the bars and dipoles was seen, with the dipole being dominant for $\delta=1$ and the bar states being dominant for $\delta\neq 1$.

In this paper we develop a low dimensional model that captures how the dominant quasi-stationary state in the stochastically forced Navier-Stokes equation is selected by the aspect ratio of the spatial domain, $\delta$. Among the existing results, those that most greatly motivate this paper can be found in [1] and [6]. The results of the latter paper [6], briefly described above, to our knowledge were the first to suggest that $\delta$ could provide such a selection mechanism. The former paper [1] was our previous work focusing on the deterministic vorticity equation, (1.2), in which we derived a finite-dimensional model that captured the selection mechanism via the parameter $\delta$. We now seek to use that same finite-dimensional model, but with the addition of noise, to numerically investigate the selection mechanism for the stochastic equation (1.3). Indeed, one can see from the sample paths of Figure 6(a), that individual sample paths exhibit transitions between $x$-bar and $y$-bar states, as it has also been observed in [6].

The rest of the paper is organized as follows. In §2 we review the finite-dimensional model originally derived in [1] and the theoretical results of that work regarding the selection mechanism in the deterministic setting. In addition, we also add noise to that model to obtain the stochastic differential equation (SDE) model that is the focus of this current work. In §3, to determine the validity of the SDE model, we compare statistics related to a direct simulation of the stochastic vorticity equation (1.3) with those of the SDE. We demonstrate numerically that the statistics of the two equations agree in all cases, $\delta>1$, $\delta<1$ and $\delta=1$. In particular, solutions to both systems evolve towards an $x$-bar, $y$-bar, and dipole in the three respective cases. In §4, we further examine the SDE model by viewing it as a perturbation in the limit as $|\delta^{2}-1|$ and $\nu$ converge to zero. We show that, after appropriate time-space rescalings, the system can be viewed as a slow-fast system and classical averaging and homogenization techniques apply. Via the backward Kolmogorov equation, a system of PDEs that governs the leading order dynamics of a key order parameter, $\mathbb{E}[Z_{red}(t)]$, defined in (2.7), is derived. This gives us an additional formal approximation to the expected value of the order parameter, which we can use to show the selection of the quasi-stationary state. Numerically solving the PDEs allows us to approximate the evolution of $\mathbb{E}[Z_{red}(t)]$ for values of $\delta$ close to 1, at least on some initial finite interval of time. Conclusions and future directions then follow in §5.

2 Fourier space representation and model reduction

Due to the form of the family of solutions (1.5), it is most convenient to express the stochastic vorticity equation in Fourier space. Hence, letting

[TABLE]

we obtain, for $\vec{j}$, $\vec{k}$ and $\vec{l}\neq(0,0)$, the following system of infinitely many coupled SDEs,

[TABLE]

where

[TABLE]

Viewing the system in Fourier space allows us to use the relative energy in certain modes to measure the proximity of solutions to an $x$-bar, $y$-bar, or dipole state. The $x$-bar states, $e^{-\frac{\nu}{\delta^{2}}t}\cos(x/\delta)$ and $e^{-\frac{\nu}{\delta^{2}}t}\sin(x/\delta)$, correspond to solutions with energy only in the $\vec{k}=(\pm 1,0)$ modes and the $y$-bar states, $e^{-\nu t}\cos(y)$ and $e^{-\nu t}\sin(y)$, correspond to solutions with energy only in the $\vec{k}=(0,\pm 1)$ modes. Solutions with energy in both the $\vec{k}=(\pm 1,0)$ and $\vec{k}=(0,\pm 1)$ modes correspond to the dipole state. These four modes are the lowest in the system and will be referred to as the “low modes”. They correspond to modes with the lowest value of $|\vec{k}|_{\delta}$ defined by (2.2). Any mode $\hat{\omega}_{\vec{k}}$ with $|\vec{k}|>\mbox{max}\{1,\delta^{2}\}$ will from here on be referred to as a “high mode”.

To measure the relative energy in the low modes, we define the stochastic order parameter,

[TABLE]

where $\hat{\omega}_{(1,0)}$ and $\hat{\omega}_{(0,1)}$ solve (2.1). Due to the condition, $\hat{\omega}_{(k_{1},k_{2})}=\bar{\hat{\omega}}_{(-k_{1},-k_{2})}$, the relative energy in all of the low modes can be captured by $Z_{vort}(t)$. The value of $Z_{vort}(t)$, bounded between 0 and 1, corresponds to the proximity of the solution to an $x$-bar, $y$-bar or dipole state. If the dynamics drive $Z_{vort}(t)$ to increase to 1, there is more energy in $\hat{\omega}_{(1,0)}$ relative to $\hat{\omega}_{(0,1)}$, indicating the system is in an $x$-bar state. Conversely if $Z_{vort}(t)$ falls toward 0, the system would be observed to be in a $y$-bar state. If $Z_{vort}(t)$ instead stays near 1/2, the system is in a dipole state with relative energy in the low modes comparable in magnitude.

The finite dimensional system that we will use to model (2.1) will be defined in terms of the lowest eight Fourier modes, which for notational convenience we denote as

[TABLE]

The variables $\omega_{1,2,3,4}$ correspond to the low modes, while $\omega_{5,6,7,8}$ represent the role of all the high modes. Since the solution $\omega(x,y)$ of (2.1) is real valued, the following complex conjugacy relationship must still hold,

[TABLE]

Thus the reduced model will be an eight dimensional approximation to the dynamics of (2.1). To derive the model, we apply a center manifold reduction to (2.1) with $\sigma_{\vec{k}}=0$ for all $\vec{k}$ to obtain an eight-dimensional deterministic ODE, which is the model studied in [1], and then add noise back to that system to obtain the final eight-dimensional SDE model we study here.

To carry out the center manifold reduction onto the lowest eight modes, assume for $\hat{\omega}_{\vec{k}}$ with $\vec{k}\notin\{(\pm 1,0),(0,\pm 1),(\pm 1,\pm 1)\}=:\mathcal{K}_{0}$, that there exists a smooth function $H(\omega_{1},\dots,\omega_{8};\vec{k})$ such that the eight-dimensional manifold defined by

[TABLE]

is invariant for the deterministic dynamics of (2.1) with $\sigma_{\vec{k}}=0$ for all $\vec{k}$. We refer to this as a center manifold because it is defined in terms of the lowest eight modes, which have the weakest linear decay rates. Based on this assumption, one can then in principle compute the coefficients of the Taylor expansion of $H(\cdot,\vec{k})$ to any order for each $\vec{k}$ by taking the derivative of each of the low modes in two ways (via the function $H$ and (2.1) with $\sigma_{\vec{k}=0}$) and equating coefficients. See [1] for the details of the derivation.

The reduction is local and will only be valid in a size $\mathcal{O}(\nu)$ neighborhood of zero due to the small spectral gaps for the operator $\nu\Delta$. Additionally, while the existence of a finite dimensional (inertial) model of the system (2.1) that describes the global dynamics cannot be expected [21], the model still provides meaningful insight into the role $\delta$ plays in selecting the dominant quasi-stationary state for small initial conditions. For additional examples in which similar reductions of the Navier-Stokes equation to a finite dimensional model have been used to understand global dynamics see [7, 15].

Adding independent (real) Brownian motions $W_{1,3,5,7}$ to each equation of the resulting ODE model leads to our final SDE model

[TABLE]

Note that (2.6) with $\sigma_{1,3,5,7}=0$ corresponds to the ODE model derived in [1]. To compare the dynamics of this model to that of $Z_{vort}(t)$, defined in (2.3), we define the analogous order parameter for the SDE model,

[TABLE]

which again is used to determine towards which quasi-stationary state the system trends. Here, $\omega_{1}(t)$ and $\omega_{3}(t)$ are solutions to the reduced system (2.6). The Monte Carlo simulation of the reduced model finds that the dominant quasi-stationary state depends on the aspect ratio of $D_{\delta}$ in the same way as the deterministic model, studied in detail in [1]. The main result there, which describes the selection of quasi-stationary states in (2.6) with $\sigma_{1,3,5,7}=0$, can be described by the following theorem.

Theorem 2.1**.**

[1, Theorem 3.4]** For $\delta\in\left(\sqrt{\frac{2}{3}},\sqrt{\frac{3}{2}}\right)$, under the dynamics of (2.6) with $\sigma_{1,3,5,7}=0$, if $\delta>1$, then $Z_{red}(t)\rightarrow 1$, indicating evolution to an $x$-bar state. Conversely if $\delta<1$, then $Z_{red}(t)\rightarrow 0$, indicating evolution to a $y$-bar state. For $\delta=1$, there exists a one-dimensional center manifold of fixed points in the phase space that determines the asymptotic limit of $Z_{red}(t)$. This center manifold is foliated with co-dimension one stable manifolds in which solutions converge to the corresponding fixed point. Exactly one of these manifolds corresponds to each of the limits $Z_{red}(t)\rightarrow 1$ and $Z_{red}(t)\rightarrow 0$. Thus, generic initial conditions are seen to evolve to the dipole state.

Remark 2.2**.**

The order parameter considered in [1] was instead the ratio $R(t)=|\omega_{1}(t)|^{2}/|\omega_{3}(t)|^{2}$. Theorem 2.1 frames the result in terms of the order parameter $Z_{red}(t)$. The choice to now consider $Z_{red}(t)$ is for convenience with regards to numerical simulation due to its being bounded between 0 and 1.

Remark 2.3**.**

A straightforward computation shows that, for any $\delta$, the set \{\mathop{\mathrm{missing}}{Im}\nolimits(\omega_{1})=\mathop{\mathrm{missing}}{Im}\nolimits(\omega_{3})=\mathop{\mathrm{missing}}{Im}\nolimits(\omega_{5})=\mathop{\mathrm{missing}}{Im}\nolimits(\omega_{7})=0\} is invariant under the dynamics of (2.6) with $\sigma_{1,3,5,7}=0$. Since the real subsystem is invariant in the deterministic setting, we simulate the reduced model where the modes, $\omega_{1,3,5,7}$, as well as the Wiener processes, $W_{1,3,5,7}$, are all real valued.

3 Numerical simulation of the vorticity equation and reduced model

This section provides simulations of the vorticity equation (2.1) and of the reduced model (2.6). Via Monte Carlo simulation, the average evolution of the order parameters $Z_{vort}(t)$ and $Z_{red}(t)$ will be plotted for several values of $\delta$ near 1. It will be seen that the reduced model captures the selection of the quasi-stationary states via the parameter $\delta$. In particular, in both models, for a particular value of $\delta\approx 1$, the system’s selection of its dominant quasi-stationary state is consistent with the motivating results, given by Theorem 2.1. In particular, the system selects, as the dominant quasi-stationary state, a dipole for $\delta=1$, an $x$-bar for $\delta>1$, and a $y$-bar for $\delta<1$.

The simulation of (2.1) is done via a spectral method which includes Fourier modes $\hat{\omega}_{\vec{k}}$ with $\vec{k}\in\mathcal{K}:=\{\vec{k}=(k_{1},k_{2})\in\mathbb{Z}^{2}:0\leq|k_{1}|,|k_{2}|\leq 64\mbox{ and }(k_{1},k_{2})\neq(0,0)\}$; see [13]. A condition of exponential decay is imposed on the noise coefficients $\sigma_{\vec{k}}$ seen in (1.4),

[TABLE]

Similar to [6], simulations are conducted with $\sum_{\{\vec{k}\in\mathcal{K}\}}e^{-\alpha_{0}|\vec{K}|^{2}}=1$. For our set $\mathcal{K}$, this means $\alpha_{0}\approx 0.349$. Time was finely discretized and a tamed semi-implicit Euler-Maruyama method was implemented to simulate the stochastically forced reduced system.

We verify the selection of the dominant quasi-stationary state using Monte Carlo simulation where the average path over $N$ trials is plotted. Individual runs will be denoted by $Z^{i}_{vort}(t)$ and $Z^{i}_{red}(t)$, for $i=1,\dots N$, with corresponding averages given by

[TABLE]

Similarly, we define the empirical variances to be

[TABLE]

It will also be useful to plot the time averages of these Monte Carlo averages. To produce a meaningful average we introduce a “burn-in time”, $t_{burn}$, and ignore the initial period during which $\bar{Z}_{vort}(t)$ and $\bar{Z}_{red}(t)$ have not yet stabilized. Define this time average for any function $f(t)$ defined on $t_{burn}\leq t\leq T$ to be

[TABLE]

3.1 Vorticity Equation

Plotted in Figures 2-4 are $\bar{Z}_{vort}(t)$, the time average $A(\bar{Z}_{vort},t_{burn})$, and the 95% confidence intervals defined via

[TABLE]

Also included are average contour plots for the vorticity. We use $N=200$ and for each trial use zero initial conditions and $\nu=0.001$. For $\delta=1$, Figure 2(a) shows $\bar{Z}_{vort}(t)$ remains near 1/2 for the duration of the simulation. We use a burn-in time of $t_{burn}=0$ when computing the time average since on the symmetric domain it is clear there is no transient initial period. In Figure 2(b), the average contour plot for each individual trial are themselves averaged over the $N=200$ trials, reflecting a dipole.

The simulations exhibited in Figures 3(a) and 3(b) show that, for $\delta=1.1$, the order parameter increases initially and the average contour plot looks like that of an $x$-bar state. In Figure 3(a), $t_{burn}=100$ is used when computing the time average.

Lastly for $\delta<1$ the simulations exhibited in Figures 4(a) and 4(b) show that, for $\delta=0.9$, the order parameter decreases over an initial period of time and the average contour plot looks like that of a $y$-bar state. Here we again set $t_{burn}=100$.

Provided in Figure 5 are plots of $\bar{Z}_{vort}(t)$ for $\delta=1.10$, $\delta=1.0$ and $\delta=0.90$ averaged over $N=1000$ trials. This is to show that as the number of trials increase, the variance is decreasing without changing the mean behavior. The variances all remain generally between 0.06-0.08. While the variance does decrease, the limiting value of $\bar{Z}_{vort}(t)$ remains relatively unchanged compared to what is seen when averaging over $N=200$ trials.

For completeness, see also the discussion in §5, we also include in Figure 6 a simulation with $\nu=0.001$ that represents a single sample path for $\delta=1.04$, the value of $\delta$ for which transitions among the quasi-stationary states were observed in [6].

Figure 6(a) shows that individual trajectories exhibit transitions between quasi-stationary states, visiting the dipole and both bar states, as also observed in [6] for the same value of $\delta$. However Figure 6(b) shows that $\mathbb{E}[Z_{vort}(t)]$ picks the dominant state. We now compute the time average of a randomly selected individual trial, given by $A(Z_{vort},t_{burn})$, to confirm that it tracks the Monte Carlo average, $\bar{Z}_{vort}(t)$. Figure 7 shows two things. First, for the given values of $\delta$, a sample path may experience many transitions among the quasi-stationary states. Second, the time average of the sample path does eventually track the Monte Carlo average.

3.2 Reduced Model

We now turn our attention to the reduced model (2.6). We confirm numerically that the reduced model captures the qualitative dynamics of the full vorticity equation with regard to the dominant quasi-stationary state.

As we mentioned in Remark 2.3 we will be working with the real system in which $\omega_{1,3,5,7}$, as well as the Wiener processes, $W_{1,3,5,7}$, are all real valued. This leads to the following system which serves as the acting reduced model in the upcoming simulations.

[TABLE]

To be consistent with the spatial decay of the noise in the simulations of the stochastically forced vorticity equation (2.1), given by (3.1), we choose

[TABLE]

First we aim to establish that the reduced model (3.2) can serve as a good approximation to the vorticity equation with noise, (2.1), for $\delta\approx 1$. Second, it will be established that the selection of the bar or dipole state that dominates is consistent with the results of [1] for the deterministic equation: $x$-bar for $\delta>1$, $y$-bar for $\delta<1$, and dipole for $\delta=1$.

Figure 8 shows numerical evidence supporting that the dynamics of the order parameter, governed by the reduced system (3.2), follows the same trend as when the full vorticity equation is simulated.

The plots of these Monte Carlo simulations (averaged over $N=200$ trials) show that the trend toward the appropriate quasi-stationary state is captured by the reduced model. Starting with zero initial conditions, when the noise is added, the simulations show that for $\delta>1$, the order parameter increases toward 1, indicating evolution to an $x$-bar state. Conversely, for $\delta<1$, the order parameter decreases toward a value corresponding to a $y$-bar state. Finally, when $\delta=1$, $\bar{Z}_{red}(t)$ remains near 1/2 indicating the system is in a dipole state. Figures 9(a)-9(e) serve to compare the evolution of $\bar{Z}_{vort}(t)$ and $\bar{Z}_{red}(t)$taken over $N=200$ trials, for values of $\delta$ close to 1. The bars denote the error for the 95% confidence intervals for $\bar{Z}_{vort}$ (bold) and $\bar{Z}_{red}$ (thin).

One can see that $\bar{Z}_{vort}(t)$ and $\bar{Z}_{red}(t)$ both trend in the same direction, with similar variances (typically between 0.06-0.08 for $0\leq t\leq 2000$). Furthermore, their respective confidence intervals begin to converge until they overlap. Indeed the model can be used to determine towards which quasi-stationary state the system evolves for a given value of $\delta.$

4 Perturbation Analysis

Motivated by the numerics from §3, this section investigates the expected behavior of $Z_{red}(t)$ as $\delta\rightarrow 1$ while viewing the problem as a perturbation from the $\delta=1$ and $\nu=0$ case. Using the backward Kolmogorov equation associated to (3.2), the goal is to derive a system of PDE that will provide insight on how the expected value of $Z_{red}(t)$, to leading order, depends on values of $\delta$ close to 1. To do this we pose the problem as a perturbation of the spatial domain, setting $\delta^{2}=1+\epsilon_{0}\epsilon$. Here, $0<\epsilon\ll 1$ acts as the small perturbation parameter and $\epsilon_{0}=\pm 1$ determines which dimension of the torus is longer. Following known homogenization techniques, see for example [18], we scale (3.2) in a way that reveals a slow-fast system of SDE. Then, we write the backward Kolmogorov equation to reach the ultimate goal of determining equations that govern the limiting evolution of $\mathbb{E}[Z_{red}(t)]$ as $\epsilon\rightarrow 0$ once the fast variables are averaged out. First, for ease of notation, rename the dependent variables as follows,

[TABLE]

Now (3.2) can be expressed as

[TABLE]

Before inserting the Taylor expansions in $\epsilon$ for the coefficients with $\delta^{2}=1+\epsilon_{0}\epsilon$, we first scale (4.1) appropriately to obtain a clear slow-fast system. As in [1], the low modes represented by $\tilde{p}$ and $\tilde{q}$ correspond to the slow variables while the high modes, $\tilde{r}$ and $\tilde{s}$, represent the fast variables. Below, we give a more general version of the scaled equations for just the $\tilde{p}$ (analogous to $\tilde{q}$) and $\tilde{r}$ (analogous to $\tilde{s}$) equations. We use the following space-time and parameter scalings: $\nu=\epsilon^{\mu}\nu_{0}$, $\tilde{p}=\epsilon^{\xi}p$, $\tilde{q}=\epsilon^{\xi}q$, $\tilde{r}=\epsilon^{\eta}r$, $\tilde{s}=\epsilon^{\eta}s$, and $\tau=\epsilon^{\gamma}t$. To simplify the scaled equations, we will relate $\mu$, $\xi$, $\eta$ and $\gamma$ to put the resulting system in a more desirable form. We neglect the $\epsilon$ dependence of $p$, $q$, $r$ and $s$ for readability. Below, the “prime” notation denotes differentiation with respect to the scaled time variable, $\tau$.

[TABLE]

Now set $2\eta=\mu+\gamma$, $2\xi=\mu-\gamma\Rightarrow\gamma=\mu-2\xi$, with $0<\gamma<\xi<\frac{\mu}{2}<\eta<\mu.$ Then the fully scaled system (still neglecting Taylor expansions of $\delta$ in $\epsilon$ for now) becomes

[TABLE]

For the scaled SDE (4.2), let $b^{\epsilon}=(b_{p}^{\epsilon},b_{q}^{\epsilon},b_{r}^{\epsilon},b_{s}^{\epsilon})$ denote the drift vector and $\Sigma^{\epsilon}(p,q,r,s;\sigma_{1},\sigma_{3},\sigma_{5},\sigma_{7})$ denote the diffusion matrix so that (4.2) can be expressed, for $\mathcal{X}^{\epsilon}=(p,q,r,s)$ and $\frac{d\mathcal{W}}{d\tau}=\left(\frac{dW_{1}}{d\tau},\frac{dW_{3}}{d\tau},\frac{dW_{5}}{d\tau},\frac{dW_{7}}{d\tau}\right)$, as

[TABLE]

Now replacing the $\delta$ coefficients appearing in (4.2) with their Taylor expansions for $\delta^{2}=1+\epsilon_{0}\epsilon$ up to $\mathcal{O}(\epsilon^{3})$, the drift vector is given by (still supressing $\epsilon$ dependence of $p,q,r,s$),

[TABLE]

and the diffusion matrix by

[TABLE]

With $\mbox{H}(u^{\epsilon})$ denoting the Hessian matrix of $u^{\epsilon}$, we now write the backward Kolmogorov equation for (4.3), which is defined as

[TABLE]

The backward Kolmogorov equation has the useful property that the evolution of $u^{\epsilon}(\mathcal{X}^{\epsilon},\tau)$ gives

[TABLE]

Thus one ultimately is interested in initializing (4.10) with $\phi(p,q)=Z_{red}=\frac{p^{2}}{p^{2}+q^{2}}$, but for now we proceed with a general initial condition, $\phi$. We seek a solution to (4.10) that takes the form

[TABLE]

and wish to find the limiting dynamics, $u^{\epsilon}$ as $\epsilon\to 0$. As the goal is to identify the leading order expansion for $u^{\epsilon}$ we determine a system of PDEs for $u_{0}$, $u_{1}$, and $u_{2}$. We present now the calculations that lead to the characterization of $u_{0}$, $u_{1}$, and $u_{2}$, see (4.20).

Define $\mathcal{L}^{\epsilon}$ to be the operator acting on the right hand side of (4.10), so that $\frac{\partial u^{\epsilon}}{\partial\tau}=\mathcal{L}^{\epsilon}u^{\epsilon}$. Decomposing $\mathcal{L}^{\epsilon}$ by powers of $\epsilon$ using the expressions for $b^{\epsilon}$ and $\Sigma^{\epsilon}$ given in (4.4) and (4.9), we write

[TABLE]

A select few of the operators, $\mathcal{L}_{i}$, $i=1,\dots,19$, that are most important in computing the leading order equations is provided in (4.12). A complete list of the expressions of the 19 operators can be found in the appendix.

[TABLE]

Now we choose explicit values for $\xi$ and $\eta$ to obtain a simpler, but still representative system: $\eta=2$, $\xi=1$, hence $\mu=3$ ($\nu=\epsilon^{\mu}\nu_{0}=\epsilon^{3}\nu_{0}$). Note that this choice of parameters is not unique. Our goal is to make a choice that makes the computations that follow tractable, preserves a clear slow-fast system in the scaled equations (4.2) and results in a system that exhibits qualitatively the same behavior as the original system in terms of the selection mechanism to $x$-bar, $y$-bar and dipole states depending on the values of $\delta$. With this choice, the generator $\mathcal{L}^{\epsilon}$ of the PDE becomes

[TABLE]

The ansatz given in (4.11) can now be inserted into the backward Kolmogorov equation (4.10) using the expression of $\mathcal{L}^{\epsilon}$ given above in (4). Matching coefficients on both sides of the equation yields the following leading order equations,

[TABLE]

Equation (4.14a) implies $u_{0}$ lies in the kernel of $\mathcal{L}_{0}$, which elliptic PDE theory tells us contains only functions constant in $r$ and $s$. Since $\mathcal{L}_{1}$ is also a differential operator in $r$ and $s$ only, (4.14b) implies that $u_{1}$ is constant in $r$ and $s$ as well. One can also see that $u_{0}$ and $u_{1}$ are in the kernel of each of $\mathcal{L}_{2,3,16,17}$ (see appendix). Hence the leading order system given by (4.14a)-(4.14d) can be reduced to

[TABLE]

where $\mathcal{L}_{0}$, $\mathcal{L}_{1}$, $\mathcal{L}_{4}$, $\mathcal{L}_{5}$, and $\mathcal{L}_{8}$ are presented in (4.12). Let $\rho^{\infty}(r,s;p,q)$ be the stationary density that satisfies the adjoint problem

[TABLE]

Once $\rho^{\infty}$ is known we can integrate against the invariant measure to obtain the solvability conditions for equations (4.15c) and (4.15d)

[TABLE]

Before we consider the integrals in (4.16), $\rho^{\infty}$ must be identified. The operator,

[TABLE]

corresponds to the backward Kolmogorov equation for the following system, parameterized by the fixed (slow) variables $p$ and $q$.

[TABLE]

These processes are independent Ornstein-Uhlenbeck processes and are therefore Gaussian. The equilibrium (stationary) density which corresponds to $\rho^{\infty}(r,s;p,q)$ is that of the bivariate Gaussian distribution with

[TABLE]

Therefore the invariant joint density is

[TABLE]

To aid in the computations of the integrals given in (4.16), the following integral evaluations will be useful and can be simply obtained through the mean and variance of the stationary distribution.

[TABLE]

Next consider the solvability conditions (4.16) one at a time. From the $u_{0}$ equation and the integral (4.17a),

[TABLE]

From this we obtain the effective equations for $p$ and $q$ for small $\epsilon$ after the fast variables $r$ and $s$ are averaged out. The slow motion can hence be approximated, for $0<\epsilon\ll 1$, by $\bar{p}$ and $\bar{q}$ governed by,

[TABLE]

Since $\epsilon_{0}$ dependence does not appear in the first order equations, we will need to determine $u_{1}$ to see its effects. Consider the solvability condition for $u_{1}$ in (4.16). Computing this integral requires us to evaluate the following integrals.

[TABLE]

In evaluating these, we see

[TABLE]

Since $u_{2}$ depends on the fast variables, this final integral cannot yet be computed. It will eventually be handled numerically. Thus, formally, we have $u^{\epsilon}=u_{0}+\epsilon u_{1}+\epsilon^{2}u_{2}+\dots$ satisfying,

[TABLE]

We shall consider the system (4.20) together with the initial conditions

[TABLE]

The PDE for $u_{0}$ immediately stands out as the backward Kolmogorov equation corresponding to the system given in (4.18). Despite its simple looking form, the regularity at the origin of the coefficients on the first derivative terms turn out to be a borderline case with regards to well posedness, see for example Chapter III, Section 1 of [11]. Nevertheless, we proceed formally and solve for $u_{0}$, $u_{1}$, and $u_{2}$ numerically after providing the initial condition

[TABLE]

so that $u^{\epsilon}(p,q,r,s,\tau)=\mathbb{E}[Z_{red}(\tau)|p_{0}=p,\mbox{ }q_{0}=q,\mbox{ }r_{0}=r,\mbox{ }s_{0}=s]$. The simulations of the system (4.20) provide an approximation to the deterministic evolution of $\mathbb{E}[Z_{red}(\tau)]$ after averaging out the fast motion. The simulations of the system (4.20) provided in this section were conducted via finite differences on the domain $(p,q,r,s,\tau)\in[-5,5]^{4}\times[0,T]$ with Neumann boundary conditions.

When $\epsilon=0$, which implies $\delta=1$, the system is unperturbed and $u_{0}(p,q,t)=\mathbb{E}[Z_{red}(\tau)|p_{0}=p,\mbox{ }q_{0}=q]$. In this case, the results of the simulation show that the expected value of the order parameter $Z_{red}(\tau)$ converges to 1/2 for any initial values $p_{0}=p$ and $q_{0}=q$ of (2.6), independent of $r$ and $s$. This indicates that the unperturbed system evolves to a dipole state, even if the initial state is close to an $x$- or $y$-bar state. Figure 10 illustrates the evolution to a dipole for $\epsilon=0$ for several initial conditions $(p,q)$ chosen within the domain.

We conclude this section by comparing the numerical approximation to $\mathbb{E}[Z_{red}(\tau)]$, given by $u^{\epsilon}(\tau)$ to the average path of the order parameter $Z_{red}(t)$, i.e. $\bar{Z}_{red}(t)$, as computed via Monte Carlo simulation in §3. Since the two models evolve on different timescales, we rescale $\tau$ so that our averaged PDE model is evolving on the original timescale. As such, suppressing the spatial $(p,q,r,s)$ dependency, let $\hat{u}^{\epsilon}(t):=u_{0}(t/\epsilon)+\epsilon u_{1}(t/\epsilon)=u_{0}(\tau)+\epsilon u_{1}(\tau)$ denote the $\mathcal{O}(\epsilon)$ approximation to $u^{\epsilon}$ in the original timescale, and let $\bar{Z}_{red}(t)$ denote the Monte Carlo average path of the order parameter under the dynamics of (2.6) obtained via the Monte Carlo simulations described in §3.

Since the Monte Carlo simulation of (2.6) used zero initial conditions, we plot the approximation $\hat{u}^{\epsilon}(p,q,t)$ for $(p,q)=(0.1,0.1)$, close to the origin. With this choice, figures 11(a) and 11(b) show that in the perturbed system, the $\mathcal{O}(\epsilon)$ approximation to $u^{\epsilon}(\tau)=\mathbb{E}[Z_{red}(\tau)]$ evolves toward 0 or 1 depending on the sign of $\hat{\epsilon}:=\epsilon_{0}\epsilon.$

Using these simulations with $\nu$ scaled, we now explore how the intervals on which $\hat{u}^{\epsilon}$ serves as a good approximation to $\bar{Z}_{red}(t)$ depend on the perturbation parameter. Figure 12 shows the relative error (RE) given by

[TABLE]

On some initial interval of time, the PDE approximation, $\hat{u}^{\epsilon}$ indeed serves as a close approximation to $\bar{Z}_{red}(t)$. Furthermore, as expected, the smaller the perturbation parameter $\epsilon$, the longer the approximation is valid.

5 Concluding remarks and future directions

In this paper, we developed a finite dimensional SDE model that can be used to elucidate the dynamics of the 2D Navier-Stokes vorticity equation with noise. Monte Carlo simulation of the reduced model showed that the major qualatative property of the system, i.e. the dominant quasi-stationary state, can be determined from the model. In particular, as has been observed numerically and rigorously, the existence and attracting nature of these quasi-stationary states play an important role in the evolution of the stochastic Navier-Stokes vorticity equation. Specifically, the aspect ratio of the periodic domain, $D_{\delta}=[0,2\pi\delta]\times[0,2\pi]$, determines whether generic solutions evolve toward an $x$-bar state ($\delta>1$), a $y$-bar state ($\delta<1)$, or a dipole state ($\delta=1$).

Perturbation analysis then shows that the proposed reduced model can be viewed as a slow-fast system, Subsequent averaging and homogenization methods show the leading order behavior as the perturbation parameter $\delta\approx 1$ goes to $\delta=1,$ in relation to how the viscocity parameter $\nu$ vanishes.

The numerical studies in §3 show that, on average, the system prefers to trend toward the appropriate quasi-stationary state as determined by $\delta$, see Figure 6(b). However, one can see from the sample path plotted in Figure 6(a), individual sample paths do exhibit transitions between $x$-bar and $y$-bar states, as it has also been observed in [6].

In regards to future directions, there are a number of interesting questions that one can ask and hope to answer. To begin with, the perturbation analysis of §4 is formal and one would like to prove both well-posedness of (4.20) and validity of the perturbation expansion.

In addition, the numerical studies of §3 suggest that while there are transitions at the individual sample path level, the system tends to converge to the preferred state depending on whether $\delta<1$ or $\delta>1$. One would like to make this mathematically rigorous. Furthermore, one could potentially use the reduced model of §2 to build a related large deviations theory describing probabilities of the system being in, and exit times for leaving, one of the quasi-stationary states.

Furthermore, the form of the noise used in (1.4) and (4.1), with scaling factor $\sqrt{\nu}$, is common in the literature, see for example [16] and the references therein. One of the important questions in the literature is the investigation of the convergence of the corresponding invariant measure as $\nu\to 0$ to that of the 2D Euler equation, see [10]. The support of the limiting invariant measure (i.e. as $\nu\to 0$), is in general still an open question, see [16], and it has been resolved in special cases in [3] and [4]. In our work, we study the long time behavior of the vorticity equation as $\nu\to 0$ and in particular the selection mechanism for the dominant quasi-stationary states. It is reasonable to expect that there is a relation between the selection mechanism and the support of the limiting invariant measure. This is an intriguing question that is left for future work and it is beyond the scope of this work.

6 Appendix

The complete list of operators in the Kolmogorov equation (4.10) is given by

[TABLE]

[TABLE]