Selection of quasi-stationary states in the stochastically forced Navier-Stokes equation on the torus
Margaret Beck, Eric Cooper, Gabriel Lord, Konstantinos Spiliopoulos
SSelection of quasi-stationary states in the stochastically forcedNavier-Stokes equation on the torus
Margaret Beck ∗ , Eric Cooper ∗ , Gabriel Lord † , Konstantinos Spiliopoulos ∗ Abstract
The stochastically forced vorticity equation associated with the two dimensional incompressibleNavier-Stokes equation on D δ := [0 , πδ ] × [0 , π ] is considered for δ ≈
1, periodic boundary conditions,and viscocity 0 < ν (cid:28)
1. An explicit family of quasi-stationary states of the deterministic vorticityequation is known to play an important role in the long-time evolution of solutions both in the presenceof and without noise. Recent results show the parameter δ plays a central role in selecting which of thequasi-stationary states is most important. In this paper, we aim to develop a finite dimensional modelthat captures this selection mechanism for the stochastic vorticity equation. This is done by projectingthe vorticity equation in Fourier space onto a center manifold corresponding to the lowest eight Fouriermodes. Through Monte Carlo simulation, the vorticity equation and the model are shown to be in agree-ment regarding key aspects of the long-time dynamics. Following this comparison, perturbation analysisis performed on the model via averaging and homogenization techniques to determine the leading orderdynamics for statistics of interest for δ ≈ Consider the 2D incompressible Navier-Stokes equation, ∂ u ∂t = ν ∆ u − ( u · ∇ ) u − ∇ p ∇ · u = 0 , (1.1)on the possibly asymmetric torus ( x, y ) ∈ D δ := [0 , πδ ] × [0 , π ] with δ ≈
1, periodic boundary conditions,and viscosity 0 < ν (cid:28)
1. To obtain the equivalent vorticity formulation of the equation, take the curl ofthe vector field u and set ω = (0 , , · ( ∇ × u ) to find ∂ t ω = ν ∆ ω − u · ∇ ω, u = (cid:32) ∂ y ( − ∆ − ) − ∂ x ( − ∆ − ) (cid:33) ω. (1.2)The relation between u and ω is known as the Biot-Savart law. The periodic boundary conditionsinsure (cid:82) D δ ω = 0, and therefore ∆ − ω is well-defined. ∗ Department of Mathematics and Statistics, Boston University, Boston, MA 02215. E-mail: [email protected],[email protected], [email protected]. M.B. was partially supported by National Science Foundation (NSF) grant DMS1411460 and K.S. was partially supported by NSF DMS 1550918. † Department of Mathematics, Heriot-Watt University, Edinburgh, EH14 4AS, UK. Email: [email protected] a r X i v : . [ m a t h . D S ] J a n dding random forcing to the system allows to account for stochasticity/genericity in the system, seefor example [17, 5, 16]. In particular, we add a stochastic forcing term to (1.2) to obtain the stochastic 2Dvorticity equation, ∂ t ω = ν ∆ ω − u · ∇ ω + ∂ W ∂t , u = (cid:32) ∂ y ( − ∆ − ) − ∂ x ( − ∆ − ) (cid:33) ω. (1.3)The noise is white in time, colored in space, and takes the form, for (cid:126)k = ( k , k ) (cid:54) = (0 , W ( t, x, y ) = √ ν (cid:88) (cid:126)k ∈K⊂ Z \{ (0 , } σ (cid:126)k e i( k x/δ + k y ) β (cid:126)k ( t ) , (1.4)with spatial correlation σ (cid:126)k and K to be commented on below. Here β ( t ) = { β (cid:126)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 vorticityremains real valued for all times t ≥
0, the following complex conjugacy conditions are imposed, ¯ σ (cid:126)k = σ − (cid:126)k and ¯ β (cid:126)k = β − (cid:126)k . Additional assumptions are often placed on the noise coefficients, σ (cid:126)k , to insure certainsmoothness properties of solutions. In particular, we assume that there exist fixed positive constants C and α such that | σ (cid:126)k | ≤ C e − α | (cid:126)k | so that solutions will then be analytic in space [14]. Since the boundaryconditions force solutions of the deterministic equation to satisfy (cid:82) D δ ω = 0, we choose σ (0 , = 0 so thisproperty is preserved. Note that if σ (cid:126)k = 0 for all (cid:126)k ∈ Z then (1.3) reduces to the deterministic vorticityequation.Although an L 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 correspondto 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, ω ( x, y, t ) = e − νδ t [ a cos( x/δ ) + a sin( x/δ )] + e − νt [ a cos( y ) + a sin( y )] . (1.5)If δ = 1, then any member of this family is an exact solution to the deterministic vorticity equation.If δ (cid:54) = 1, then (1.5) remains a solution if and only if a = a = 0 or if a = a = 0. These members, whichonly depend on one spatial variable, are called bar states, and they are also known as unidirectional orKolmogorov flow. The x - and y -bar states are members of this family given by ω xbar ( x, t ) = e − νδ t sin( x/δ ) , ω ybar ( y, t ) = e − νt sin y, or similarly with sine replaced by cosine. The associated velocity fields are given by u xbar ( x, t ) = − δe − νδ t (cid:32) x/δ ) (cid:33) , u ybar ( y, t ) = e − νt (cid:32) cos y (cid:33) , respectively. The dipoles are also members of the family (1.5) and are given by ω dipole ( x, y, t ) = e − νδ t sin( x/δ ) + e − νt sin y, or similarly with sine replaced by cosine, with velocity field u dipole ( x, y, t ) = (cid:32) e − νt cos y − δe − νδ t cos( x/δ ) (cid:33) . t = 0 on the symmetric torus( δ = 1) are shown in Figure 1. (a) x-bar: ω xbar = sin( x ) (b) y-bar: ω ybar = sin( y ) (c) Dipole: ω dipole = sin( x ) + sin( y ) Figure 1: Contour plots of the three quasi-stationary states on the symmetric torusWhen ν = 0, equation (1.1) becomes the Euler equation. It is reasonable to expect that stationarysolutions of the Euler equation could play an critical role in the evolution of the Navier-Stokes equation for0 < ν (cid:28)
1. However, there are infinitely many stationary solutions, including the bars and dipoles, and soit is not immediately clear how to determine which would be most important. In [20], entropy argumentsand extensive numerical studies were conducted in the case δ = 1 and suggested that the bars and dipolesshould be the two most important stationary solutions of the Euler equations. Although both states wereobserved after initial transient periods in the evolution of the Navier-Stokes equation, interestingly thedipole seemed to emerge for a large class of initial data, whereas the bar states only emerged for a specialclass of initial data. Subsequent work, again for the deterministic system, showed that indeed the barstates attract nearby solutions at a rate much faster than the background global decay rate, confirmingtheir importance as quasi-stationary states. Results in the case δ = 1 can be found in [2, 9] and resultsfor more general values of δ 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, withthe dipole being dominant for δ = 1 and the bar states being dominant for δ (cid:54) = 1.In this paper we develop a low dimensional model that captures how the dominant quasi-stationarystate in the stochastically forced Navier-Stokes equation is selected by the aspect ratio of the spatialdomain, δ . 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 tosuggest that δ could provide such a selection mechanism. The former paper [1] was our previous workfocusing on the deterministic vorticity equation, (1.2), in which we derived a finite-dimensional model thatcaptured the selection mechanism via the parameter δ . We now seek to use that same finite-dimensionalmodel, but with the addition of noise, to numerically investigate the selection mechanism for the stochasticequation (1.3). Indeed, one can see from the sample paths of Figure 6a, that individual sample pathsexhibit 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 § §
3, to determine the validity of the SDE model, we comparestatistics 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, δ > δ < δ = 1. In particular, solutions to both systems evolve towards an x -bar, y -bar, and dipole in the threerespective cases. In §
4, we further examine the SDE model by viewing it as a perturbation in the limitas | δ − | and ν converge to zero. We show that, after appropriate time-space rescalings, the system canbe viewed as a slow-fast system and classical averaging and homogenization techniques apply. Via thebackward Kolmogorov equation, a system of PDEs that governs the leading order dynamics of a key orderparameter, E [ Z red ( t )], defined in (2.7), is derived. This gives us an additional formal approximation tothe expected value of the order parameter, which we can use to show the selection of the quasi-stationarystate. Numerically solving the PDEs allows us to approximate the evolution of E [ Z red ( t )] for values of δ close to 1, at least on some initial finite interval of time. Conclusions and future directions then follow in § Due to the form of the family of solutions (1.5), it is most convenient to express the stochastic vorticityequation in Fourier space. Hence, letting ω ( x, y ) = (cid:88) (cid:126)k (cid:54) =(0 , ˆ ω (cid:126)k e i( k x/δ + k y ) , ˆ ω (cid:126)k = 14 π δ (cid:90) D δ ω ( x, y ) e − i( k x/δ + k y ) d x d y, we obtain, for (cid:126)j , (cid:126)k and (cid:126)l (cid:54) = (0 , ω (cid:126)k = − νδ | (cid:126)k | δ ˆ ω (cid:126)k − δ (cid:88) (cid:126)l (cid:104) (cid:126)k ⊥ ,(cid:126)l (cid:105)| (cid:126)l | δ ˆ ω (cid:126)k − (cid:126)l ˆ ω (cid:126)l + √ νσ (cid:126)k ˙ β (cid:126)k = − νδ | (cid:126)k | δ ˆ ω (cid:126)k − δ (cid:88) (cid:126)j + (cid:126)l = (cid:126)k (cid:104) (cid:126)j ⊥ ,(cid:126)l (cid:105) (cid:32) | (cid:126)l | δ − | (cid:126)j | δ (cid:33) ˆ ω (cid:126)j ˆ ω (cid:126)l + √ νσ (cid:126)k ˙ β (cid:126)k , (2.1)where | (cid:126)k | δ = k + δ k , (cid:126)k ⊥ = ( k , − k ) . (2.2)Viewing the system in Fourier space allows us to use the relative energy in certain modes to mea-sure the proximity of solutions to an x -bar, y -bar, or dipole state. The x -bar states, e − νδ t cos( x/δ ) and e − νδ t sin( x/δ ), correspond to solutions with energy only in the (cid:126)k = ( ± ,
0) modes and the y -bar states, e − νt cos( y ) and e − νt sin( y ), correspond to solutions with energy only in the (cid:126)k = (0 , ±
1) modes. Solutionswith energy in both the (cid:126)k = ( ± ,
0) and (cid:126)k = (0 , ±
1) modes correspond to the dipole state. These fourmodes are the lowest in the system and will be referred to as the “low modes”. They correspond to modeswith the lowest value of | (cid:126)k | δ defined by (2.2). Any mode ˆ ω (cid:126)k with | (cid:126)k | > max { , δ } will from here on bereferred to as a “high mode”.To measure the relative energy in the low modes, we define the stochastic order parameter, Z vort ( t ) := | ˆ ω (1 , ( t ) | | ˆ ω (1 , ( t ) | + | ˆ ω (0 , ( t ) | , (2.3)4here ˆ ω (1 , and ˆ ω (0 , solve (2.1). Due to the condition, ˆ ω ( k ,k ) = ¯ˆ ω ( − k , − k ) , the relative energy in all ofthe low modes can be captured by Z vort ( t ). The value of Z vort ( t ), bounded between 0 and 1, correspondsto the proximity of the solution to an x -bar, y -bar or dipole state. If the dynamics drive Z vort ( t ) to increaseto 1, there is more energy in ˆ ω (1 , relative to ˆ ω (0 , , indicating the system is in an x -bar state. Converselyif Z vort ( t ) falls toward 0, the system would be observed to be in a y -bar state. If Z vort ( t ) instead staysnear 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 lowesteight Fourier modes, which for notational convenience we denote as ω := ˆ ω (1 , , ω := ˆ ω ( − , , ω := ˆ ω (0 , , ω := ˆ ω (0 , − ,ω := ˆ ω (1 , , ω := ˆ ω ( − , , ω := ˆ ω (1 , − , ω := ˆ ω ( − , − . (2.4)The variables ω , , , correspond to the low modes, while ω , , , represent the role of all the highmodes. Since the solution ω ( x, y ) of (2.1) is real valued, the following complex conjugacy relationshipmust still hold, ω = ¯ ω , ω = ¯ ω , ω = ¯ ω , ω = ¯ ω . (2.5)Thus the reduced model will be an eight dimensional approximation to the dynamics of (2.1). Toderive the model, we apply a center manifold reduction to (2.1) with σ (cid:126)k = 0 for all (cid:126)k to obtain an eight-dimensional deterministic ODE, which is the model studied in [1], and then add noise back to that systemto obtain the final eight-dimensional SDE model we study here.To carry out the center manifold reduction onto the lowest eight modes, assume for ˆ ω (cid:126)k with (cid:126)k / ∈{ ( ± , , (0 , ± , ( ± , ± } =: K , that there exists a smooth function H ( ω , . . . , ω ; (cid:126)k ) such that theeight-dimensional manifold defined by M = { ˆ ω : ˆ ω (cid:126)k = H ( ω , . . . , ω ; (cid:126)k ) , (cid:126)k / ∈ K } is invariant for the deterministic dynamics of (2.1) with σ (cid:126)k = 0 for all (cid:126)k . We refer to this as a centermanifold because it is defined in terms of the lowest eight modes, which have the weakest linear decayrates. Based on this assumption, one can then in principle compute the coefficients of the Taylor expansionof H ( · , (cid:126)k ) to any order for each (cid:126)k by taking the derivative of each of the low modes in two ways (via thefunction H and (2.1) with σ (cid:126)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 O ( ν ) neighborhood of zero due to the smallspectral gaps for the operator ν ∆. Additionally, while the existence of a finite dimensional (inertial) modelof the system (2.1) that describes the global dynamics cannot be expected [21], the model still providesmeaningful insight into the role δ plays in selecting the dominant quasi-stationary state for small initialconditions. For additional examples in which similar reductions of the Navier-Stokes equation to a finitedimensional model have been used to understand global dynamics see [7, 15].Adding independent (real) Brownian motions W , , , to each equation of the resulting ODE model5eads to our final SDE model˙ ω = − νδ ω + 1 δ (1 + δ ) [ ω ω − ¯ ω ω ] + 3 δ ν (4 + δ )(1 + δ ) ω ( | ω | + | ω | ) + √ νσ ˙ W ˙ ω = − νω + δ (1 + δ ) [¯ ω ω − ω ¯ ω ] + 3 δ ν (1 + 4 δ )(1 + δ ) ω ( | ω | + | ω | ) + √ νσ ˙ W ˙ ω = − ν δ δ ω − δ − δ ω ω − δ (3 + δ )2 ν (4 + δ )(1 + δ ) ω | ω | − δ νδ (1 + 4 δ )(1 + δ ) ω | ω | + √ νσ ˙ W ˙ ω = − ν δ δ ω + δ − δ ω ¯ ω − δ (3 + δ )2 ν (4 + δ )(1 + δ ) ω | ω | − δ νδ (1 + 4 δ )(1 + δ ) ω | ω | + √ νσ ˙ W . (2.6)Note that (2.6) with σ , , , = 0 corresponds to the ODE model derived in [1]. To compare thedynamics of this model to that of Z vort ( t ), defined in (2.3), we define the analogous order parameter forthe SDE model, Z red ( t ) := | ω ( t ) | | ω ( t ) | + | ω ( t ) | , (2.7)which again is used to determine towards which quasi-stationary state the system trends. Here, ω ( t )and ω ( t ) are solutions to the reduced system (2.6). The Monte Carlo simulation of the reduced modelfinds that the dominant quasi-stationary state depends on the aspect ratio of D δ in the same way asthe deterministic model, studied in detail in [1]. The main result there, which describes the selection ofquasi-stationary states in (2.6) with σ , , , = 0, can be described by the following theorem. Theorem 2.1. [1, Theorem 3.4] For δ ∈ (cid:16)(cid:113) , (cid:113) (cid:17) , under the dynamics of (2.6) with σ , , , = 0 , if δ > , then Z red ( t ) → , indicating evolution to an x -bar state. Conversely if δ < , then Z red ( t ) → ,indicating evolution to a y -bar state. For δ = 1 , there exists a one-dimensional center manifold of fixedpoints in the phase space that determines the asymptotic limit of Z red ( t ) . This center manifold is foliatedwith co-dimension one stable manifolds in which solutions converge to the corresponding fixed point. Exactlyone of these manifolds corresponds to each of the limits Z red ( t ) → and Z red ( t ) → . Thus, generic initialconditions are seen to evolve to the dipole state. Remark 2.2.
The order parameter considered in [1] was instead the ratio R ( t ) = | ω ( t ) | / | ω ( t ) | . The-orem 2.1 frames the result in terms of the order parameter Z red ( t ) . The choice to now consider Z red ( t ) isfor 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 δ , the set { Im( ω ) = Im( ω ) = Im( ω ) =Im( ω ) = 0 } is invariant under the dynamics of (2.6) with σ , , , = 0 . Since the real subsystem is invariantin the deterministic setting, we simulate the reduced model where the modes, ω , , , , as well as the Wienerprocesses, W , , , , are all real valued. This section provides simulations of the vorticity equation (2.1) and of the reduced model (2.6). Via MonteCarlo simulation, the average evolution of the order parameters Z vort ( t ) and Z red ( t ) will be plotted forseveral values of δ near 1. It will be seen that the reduced model captures the selection of the quasi-stationary states via the parameter δ . In particular, in both models, for a particular value of δ ≈
1, the6ystem’s selection of its dominant quasi-stationary state is consistent with the motivating results, given byTheorem 2.1. In particular, the system selects, as the dominant quasi-stationary state, a dipole for δ = 1,an x -bar for δ >
1, and a y -bar for δ < ω (cid:126)k with (cid:126)k ∈ K := { (cid:126)k = ( k , k ) ∈ Z : 0 ≤ | k | , | k | ≤
64 and ( k , k ) (cid:54) = (0 , } ; see [13]. A condition of exponential decay isimposed on the noise coefficients σ (cid:126)k seen in (1.4), | σ (cid:126)k | ≤ e − α | (cid:126)k | . (3.1)Similar to [6], simulations are conducted with (cid:80) { (cid:126)k ∈K} e − α | (cid:126)K | = 1. For our set K , this means α ≈ . N trials is plotted. Individual runs will be denoted by Z ivort ( t ) and Z ired ( t ), for i = 1 , . . . N , with corresponding averages given by¯ Z vort ( t ) = 1 N N (cid:88) i =1 Z ivort ( t ) , ¯ Z red ( t ) = 1 N N (cid:88) i =1 Z ired ( t ) . Similarly, we define the empirical variances to be V vort ( t ) = 1 N − N (cid:88) i =1 ( Z ivort ( t ) − ¯ Z vort ( t )) , V red ( t ) = 1 N − N (cid:88) i =1 ( Z ired ( t ) − ¯ Z red ( t )) . It will also be useful to plot the time averages of these Monte Carlo averages. To produce a meaningfulaverage we introduce a “burn-in time”, t burn , and ignore the initial period during which ¯ Z vort ( t ) and ¯ Z red ( t )have not yet stabilized. Define this time average for any function f ( t ) defined on t burn ≤ t ≤ T to be A ( f, t burn ) := 1 T − t burn (cid:90) Tt burn f ( t ) d t. Plotted in Figures 2-4 are ¯ Z vort ( t ), the time average A ( ¯ Z vort , t burn ), and the 95% confidence intervalsdefined via CI ± ( t ) = ¯ Z vort ( t ) ± . ∗ (cid:114) V vort ( t ) N .
Also included are average contour plots for the vorticity. We use N = 200 and for each trial use zeroinitial conditions and ν = 0 . δ = 1, Figure 2a shows ¯ Z vort ( t ) remains near 1/2 for the durationof the simulation. We use a burn-in time of t burn = 0 when computing the time average since on thesymmetric domain it is clear there is no transient initial period. In Figure 2b, the average contour plot foreach individual trial are themselves averaged over the N = 200 trials, reflecting a dipole.7 a) ¯ Z vort ( t ) with 95% confidence interval. (b) Average contour plot of vorticity. Figure 2: Vorticity aligns on average as a dipole for δ = 1.The simulations exhibited in Figures 3a and 3b show that, for δ = 1 .
1, the order parameter increasesinitially and the average contour plot looks like that of an x -bar state. In Figure 3a, t burn = 100 is usedwhen computing the time average. (a) ¯ Z vort ( t ) with 95% confidence interval. (b) Contour plot of vorticity. Figure 3: Vorticity aligns on average as an x -bar for δ = 1 . δ < δ = 0 .
9, the orderparameter decreases over an initial period of time and the average contour plot looks like that of a y -barstate. Here we again set t burn = 100. 8 a) ¯ Z vort ( t ) with 95% confidence interval. (b) Average contour plot of vorticity. Figure 4: Vorticity aligns on average as a y -bar for δ = 0 . Z vort ( t ) for δ = 1 . δ = 1 . δ = 0 .
90 averaged over N = 1000trials. This is to show that as the number of trials increase, the variance is decreasing without changing themean behavior. The variances all remain generally between 0.06-0.08. While the variance does decrease,the limiting value of ¯ Z vort ( t ) remains relatively unchanged compared to what is seen when averaging over N = 200 trials. (a) δ = 1 .
10 (b) δ = 1 . δ = 0 . Figure 5: Plot of ¯ Z vort ( t ) and of 95% confidence level error bars with N = 1000 trials and ν = 0 . §
5, we also include in Figure 6 a simulation with ν = 0 . δ = 1 .
04, the value of δ for which transitions among the quasi-stationary states were observed in [6]. 9 a) An individual trajectory transitions among quasi-stationary states. (b) On average, the system is close to an x -bar state. Figure 6: A single trajectory and the Monte Carlo average ¯ Z vort ( t ) for δ = 1 . δ . However Figure 6bshows that E [ Z vort ( t )] picks the dominant state. We now compute the time average of a randomly selectedindividual trial, given by A ( Z vort , t burn ), to confirm that it tracks the Monte Carlo average, ¯ Z vort ( t ). Figure7 shows two things. First, for the given values of δ , a sample path may experience many transitions amongthe quasi-stationary states. Second, the time average of the sample path does eventually track the MonteCarlo average. (a) δ = 1 .
10 (b) δ = 1 . δ = 0 . Figure 7: Comparing individual time average of a sample path with Monte Carlo average
We now turn our attention to the reduced model (2.6). We confirm numerically that the reduced modelcaptures the qualitative dynamics of the full vorticity equation with regard to the dominant quasi-stationarystate.As we mentioned in Remark 2.3 we will be working with the real system in which ω , , , , as well asthe Wiener processes, W , , , , are all real valued. This leads to the following system which serves as the10cting reduced model in the upcoming simulations.˙ ω = − νδ ω + 1 δ (1 + δ ) [ ω ω − ω ω ] + 3 δ ν (4 + δ )(1 + δ ) ω ( ω + ω ) + √ νσ ˙ W ˙ ω = − νω + δ (1 + δ ) [ ω ω − ω ω ] + 3 δ ν (1 + 4 δ )(1 + δ ) ω ( ω + ω ) + √ νσ ˙ W ˙ ω = − ν δ δ ω − δ − δ ω ω − δ (3 + δ )2 ν (4 + δ )(1 + δ ) ω ω − δ νδ (1 + 4 δ )(1 + δ ) ω ω + √ νσ ˙ W ˙ ω = − ν δ δ ω + δ − δ ω ω − δ (3 + δ )2 ν (4 + δ )(1 + δ ) ω ω − δ νδ (1 + 4 δ )(1 + δ ) ω ω + √ νσ ˙ W . (3.2)To be consistent with the spatial decay of the noise in the simulations of the stochastically forcedvorticity equation (2.1), given by (3.1), we choose σ , = e − α and σ , = e − α . First we aim to establish that the reduced model (3.2) can serve as a good approximation to thevorticity equation with noise, (2.1), for δ ≈
1. Second, it will be established that the selection of the baror dipole state that dominates is consistent with the results of [1] for the deterministic equation: x -bar for δ > y -bar for δ <
1, and dipole for δ = 1.Figure 8 shows numerical evidence supporting that the dynamics of the order parameter, governedby the reduced system (3.2), follows the same trend as when the full vorticity equation is simulated.Figure 8: Simulation of ¯ Z red ( t ) with noise for ν = 0 . N = 200 trials) show that the trendtoward the appropriate quasi-stationary state is captured by the reduced model. Starting with zero initialconditions, when the noise is added, the simulations show that for δ >
1, the order parameter increasestoward 1, indicating evolution to an x -bar state. Conversely, for δ <
1, the order parameter decreases11oward a value corresponding to a y -bar state. Finally, when δ = 1, ¯ Z red ( t ) remains near 1/2 indicatingthe system is in a dipole state. Figures 9a-9e serve to compare the evolution of ¯ Z vort ( t ) and ¯ Z red ( t )takenover N = 200 trials, for values of δ close to 1. The bars denote the error for the 95% confidence intervalsfor ¯ Z vort (bold) and ¯ Z red (thin). (a) δ = 0 .
90 (b) δ = 0 .
95 (c) δ = 1 . δ = 1 .
05 (e) δ = 1 . Figure 9: Comparing ¯ Z vort ( t ) and ¯ Z red ( t ) averaged over N=200 trials with ν = 0 . Z vort ( t ) and ¯ Z red ( t ) both trend in the same direction, with similar variances (typicallybetween 0.06-0.08 for 0 ≤ t ≤ δ. Motivated by the numerics from §
3, this section investigates the expected behavior of Z red ( t ) as δ → δ = 1 and ν = 0 case. Using the backward Kolmogorovequation associated to (3.2), the goal is to derive a system of PDE that will provide insight on how theexpected value of Z red ( t ), to leading order, depends on values of δ close to 1. To do this we pose theproblem as a perturbation of the spatial domain, setting δ = 1 + (cid:15) (cid:15) . Here, 0 < (cid:15) (cid:28) (cid:15) = ± E [ Z red ( t )] as (cid:15) → p := Re( ω ) , ˜ q := Re( ω ) , ˜ r := Re( ω ) , ˜ s := Re( ω ) . Now (3.2) can be expressed as˙˜ p = − νδ ˜ p + 1 δ (1 + δ ) ˜ q (˜ s − ˜ r ) + 3 δ ν (4 + δ )(1 + δ ) ˜ p ( | ˜ r | + | ˜ s | ) + √ νσ ˙ W ˙˜ q = − ν ˜ q + δ (1 + δ ) ˜ p (˜ r − ˜ s ) + 3 δ ν (1 + 4 δ )(1 + δ ) ˜ q ( | ˜ r | + | ˜ s | ) + √ νσ ˙ W (4.1)˙˜ r = − ν δ δ ˜ r − δ − δ ˜ p ˜ q − δ (3 + δ )2 ν (4 + δ )(1 + δ ) ˜ r | ˜ p | − δ νδ (1 + 4 δ )(1 + δ ) ˜ r | ˜ y | + √ νσ ˙ W ˙˜ s = − ν δ δ ˜ s + δ − δ ˜ p ˜ q − δ (3 + δ )2 ν (4 + δ )(1 + δ ) ˜ s | ˜ p | − δ νδ (1 + 4 δ )(1 + δ ) ˜ s | ˜ q | + √ νσ ˙ W . Before inserting the Taylor expansions in (cid:15) for the coefficients with δ = 1 + (cid:15) (cid:15) , we first scale (4.1)appropriately to obtain a clear slow-fast system. As in [1], the low modes represented by ˜ p and ˜ q correspondto the slow variables while the high modes, ˜ r and ˜ s , represent the fast variables. Below, we give a moregeneral version of the scaled equations for just the ˜ p (analogous to ˜ q ) and ˜ r (analogous to ˜ s ) equations.We use the following space-time and parameter scalings: ν = (cid:15) µ ν , ˜ p = (cid:15) ξ p , ˜ q = (cid:15) ξ q , ˜ r = (cid:15) η r , ˜ s = (cid:15) η s ,and τ = (cid:15) γ t . To simplify the scaled equations, we will relate µ , ξ , η and γ to put the resulting system ina more desirable form. We neglect the (cid:15) dependence of p , q , r and s for readability. Below, the “prime”notation denotes differentiation with respect to the scaled time variable, τ . p (cid:48) = (cid:15) µ − γ ( − ν δ p ) + (cid:15) η − γ δ (1 + δ ) q ( s − r ) + (cid:15) η − µ − γ δ ν (4 + δ )(1 + δ ) p ( r + s ) + (cid:15) µ − γ − ξ √ ν σ W (cid:48) ( τ ) r (cid:48) = (cid:15) µ − γ ( − ν δ δ r ) − (cid:15) ξ − γ − η δ − δ ( pq ) − (cid:15) ξ − µ − γ ν (cid:18) δ (3 + δ )2(4 + δ )(1 + δ ) rp + 1 + 3 δ δ )(1 + δ ) rq (cid:19) + (cid:15) µ − γ − η √ ν σ W (cid:48) ( τ )Now set 2 η = µ + γ , 2 ξ = µ − γ ⇒ γ = µ − ξ , with 0 < γ < ξ < µ < η < µ. Then the fully scaledsystem (still neglecting Taylor expansions of δ in (cid:15) for now) becomes p (cid:48) = (cid:15) ξ ( − ν δ p ) + (cid:15) ξ δ (1 + δ ) q ( s − r ) + 3 δ ν (4 + δ )(1 + δ ) p ( | s | + | s | ) + √ ν σ W (cid:48) ( τ ) q (cid:48) = (cid:15) ξ ( − ν q ) + (cid:15) ξ δ (1 + δ ) p ( r − s ) + 3 δ ν (1 + 4 δ )(1 + δ ) q ( | r | + | s | ) + √ ν σ W (cid:48) ( τ ) (4.2) r (cid:48) = (cid:15) ξ (cid:18) − ν δ δ r (cid:19) − (cid:15) ξ − η δ − δ pq − (cid:15) ξ − η ) (cid:18) δ (3 + δ )2 ν (4 + δ )(1 + δ ) r | p | + 1 + 3 δ ν δ (1 + 4 δ )(1 + δ ) r | q | (cid:19) + (cid:15) ξ − η √ ν σ W (cid:48) ( τ ) s (cid:48) = (cid:15) ξ (cid:18) − ν δ δ s (cid:19) + (cid:15) ξ − η δ − δ pq − (cid:15) ξ − η ) (cid:18) δ (3 + δ )2 ν (4 + δ )(1 + δ ) s | p | − δ ν δ (1 + 4 δ )(1 + δ ) s | q | (cid:19) + (cid:15) ξ − η √ ν σ W (cid:48) ( τ ) . For the scaled SDE (4.2), let b (cid:15) = ( b (cid:15)p , b (cid:15)q , b (cid:15)r , b (cid:15)s ) denote the drift vector and Σ (cid:15) ( p, q, r, s ; σ , σ , σ , σ )denote the diffusion matrix so that (4.2) can be expressed, for X (cid:15) = ( p, q, r, s ) and d W dτ = (cid:16) dW dτ , dW dτ , dW dτ , dW dτ (cid:17) ,13s d X (cid:15) dτ = b (cid:15) + Σ (cid:15) d W dτ . (4.3)Now replacing the δ coefficients appearing in (4.2) with their Taylor expansions for δ = 1 + (cid:15) (cid:15) upto O ( (cid:15) ), the drift vector is given by (still supressing (cid:15) dependence of p, q, r, s ), b (cid:15)p = 1 ν (cid:18)
340 + (cid:15)(cid:15) (cid:15) − (cid:15) (cid:15) (cid:19) p ( r + s )+ (cid:15) ξ (cid:18) − (cid:15) (cid:15) + 716 (cid:15) − (cid:15) (cid:15) (cid:19) q ( s − r ) − (cid:15) ξ ν (1 − (cid:15) (cid:15) + (cid:15) − (cid:15) (cid:15) ) p (4.4) b (cid:15)q = 1 ν (cid:18) − (cid:15) (cid:15) + 1174000 (cid:15) − (cid:15) (cid:15) (cid:19) q ( r + s )+ (cid:15) ξ (cid:18)
12 + (cid:15) − (cid:15) (cid:19) p ( r − s ) − (cid:15) ξ ν qb (cid:15)r = (cid:15) ξ − η ) ν (cid:20) − r ( p + q ) − (cid:15) (cid:15) r (51 p − q ) − (cid:15) r ( p + q ) − (cid:15) (cid:15) r (cid:0) p − q (cid:1)(cid:21) − (cid:15) ξ ν (2 − (cid:15) (cid:15) + (cid:15) − (cid:15) (cid:15) ) r − (cid:15) ξ − η (cid:18) (cid:15) (cid:15) − (cid:15) + 38 (cid:15) (cid:15) (cid:19) pqb (cid:15)s = (cid:15) ξ − η ) ν (cid:20) − s ( p + q ) − (cid:15) (cid:15) s (51 p − q ) − (cid:15) s ( p + q ) − (cid:15) (cid:15) s (cid:0) p − q (cid:1)(cid:21) − (cid:15) ξ ν (2 − (cid:15) (cid:15) + (cid:15) − (cid:15) (cid:15) ) s + (cid:15) ξ − η (cid:18) (cid:15) (cid:15) − (cid:15) + 38 (cid:15) (cid:15) (cid:19) pq and the diffusion matrix byΣ (cid:15) ( p, q, r, s ; σ , σ , σ , σ ) = √ ν σ √ ν σ (cid:15) ξ − η √ ν σ
00 0 0 (cid:15) ξ − η √ ν σ . (4.5)With H( u (cid:15) ) denoting the Hessian matrix of u (cid:15) , we now write the backward Kolmogorov equation for(4.3), which is defined as ∂u (cid:15) ∂τ = b (cid:15) · ∇ u (cid:15) + 12 Tr[(Σ (cid:15) ) H( u (cid:15) )] , in R × [0 , T ] u (cid:15) ( p, q, r, s,
0) = φ ( p, q ) , on R × { } . (4.6)The backward Kolmogorov equation has the useful property that the evolution of u (cid:15) ( X (cid:15) , τ ) gives u (cid:15) ( p, q, r, s, τ ) = E [ φ ( p τ , q τ ) | p τ (0) = p, q τ (0) = q, r τ (0) = r, s τ (0) = s ] . Thus one ultimately is interested in initializing (4.6) with φ ( p, q ) = Z red = p p + q , but for now weproceed with a general initial condition, φ . We seek a solution to (4.6) that takes the form u (cid:15) ( p, q, r, s, τ ) = u ( p, q, r, s, τ ) + (cid:15)u ( p, q, r, s, τ ) + (cid:15) u ( p, q, r, s, τ ) + . . . (4.7)and wish to find the limiting dynamics, u (cid:15) as (cid:15) →
0. As the goal is to identify the leading order expansionfor u (cid:15) we determine a system of PDEs for u , u , and u . We present now the calculations that lead tothe characterization of u , u , and u , see (4.16). 14efine L (cid:15) to be the operator acting on the right hand side of (4.6), so that ∂u (cid:15) ∂τ = L (cid:15) u (cid:15) . Decomposing L (cid:15) by powers of (cid:15) using the expressions for b (cid:15) and Σ (cid:15) given in (4.4) and (4.5), we write L (cid:15) u (cid:15) = (cid:15) ξ − η ) L u (cid:15) + (cid:15) ξ − η )+1 L u (cid:15) + (cid:15) ξ − η )+2 L u (cid:15) + (cid:15) ξ − η )+3 L u (cid:15) + L u (cid:15) + (cid:15) L u (cid:15) + (cid:15) L + (cid:15) L u (cid:15) + (cid:15) ξ L u (cid:15) + (cid:15) ξ +1 L u (cid:15) + (cid:15) ξ +2 L u (cid:15) + (cid:15) ξ +3 L u (cid:15) + (cid:15) ξ L u (cid:15) + (cid:15) ξ +1 L u (cid:15) + (cid:15) ξ +2 L u (cid:15) + (cid:15) ξ +3 L u (cid:15) + (cid:15) ξ − η +1 L u (cid:15) + (cid:15) ξ − η +2 L u (cid:15) + (cid:15) ξ − η +3 L u (cid:15) . A select few of the operators, L i , i = 1 , . . . ,
19, that are most important in computing the leadingorder equations is provided in (4.8). A complete list of the expressions of the 19 operators can be found inthe appendix. L u = − ν ( p + q ) (cid:18) r ∂u∂r + s ∂u∂s (cid:19) + 2 ν (cid:18) σ ∂ u∂r + σ ∂ u∂s (cid:19) L u = − (cid:15) ν (51 p − q ) (cid:18) r ∂u∂r + s ∂u∂s (cid:19) L u = 340 ν ( r + s ) (cid:18) p ∂u∂p + q ∂u∂q (cid:19) + 2 ν (cid:18) σ ∂ ∂p + σ ∂ ∂q (cid:19) (4.8) L u = (cid:15) ν ( r + s ) (cid:18) p ∂u∂p − q ∂u∂q (cid:19) L u = −
12 ( r − s ) (cid:18) q ∂u∂p − p ∂u∂q (cid:19) Now we choose explicit values for ξ and η to obtain a simpler, but still representative system: η = 2, ξ = 1, hence µ = 3 ( ν = (cid:15) µ ν = (cid:15) ν ). Note that this choice of parameters is not unique. Our goal is tomake a choice that makes the computations that follow tractable, preserves a clear slow-fast system in thescaled equations (4.2) and results in a system that exhibits qualitatively the same behavior as the originalsystem in terms of the selection mechanism to x -bar, y -bar and dipole states depending on the values of δ . With this choice, the generator L (cid:15) of the PDE becomes L (cid:15) u (cid:15) = (cid:15) − L u (cid:15) + (cid:15) − L u (cid:15) + ( L + L + L ) u (cid:15) + (cid:15) ( L + L + L + L ) u (cid:15) (4.9)+ (cid:15) ( L + L + L + L ) u (cid:15) + (cid:15) ( L + L + L ) u (cid:15) + (cid:15) ( L + L ) u (cid:15) + (cid:15) L u (cid:15) . The ansatz given in (4.7) can now be inserted into the backward Kolmogorov equation (4.6) usingthe expression of L (cid:15) given above in (4.9). Matching coefficients on both sides of the equation yields thefollowing leading order equations, O ( (cid:15) − ) : −L u = 0 (4.10a) O ( (cid:15) − ) : −L u = L u (4.10b) O (1) : −L u = − ∂u ∂τ + L u + ( L + L + L ) u (4.10c) O ( (cid:15) ) : −L u = − ∂u ∂τ + L u + ( L + L + L ) u + ( L + L + L + L ) u . (4.10d)15quation (4.10a) implies u lies in the kernel of L , which elliptic PDE theory tells us contains onlyfunctions constant in r and s . Since L is also a differential operator in r and s only, (4.10b) implies that u is constant in r and s as well. One can also see that u and u are in the kernel of each of L , , , (seeappendix). Hence the leading order system given by (4.10a)-(4.10d) can be reduced to O ( (cid:15) − ) : −L u = 0 ⇒ u = u ( p, q, τ ) (4.11a) O ( (cid:15) − ) : −L u = L u ⇒ u = u ( p, q, τ ) (4.11b) O (1) : −L u = − ∂u ∂τ + L u (4.11c) O ( (cid:15) ) : −L u = − ∂u ∂τ + L u + L u + ( L + L ) u , (4.11d)where L , L , L , L , and L are presented in (4.8). Let ρ ∞ ( r, s ; p, q ) be the stationary density thatsatisfies the adjoint problem L ∗ ρ ∞ ( r, s ; p, q ) = 0 . Once ρ ∞ is known we can integrate against the invariant measure to obtain the solvability conditionsfor equations (4.11c) and (4.11d) ∂u ∂τ = (cid:90) R L u ρ ∞ ( r, s ; p, q )d r d s (4.12) ∂u ∂τ = (cid:90) R ( L u + L u + ( L + L ) u ) ρ ∞ ( r, s ; p, q )d r d s. Before we consider the integrals in (4.12), ρ ∞ must be identified. The operator, L = (cid:18) − ν r ( p + q ) , − ν s ( p + q ) (cid:19) · (cid:18) ∂∂r , ∂∂s (cid:19) + 12 (cid:32) ν σ
00 4 ν σ (cid:33) (cid:32) ∂ ∂r ∂ ∂s (cid:33) , corresponds to the backward Kolmogorov equation for the following system, parameterized by the fixed(slow) variables p and q . ˙ˆ r = − ν ( p + q )ˆ r + σ √ ν ˙ W ˙ˆ s = − ν ( p + q )ˆ s + σ √ ν ˙ W . These processes are independent Ornstein-Uhlenbeck processes and are therefore Gaussian. The equi-librium (stationary) density which corresponds to ρ ∞ ( r, s ; p, q ) is that of the bivariate Gaussian distributionwith r ∼ N (cid:18) , ν σ p + q (cid:19) , s ∼ N (cid:18) , ν σ p + q (cid:19) . Therefore the invariant joint density is ρ ∞ ( r, s, p, q ) = p + q πν σ σ e − p q ν ( r σ + s σ ) . To aid in the computations of the integrals given in (4.12), the following integral evaluations will beuseful and can be simply obtained through the mean and variance of the stationary distribution. (cid:90) R ( r + s ) ρ ∞ ( r, s ; p, q )d r d s = 5 ν p + q ( σ + σ ) (4.13a) (cid:90) R ( r − s ) ρ ∞ ( r, s ; p, q )d r d s = 0 . (4.13b)16ext consider the solvability conditions (4.12) one at a time. From the u equation and the integral(4.13a), ∂u ∂τ = (cid:90) R L u ρ ∞ d r d s = 2 ν (cid:18) σ ∂ u ∂p + σ ∂ u ∂q (cid:19) + 340 ν (cid:18) p ∂u ∂p + q ∂u ∂q (cid:19) (cid:90) R ( r + s ) ρ ∞ ( r, s ; p, q )d r d s = 2 ν (cid:18) σ ∂ u ∂p + σ ∂ u ∂q (cid:19) + 3 ν σ + σ ) (cid:18) pp + q ∂u ∂p + qp + q ∂u ∂q (cid:19) . From this we obtain the effective equations for p and q for small (cid:15) after the fast variables r and s areaveraged out. The slow motion can hence be approximated, for 0 < (cid:15) (cid:28)
1, by ¯ p and ¯ q governed by,¯ p (cid:48) = 3 ν σ + σ ) ¯ p ¯ p + ¯ q + σ √ ν W (cid:48) ¯ q (cid:48) = 3 ν σ + σ ) ¯ q ¯ p + ¯ q + σ √ ν W (cid:48) . (4.14)Since (cid:15) dependence does not appear in the first order equations, we will need to determine u tosee its effects. Consider the solvability condition for u in (4.12). Computing this integral requires us toevaluate the following integrals. I = (cid:90) R L u ρ ∞ d r d sI (cid:48) = (cid:90) R L u ρ ∞ d r d sI = (cid:90) R L u ρ ∞ d r d sI = (cid:90) R L u ρ ∞ d r d s. In evaluating these, we see I = (cid:90) R (cid:15) ν (cid:18) p ∂u ∂p − q ∂u ∂q (cid:19) ( r + s ) ρ ∞ ( r, s ; p, q )d r d s == 5 (cid:15) ν ( σ + σ ) (cid:18) pp + q ∂u ∂p − qp + q ∂u ∂q (cid:19) I (cid:48) = (cid:90) R − (cid:18) q ∂u ∂p − p ∂u ∂q (cid:19) ( r − s ) ρ ∞ ( r, s ; p, q )d r d s = 0 I = 2 ν (cid:18) σ ∂ u ∂p + σ ∂ u ∂q (cid:19) + 3 ν σ + σ ) (cid:18) pp + q ∂u ∂p + qp + q ∂u ∂q (cid:19) I = − (cid:15) ν (51 p − q ) (cid:90) R (cid:18) r ∂u ∂r + s ∂u ∂s (cid:19) ρ ∞ ( r, s ; p, q )d r d s. (4.15)Since u depends on the fast variables, this final integral cannot yet be computed. It will eventually17e handled numerically. Thus, formally, we have u (cid:15) = u + (cid:15)u + (cid:15) u + . . . satisfying, ∂u ∂τ = 3 ν σ + σ ) (cid:18) pp + q ∂u ∂p + qp + q ∂u ∂q (cid:19) + 2 ν (cid:18) σ ∂ u ∂p + σ ∂ u ∂q (cid:19) (4.16) ∂u ∂τ = 3 ν σ + σ ) (cid:18) pp + q ∂u ∂p + qp + q ∂u ∂q (cid:19) + 2 ν (cid:18) σ ∂ u ∂p + σ ∂ u ∂q (cid:19) + 5 (cid:15) ν ( σ + σ ) (cid:18) pp + q ∂u ∂p − qp + q ∂u ∂q (cid:19) − (cid:15) ν (51 p − q ) (cid:90) R (cid:18) r ∂u ∂r + s ∂u ∂s (cid:19) d ρ ∞ −L u = ∂u ∂τ − L u We shall consider the system (4.16) together with the initial conditions u ( p, q,
0) = φ ( p, q ) , u ( p, q,
0) = 0 , and u ( p, q, r, s,
0) = 0 . The PDE for u immediately stands out as the backward Kolmogorov equation corresponding to thesystem given in (4.14). Despite its simple looking form, the regularity at the origin of the coefficients onthe first derivative terms turn out to be a borderline case with regards to well posedness, see for exampleChapter III, Section 1 of [11]. Nevertheless, we proceed formally and solve for u , u , and u numericallyafter providing the initial condition φ ( p, q,
0) = p p + q so that u (cid:15) ( p, q, r, s, τ ) = E [ Z red ( τ ) | p = p, q = q, r = r, s = s ]. The simulations of the system (4.16)provide an approximation to the deterministic evolution of E [ Z red ( τ )] after averaging out the fast motion.The simulations of the system (4.16) provided in this section were conducted via finite differences on thedomain ( p, q, r, s, τ ) ∈ [ − , × [0 , T ] with Neumann boundary conditions.When (cid:15) = 0, which implies δ = 1, the system is unperturbed and u ( p, q, t ) = E [ Z red ( τ ) | p = p, q = q ]. In this case, the results of the simulation show that the expected value of the order parameter Z red ( τ )converges to 1/2 for any initial values p = p and q = q of (2.6), independent of r and s . This indicatesthat the unperturbed system evolves to a dipole state, even if the initial state is close to an x - or y -barstate. Figure 10 illustrates the evolution to a dipole for (cid:15) = 0 for several initial conditions ( p, q ) chosenwithin the domain. Figure 10: For (cid:15) = 0, E [ Z red ( τ )] → / E [ Z red ( τ )], given by u (cid:15) ( τ ) tothe average path of the order parameter Z red ( t ), i.e. ¯ Z red ( t ), as computed via Monte Carlo simulation in §
3. Since the two models evolve on different timescales, we rescale τ so that our averaged PDE modelis evolving on the original timescale. As such, suppressing the spatial ( p, q, r, s ) dependency, let ˆ u (cid:15) ( t ) := u ( t/(cid:15) ) + (cid:15)u ( t/(cid:15) ) = u ( τ ) + (cid:15)u ( τ ) denote the O ( (cid:15) ) approximation to u (cid:15) in the original timescale, and let¯ Z red ( t ) denote the Monte Carlo average path of the order parameter under the dynamics of (2.6) obtainedvia the Monte Carlo simulations described in § u (cid:15) ( p, q, t ) for ( p, q ) = (0 . , . O ( (cid:15) ) approximation to u (cid:15) ( τ ) = E [ Z red ( τ )] evolves toward 0 or 1 depending on thesign of ˆ (cid:15) := (cid:15) (cid:15). (a) Approximation evolves to an y -bar state for ˆ (cid:15) = 0 .
1. (b) Approximation evolves to an x -bar state for ˆ (cid:15) = − . Figure 11: Approximation follows the evolution of the Z red ( t ) for small (cid:15) .Using these simulations with ν scaled, we now explore how the intervals on which ˆ u (cid:15) serves as a goodapproximation to ¯ Z red ( t ) depend on the perturbation parameter. Figure 12 shows the relative error (RE)given by RE = | ˆ u (cid:15) ( t ) − ¯ Z red ( t ) | ¯ Z red ( t ) . | ˆ u (cid:15) ( t ) − ¯ Z red ( t ) | ¯ Z red ( t ) .On some initial interval of time, the PDE approximation, ˆ u (cid:15) indeed serves as a close approximation to¯ Z red ( t ). Furthermore, as expected, the smaller the perturbation parameter (cid:15) , the longer the approximationis valid. In this paper, we developed a finite dimensional SDE model that can be used to elucidate the dynamicsof the 2D Navier-Stokes vorticity equation with noise. Monte Carlo simulation of the reduced modelshowed 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, theexistence and attracting nature of these quasi-stationary states play an important role in the evolutionof the stochastic Navier-Stokes vorticity equation. Specifically, the aspect ratio of the periodic domain, D δ = [0 , πδ ] × [0 , π ], determines whether generic solutions evolve toward an x -bar state ( δ > y -barstate ( δ < δ = 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 perturbationparameter δ ≈ δ = 1 , in relation to how the viscocity parameter ν vanishes.The numerical studies in § δ , see Figure 6b. However, one can see from the sample path plottedin Figure 6a, individual sample paths do exhibit transitions between x -bar and y -bar states, as it has alsobeen observed in [6].In regards to future directions, there are a number of interesting questions that one can ask and hopeto answer. To begin with, the perturbation analysis of § § δ < δ >
1. One would like to make this mathematically rigorous. Furthermore, one could potentially use the20educed model of § √ ν , is common in theliterature, see for example [16] and the references therein. One of the important questions in the literatureis the investigation of the convergence of the corresponding invariant measure as ν → ν → ν → The complete list of operators in the Kolmogorov equation (4.6) is given by L u = − ν ( p + q ) (cid:18) r ∂u∂r + s ∂u∂s (cid:19) + 2 ν (cid:18) σ ∂ u∂r + σ ∂ u∂s (cid:19) L u = − (cid:15) ν (51 p − q ) (cid:18) r ∂u∂r + s ∂u∂s (cid:19) L u = − ν ( p + q ) (cid:18) r ∂u∂r + s ∂u∂s (cid:19) L u = − (cid:15) ν (379 p − q ) (cid:18) r ∂u∂r + s ∂u∂s (cid:19) L u = 340 ν ( r + s ) (cid:18) p ∂u∂p + q ∂u∂q (cid:19) + 2 ν (cid:18) σ ∂ ∂p + σ ∂ ∂q (cid:19) L u = (cid:15) ν ( r + s ) (cid:18) p ∂u∂p − q ∂u∂q (cid:19) L u = 1174000 ν ( r + s ) (cid:18) p ∂u∂p + q ∂u∂q (cid:19) L u = − (cid:15) ν ( r + s ) (cid:18) p ∂u∂p + 9320000 q ∂u∂q (cid:19) L u = −
12 ( r − s ) (cid:18) q ∂u∂p − p ∂u∂q (cid:19) L u = (cid:15) r − s ) (cid:18) q ∂u∂p + p ∂u∂q (cid:19) (6.1) L u = 18 ( r − s ) (cid:18) q ∂u∂p − p ∂u∂q (cid:19) L u = (cid:15) q ( r − s ) ∂u∂p u = − ν (cid:18) p ∂u∂p + q ∂u∂q + 2 (cid:18) r ∂u∂r + s ∂u∂s (cid:19)(cid:19) L u = ν (cid:15) (cid:18) p ∂u∂p + r ∂u∂r + s ∂u∂s (cid:19) L u = − ν (cid:18) p ∂u∂p + r ∂u∂r + s ∂u∂s (cid:19) L u = ν (cid:15) (cid:18) p ∂u∂p + r ∂u∂r + s ∂u∂s (cid:19) L u = − (cid:15) pq (cid:18) ∂u∂r − ∂u∂s (cid:19) L u = 12 pq (cid:18) ∂u∂r − ∂u∂s (cid:19) L u = − (cid:15) pq (cid:18) ∂u∂r − ∂u∂s (cid:19) References [1] M. Beck, E. Cooper, and K. Spiliopoulos. Selection of quasi-stationary states in the navier-stokesequation on the torus.
Nonlinearity , 32(1):209–237, 2019.[2] M. Beck and C. E. Wayne. Metastability and rapid convergence to quasi-stationary bar states for thetwo-dimensional Navier–Stokes equations.
Proc. Roy. Soc. Edinburgh Sect. A , 143(5):905–927, 2013.[3] J. Bedrossian and M. Coti Zelati. Enhanced dissipation, hypoellipticity, and anomalous small noiseinviscid limits in shear flows.
Archive for Rational Mechanics and Analysis , 224(3):1161–1204, 2017.[4] J. Bedrossian, M. Coti Zelati, and N. Glatt-Holtz. Invariant measures for passive scalars in the smallnoise inviscid limit.
Communications in Mathematical Physics , 348(1):101–127, 2016.[5] A. Bensoussan and R. Temam. ´Equations stochastiques du type navier-stokes.
J. Functional Analysis ,13:195–222, 1973.[6] F. Bouchet and E. Simonnet. Random changes of flow topology in two-dimensional and geophysicalturbulence.
Phys. Rev. Lett. , 102(094504), 2009.[7] W. E and J. C. Mattingly. Ergodicity for the Navier-Stokes equation with degenerate random forcing:finite-dimensional approximation.
Comm. Pure Appl. Math. , 54(11):1386–1402, 2001.[8] E. Grenier, T.T. Nguyen, F. Rousset, and A. Soffer. Linear inviscid damping and enhanced viscousdissipation of shear flows by using the conjugate operator method. preprint , 2017.[9] S. Ibrahim, Y. Maekawa, and N. Masmoudi. On pseudospectral bound for non-selfadjoint operatorsand its application to stability of kolmogorov flows.
Preprint , https://arxiv.org/pdf/1710.05132.pdf,2017.[10] S. Kuksin and A. Shirikyan.
Mathematics of Two-Dimensional Turbulence . Cambridge Tracts inMathematics. Cambridge University Press, 2012.2211] O.A. Ladyzenskaja, V.A. Solonnikov, and N.N. Ural’ceva.
Linear and Quasi-linear Equations ofParabolic Type . American Mathematical Society, 1968.[12] Z. Lin and M. Xu. Metastability of kolmogorov flows and inviscid damping of shear flows.
Archivefor Rational Mechanics and Analysis , 231:1811–1852, 2019.[13] G.L. Lord, C.E. Powell, and T. Shardlow.
An Introduction to Computational Stochastic PDEs . Cam-bridge University Press, 2014.[14] J. C. Mattingly. The dissipative scale of the stochastics Navier-Stokes equation: regularization andanalyticity.
J. Statist. Phys. , 108(5-6):1157–1179, 2002. Dedicated to David Ruelle and Yasha Sinaion the occasion of their 65th birthdays.[15] J. C. Mattingly and E. Pardoux. Invariant measure selection by noise: an example.
Discrete andContinuous Dynamical Systems , 34(10):4223–4257, 2014.[16] Glatt-Holtz N., V. ˇSver´ak, and V. Vicol. On inviscid limits for the stochastic navier–stokes equationsand related models.
Archive for Rational Mechanics and Analysis , 217(2):619–649, 2015.[17] E.A. Novikov. Functionals and the random-force method in turbulence theory.
Soviet Physics JETP ,20:1290–1294, 1965.[18] G. Pavliotis and A. Stuart.
Multiscale Methods: Averaging and Homogenization . Springer, 2008.[19] D. Wei, Z. Zhang, and W. Zhao. Linear inviscid damping and enhanced dissipation for the kolmogorovflow. preprint , 2017.[20] Z. Yin, D.C. Montgomery, and H.J.H. Clercx. Alternative statistical-mechanical descriptions of de-caying two-dimensional turbulence in terms of “patches” and “points”.
Phys. Fluids , 15:1937–1953,2003.[21] S. Zelik. Inertial manifolds and finite-dimensional reduction for dissipative pdes.