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

TL;DR

This paper investigates how the aspect ratio parameter elta in the 2D Navier-Stokes equations on a torus influences the selection between bar and dipole quasi-stationary states, using a Fourier space center manifold model.

Contribution

It provides a mathematical analysis showing elta's role in state selection, supported by energy estimates, invariant manifolds, asymptotic expansions, and numerical evidence.

Findings

elta determines the likelihood of observing bar or dipole states.

High Fourier modes decay rapidly initially, then low modes evolve slowly.

The model confirms the elta-dependent selection mechanism.

Abstract

The two dimensional incompressible Navier-Stokes equation on $D_\delta := [0, 2\pi\delta] \times [0, 2\pi]$ with $\delta \approx 1$, periodic boundary conditions, and viscosity $0 < \nu \ll 1$ is considered. Bars and dipoles, two explicitly given quasi-stationary states of the system, evolve on the time scale $\mathcal{O}(e^{-\nu t})$ and have been shown to play a key role in its long-time evolution. Of particular interest is the role that $\delta$ plays in selecting which of these two states is observed. Recent numerical studies suggest that, after a transient period of rapid decay of the high Fourier modes, the bar state will be selected if $\delta \neq 1$, while the dipole will be selected if $\delta = 1$. Our results support this claim and seek to mathematically formalize it. We consider the system in Fourier space, project it onto a center manifold consisting of the lowest eight Fourier modes, and use this as a model to study the selection of bars and dipoles. It is shown for this ODE model that the value of $\delta$ controls the behavior of the asymptotic ratio of the low modes, thus determining the likelihood of observing a bar state or dipole after an initial transient period. Moreover, in our model, for all $\delta \approx 1$, there is an initial time period in which the high modes decay at the rapid rate $\mathcal{O}(e^{-t/\nu})$, while the low modes evolve at the slower $\mathcal{O}(e^{-\nu t})$ rate. The results for the ODE model are proven using energy estimates and invariant manifolds and further supported by formal asymptotic expansions and numerics.

Figures

Equations

$${\partial }_t\mathbf{u}=\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 .$$

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

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

$$\mathbf{u}_{bar}(x,t)=-\delta e^{-\frac{\nu}{\delta^{2}}t}\left(\begin{array}[]{c}0\\ \cos(x/\delta)\\ \end{array}\right),\qquad\mathbf{u}_{bar}(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),$$

$${\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}},$$

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

$$\omega (x,y)=\sum _{\vec{k}\ne 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\omega (x,y)e^{-\mathrm{i}(k_{1}x/\delta +k_{2}y)}\mathrm{d}x\mathrm{d}y.$$

$${\omega }_1$$

$${\omega }_5$$

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

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

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

$$+(d_{34}{\omega }_3{\omega }_4+\cdots +d_{38}{\omega }_3{\omega }_8)+(e_{45}{\omega }_4{\omega }_5+\dots e_{48}{\omega }_4{\omega }_8)+(f_{56}{\omega }_5{\omega }_6+\cdots +f_{58}{\omega }_5{\omega }_8)$$

$$+(g_{67}{\omega }_6{\omega }_7+g_{68}{\omega }_6{\omega }_8)+h_{78}{\omega }_7{\omega }_8+\Theta (3)$$

$$=:H({\omega }_1,\dots ,{\omega }_8;\vec{k}),$$

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

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

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

$$+\frac{\delta (3-{\delta }^2)}{2(1+{\delta }^2)}\hat{\omega }(2,0){\bar{\omega }}_7+\frac{1-3{\delta }^2}{2\delta (1+{\delta }^2)}\hat{\omega }(0,2){\omega }_7+\Theta (4)$$

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

$$+\frac{3{\delta }^2-1}{2\delta (1+{\delta }^2)}\hat{\omega }(0,-2){\omega }_5+\frac{\delta ({\delta }^2-3)}{2(1+{\delta }^2)}\hat{\omega }(2,0){\bar{\omega }}_5+\Theta (4)$$

$${\dot{\omega }}_{\vec{k}}=\nabla H({\omega }_1,\dots ,{\omega }_8;\vec{k})\cdot \left[-\nu \left(\frac{{\omega }_1}{{\delta }^2},\frac{{\omega }_2}{{\delta }^2},{\omega }_3,{\omega }_4,\frac{1+{\delta }^2}{{\delta }^2}{\omega }_5,\dots ,\frac{1+{\delta }^2}{{\delta }^2}{\omega }_8\right)+\Theta (2)\right].$$

$$-\frac{\nu }{{\delta }^2}{b}_{15}-\frac{\nu }{{\delta }^2}(1+{\delta }^2)b_{15}=-\frac{\nu }{{\delta }^2}(4+{\delta }^2)b_{15}-\frac{{\delta }^3}{1+{\delta }^2}\Rightarrow b_{15}=-\frac{{\delta }^5}{2\nu (1+{\delta }^2)}.$$

$$\hat{\omega }(2,1)$$

$$\hat{\omega }(2,-1)$$

$$\hat{\omega }(1,2)$$

$$\hat{\omega }(1,-2)$$

$$-2\nu \frac{(1+{\delta }^2)}{{\delta }^2}{a}_1=-\nu \frac{(4+{\delta }^2)}{{\delta }^2}{a}_1$$

$$\hat{\omega }(2,0)=G_{(2,0)}({\omega }_5,{\omega }_7)$$

$$+{\alpha }_7{\omega }_5{\omega }_7+{\alpha }_8{\omega }_5{\bar{\omega }}_7+{\alpha }_9{\bar{\omega }}_5{\omega }_7+{\alpha }_{10}{\bar{\omega }}_5{\bar{\omega }}_7+\Theta (3)$$

$$\hat{\omega }(0,2)=G_{(0,2)}({\omega }_5,{\omega }_7)$$

$$+{\beta }_7{\omega }_5{\omega }_7+{\beta }_8{\omega }_5{\bar{\omega }}_7+{\beta }_9{\bar{\omega }}_5{\omega }_7+{\beta }_{10}{\bar{\omega }}_5{\bar{\omega }}_7+\Theta (3)$$

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

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

Full text

1 Introduction

In this paper we 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$. Defining $\omega=(0,0,1)\cdot(\nabla\times\mathbf{u})$, one obtains the 2D vorticity equation

[TABLE]

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

Because the viscosity is small, it is reasonable to expect that stationary solutions of the Euler equation ((1.1) or (1.2) with $\nu=0$) would play a role in the long-time evolution of the Navier-Stokes equation. However, the Euler equation has infinitely many stationary solutions, so it is not obvious which such solutions are important. In [16], entropy arguments and extensive numerical studies were conducted in the case $\delta=1$ and suggested that the so-called bar states 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 state only emerged for a special class of initial data. A later study [3] numerically analyzed (1.1) on $D_{\delta}$ with the addition of a certain type of stochastic forcing. There, after an initial transient period, a 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$. Related analytical work was conducted in [2, 12] where the rate of convergence to a bar state for appropriate initial conditions was shown to be $\mathcal{O}(e^{-\sqrt{\nu}t})$, while the bar state itself decayed at the $\mathcal{O}(e^{-\nu t})$ background rate. In this work, we will analyze the selection of bars and dipoles, based on the parameter $\delta$. At this point, we also refer the interested reader to [1, 7, 8, 9, 10, 13] for steady state results and results in the asymptotic regime as time goes to infinity.

If $\delta=1$, any function of the form

[TABLE]

is an exact solution to (1.2). If $\delta\neq 1$, then (1.3) is an exact solution to (1.2) if and only if $a_{1}=a_{2}=0$ or if $a_{3}=a_{4}=0$. Bar states, also known as unidirectional or Komogorov flow, are members of this family for $m=1$ 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 for $m=1$ and are given by

[TABLE]

with velocity field

[TABLE]

or similarly with sine replaced by cosine. The bar states are exact solutions of (1.2) for all $\delta\approx 1$, while the dipoles are only exact solutions for $\delta=1$. In addition to the references mentioned above, the bar states were also studied analytically in [15]. Although the setting was slightly different, their results suggest that, when $\delta=1$ an $m$-bar state $e^{-\nu m^{2}t}\cos(my)$ (or similarly with sine replaced by cosine or $y$ replaced by $x$) is attracting if and only if $m=1$. Because the dipoles are only approximate solutions for $\delta\neq 1$, it may be intuitive that they would not play a key role in the long-time evolution in that case. However, they were still observed in the metastable switching in the appropriately stochastically forced Navier-Stokes equation for $\delta\neq 1$ [3].

Because of the form of the bar states and dipoles, it is useful to study (1.2) in Fourier space, in which it can be written

[TABLE]

where

[TABLE]

and

[TABLE]

In terms of these variables, 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, 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 dipoles correspond to solutions with energy in both the $\vec{k}=(0,\pm 1)$ and $\vec{k}=(\pm 1,0)$ modes. These four modes are the lowest modes in the system, in that they correspond to the modes with the lowest values of $|\vec{k}|_{\delta}$, with $|\vec{k}|_{\delta}=1$ or $\delta^{2}$. We will refer to any modes with $|\vec{k}|_{\delta}>\max\{1,\delta^{2}\}$ to be high modes.

When $\delta=1$, the set $\{\hat{\omega}_{\vec{k}}=0\mbox{ if }|\vec{k}|>1\}$ is an exact global invariant manifold for (1.4). However, the dynamics on it are trivial, determined by the linear terms. Therefore, even though both the bars and dipoles lie within this manifold, if we want to understand how the system selects between them, we must include at least some of the higher modes. To do so, we conduct a center manifold reduction on (1.4) and project onto the lowest eight modes, which we denote by

[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 (1.2) is real valued,

[TABLE]

Thus, the resulting ODE, which is derived in detail in §2, will be eight dimensional.

The reduction to the eight-dimensional ODE is local. In fact, since the size of the spectral gaps for the linear operator, $\nu\Delta$, is $\mathcal{O}(\nu)$, this reduction will only be valid in a small neighborhood of [math] of size $\mathcal{O}(\nu)$. Moreover, one cannot expect to obtain a finite-dimensional model of the full system (1.4) that describes the global dynamics [17]. However, we will still use this finite-dimensional model to provide insight into the potential role that $\delta$ plays in the selection of bars and dipoles. For other examples in which finite-dimensional models have been used to study the dynamics of the Navier-Stokes equation, see [5, 14].

The ODE derived in §2 will be analyzed in §3-§5. In §3 we focus on the case $\delta=1$, which corresponds to the symmetric torus, and in §4 we focus on the case $\delta\neq 1$. In both cases, to study the relative importance of the bar states versus the dipoles, we consider the evolution of the ratio $R(t):=\frac{|\omega_{1}(t)|^{2}}{|\omega_{3}(t)|^{2}}$. (In the case $\omega_{3}=0$, one can study the inverse of this quantity.) Note that asymptotic convergence of $R(t)\to 0$ or $\infty$ would correspond to convergence to a $y$- or $x$-bar state, respectively, while convergence to some finite, nonzero value would correspond to convergence to a dipole. (We note, however, that if $R(t)\to R_{\infty}\ll 1$ for example, then such a state would qualitatively appear to be a $y$-bar state, even though there would be nonzero variation in $x$.)

For the case $\delta=1$, in Theorem 3.4 we show that there is a family of co-dimension one stable manifolds in the phase space of the ODE that determines the asymptotic limit of $R(t)$. The limit $R(t)\to 0$ corresponds to exactly one of these manifolds, and hence a $y$-bar state would only be observed for the special class of initial conditions starting on this manifold. (A similar result holds for the $x$-bar states.) Therefore, we conclude that, for the symmetric torus, general initial conditions will typically lead to the emergence of a dipole as the dominant quasi-stationary state. For the case $\delta\neq 1$, the single center direction that had been present in the system for $\delta=1$ becomes hyperbolic, with the sign of $\delta-1$ determining if it is expanding or contracting. Thus, this selects the limit $R(t)\to 0$ or $R(t)\to\infty$, selecting an $y$-bar state or $x$-state respectively. These results are found in Theorem 4.4.

In both cases, $\delta=1$ and $\delta\neq 1$, we additionally show that the high modes decay at the rate $\mathcal{O}(e^{-t/\nu})$, while the low modes decay at the rate $\mathcal{O}(e^{-\nu t})$. These results can be found in Lemma 3.1 and Proposition 4.2. This allows for the rapid convergence to a metastable state as seen for the bar states in [2, 12]. We note however, that the rapid decay in [2, 12] was $\mathcal{O}(e^{-\sqrt{\nu}t})$, whereas here we obtain decay of the high modes at $\mathcal{O}(e^{-t/\nu})$. We expect that this discrepancy is due to the fact that the ODE we study is only a model of the full PDE. The main point is the qualitative prediction of a separation in time scales, rather than the specific rate. In addition, we comment here that the reason we do not scale the $\nu$ out of the equation

In §5, we reframe the problem as a perturbation problem, so as to confirm the results of the earlier sections using a different method. After setting $\delta=1+\epsilon_{0}\epsilon$, where $\epsilon_{0}=\pm 1$ will determine if $\delta$ is less than or greater than 1, relating the parameters $\nu$ and $\delta$ via $\epsilon$, and scaling the system in an appropriate way, a slow-fast system emerges. Perturbation expansions are then used to illustrate the claims made in §3 and §4. These expansions also reveal that evolution to a bar state accelerates as $\delta$ is moved slightly farther from 1. This result is consistent with the work of [3] in that the simulations done there suggest that a bar state dominates the metastable stochastic transitions only when $\delta$ is sufficiently far from 1.

Lastly, we conclude the introduction by pointing out that one could in principle scale time and velocity from the very beginning in order to make $\nu=1$. However, our goal here is to investigate the selection mechanism based on the value of $\delta$ and it turns out that there is an interplay with $\nu$ as well. As a matter of fact, estimates like those in Lemma 3.1 and Proposition 4.2 depend on $\nu$ in very precise ways and if one had scaled out $\nu$ in the beginning, one would have to undo the scalings later on. Hence, we chose to keep the dependence on $\nu$ as it is given by the equation originally.

2 Center Manifold

In this section, we carry out a center manifold reduction of (1.4) onto the eight modes listed in (1.5). This is a standard calculation that can be found, for example, in [11].

The basic idea is, for any $\omega_{\vec{k}}$ with $\vec{k}\notin\{(\pm 1,0),(0,\pm 1),(\pm 1,\pm 1)\}=:\mathcal{K}_{0}$, to assume 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 dynamics of (1.4). 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}$. To make this precise, define

[TABLE]

where $\Theta(3)$ is defined in Definition 2.1.

Definition 2.1**.**

Let $n\in\mathbb{N}$ and let $\omega_{1},\dots,\omega_{8}$ be defined as in (1.5). We define $\Theta(n)$ to be the set of terms that are of the form $\omega_{1}^{n_{1}}\omega_{2}^{n_{2}}...\omega_{8}^{n_{8}}$ where $n_{1}+...+n_{8}\geq n$ for $n_{1},\dots,n_{8}\in\mathbb{N}$. Thus, $\Theta(n)$ is the set of monomials in $\omega_{1},\dots,\omega_{8}$ of degree $n$ or larger.

In order to determine the values of the unknown coefficients that appear in (2.1), we will compute $\partial_{t}\omega_{j}$, for $j=1,\dots,8$ in two ways: 1) by using the right hand side of equation (1.4), where we substitute $H(\cdot,\vec{k})$ in for $\omega_{\vec{k}}$ whenever $\vec{k}\notin\mathcal{K}_{0}$; and 2) by computing the time derivative of the expansion in (2.1), substituting in the equation for $\dot{\omega}_{j}$, $j=1,\dots,8$ given in (1.4) as necessary. Equating these two results, and comparing terms with equal order in powers of $\omega_{1},\dots,\omega_{8}$, will lead to equations that should in theory determine the values of the coefficients. Note that the coefficients in (2.1) will depend on $\vec{k}$, but we have suppressed this for notational convenience.

We are only interested in computing the expansion for $H(\cdot,\vec{k})$ up to and including terms of $\Theta(2)$. Therefore, when carrying out the above-described calculation, we will only need to retain terms up to any including $\Theta(3)$. This means any product of the form $\omega_{\vec{k}}\omega_{\vec{j}}$ with $\vec{k},\vec{j}\notin\mathcal{K}_{0}$ will be of higher order, and thus we can discard it. Carrying out step 1) above, leads to

[TABLE]

Note we have listed only four of the equations, due to (1.6). Therefore, we need only focus on determining the coefficients of $H(\cdot,\vec{k})$ for $\vec{k}\in\mathcal{K}:=\{(\pm 1,\pm 2),(\pm 2,\pm 1),(\pm 2,0),(0,\pm 2)\}$. Carrying out step 2) leads to

[TABLE]

Equating (2.2) and (2.3) leads to, for example, the following for the $b_{15}$ coefficient of $H(\omega_{1},\dots,\omega_{8};(2,1))$:

[TABLE]

Continuing in this manner, we find

[TABLE]

Conveniently, most of the coefficients are zero. In order for the above equations to be the unique expansion for the $H$’s that satisfy the invariance condition, we need to restrict the values of $\delta$ to be sufficiently close to $1$. To see why this is the case, consider, for example, the $a_{1}$ coefficient of $\vec{k}=(2,1)$. Its defining equation is given by

[TABLE]

If $\delta=\sqrt{2}$, then $a_{1}=0$ is not the unique solution. To prevent such ambiguities and obtain a unique expansion for $\hat{\omega}(2,1)$, $\hat{\omega}(2,-1)$, $\hat{\omega}(1,2)$, and $\hat{\omega}(1,-2)$, it turns out that we should restrict the value of $\delta$ to $(\frac{1}{\sqrt{2}},\sqrt{2})$.

Computing the expansions for $\hat{\omega}(\pm 2,0)$ and $\hat{\omega}(0,\pm 2)$ turns out to be somewhat different. When $\delta\neq 1$, one can check that, if $\sqrt{\frac{2}{3}}<\delta<\sqrt{\frac{3}{2}}$, the coefficients are all unique and equal to 0. However, when $\delta=1$, a unique set of coefficients cannot be determined. In this case, for any constants $\alpha_{i}$ and $\beta_{i}$, the following functions will satisfy the invariance condition

[TABLE]

Note that we have relabeled these functions using the letter $G$, rather than $H$ as above, to emphasize that this is a special case only when $\delta=1$. Moreover, we have used (1.6) to write these as functions of $\omega_{5,7}$ only, for notational convenience and to highlight the fact that they depend only on the high modes. Similar equations for $\hat{\omega}(-2,0)$ and $\hat{\omega}(0,-2)$ can be found, based on (2.5), using (1.6). We will comment more on this issue of nonuniqueness in the $\delta=1$ case in Remark 2.2, below. Thus, we arrive at the following eight-dimensional ODE model, for $\delta\neq 1$ and $\sqrt{\frac{2}{3}}<\delta<\sqrt{\frac{3}{2}}$

[TABLE]

and for $\delta=1$ we find the following eight-dimensional ODE model

[TABLE]

Remark 2.2**.**

We suspect that the nonuniqueness of the expansions in the $\delta=1$ case is related to the fact that the above calculation is only local and valid in a small neighborhood of $\mathcal{O}(\nu)$. In order to obtain continuity of our ODE model for $\delta\approx 1$, it makes sense to chose $G_{(0,2)}=G_{(2,0)}=0$ in the $\delta=1$ case. This makes (2.6) equal to (2.7) in the limit $\delta\to 1$. Therefore, we make this choice in the following sections, and consider only equation (2.6) for all $\delta\in(\sqrt{\frac{2}{3}},\sqrt{\frac{3}{2}})$.

We conclude this section with a quick observation that reflects a symmetry in (1.4) and that can be used to simplify some of the proofs in the following sections. Its proof is omitted as it is relatively straighforward.

Lemma 2.3**.**

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 for (2.6).

3 Symmetric Torus

We now focus on (2.6) with $\delta=1$, which corresponds to the symmetric torus. The goal will be to prove two results: Lemma 3.1, which states that the high modes decay much more rapidly than the low modes, and Theorem 3.4, which states that most initial conditions will evolve toward a dipole, rather than a bar state. We begin by rewriting system (2.6) with $\delta=1$:

[TABLE]

It will be helpful to study the evolution of the quantities $A:=|\omega_{1}|^{2}+|\omega_{3}|^{2}$ and $B:=|\omega_{5}|^{2}+|\omega_{7}|^{2}$, in order to separate the evolution of the low modes from that of the high modes, and the quantity $R(t)=\frac{|\omega_{1}(t)|^{2}}{|\omega_{3}(t)|^{2}}$ to study whether it is a bar state or dipole that can be expected to be observed for large time. Recall that $R(t)\to 0$ corresponds to evolution towards a $y$-bar state, $R(t)\to\infty$ corresponds to evolution to an $x$-bar state, and $R(t)\to R_{\infty}\in(0,\infty)$ corresponds to convergence to a dipole, as $t\to\infty$.

Simulations of (3.1) shown in Figure 1 below suggest that there is a separation in time scales between the evolution of $A$ and $B$, consistent with previous numerical studies of the full Navier-Stokes equation [16]. Moreover, $A$ exhibits an initial period of growth, before it begins to decay. Furthermore, in Figure 2 we see that a variety of initial values for $R$ lead to solutions that converge to a dipole. These behaviors will be made more precise below, in Lemma 3.1 and Theorem 3.4.

Lemma 3.1**.**

Define $A(t)$ := $|\omega_{1}(t)|^{2}+|\omega_{3}(t)|^{2}$ and $B(t):=|\omega_{5}(t)|^{2}+|\omega_{7}(t)|^{2}$. Let $t_{0}=1/\nu$, $\delta=1$, and denote the initial data by $A(0)=A_{0}$ and $B(0)=B_{0}$. We have

[TABLE]

Moreover, for all $0\leq t\leq t_{0}$, $A(t)\geq A_{0}e^{-2}$ and $B(t)\leq B_{0}e^{-\frac{2A_{0}}{5\nu e^{2}}t}$. Finally, for all $t\geq t_{0}$, $B(t)\leq B_{0}e^{-\frac{2A_{0}}{5\nu^{2}e^{2}}}$.

Proof.

Using (3.1) we find that the dynamics of $A$ and $B$ are governed by

[TABLE]

The first claim follows from the fact that, since $A$ and $B$ are both nonnegative,

[TABLE]

Also, since $A$ and $B$ are nonnegative, $\dot{A}\geq-2\nu A$, and so $A(t)\geq A_{0}e^{-2\nu t}$. As a result, for all $0\leq t\leq t_{0}$, $A(t)\geq A_{0}e^{-2}$. With this estimate, we then see that for all $0\leq t\leq t_{0}$

[TABLE]

from which the claim about $B$ during this time period follows by Gronwall’s Inequality. Finally, note that $\dot{B}\leq 0$ so $B$ is monotonically decreasing (in fact strictly decreasing as long as $B>0$), which proves the final claim. ∎

Remark 3.2**.**

Regarding the timescales predicted by the above lemma, recall that the center manifold reduction is expected to be valid only in an $\mathcal{O}(\nu)$ neighborhood of $\omega=0$. If we correspondingly scale $A=\nu^{2}\tilde{A}$ and $B=\nu^{2}\tilde{B}$, then (3.2) becomes

[TABLE]

Note, however, that despite this rescaling with $\nu$, and due to the difference in sign of the nonlinear terms, there will (for most initial data) be an transient period when $\tilde{A}$ grows, while $\tilde{B}$ decays, before both $\tilde{A}$ and $\tilde{B}$ decay to zero. We thank an anonymous referee for pointing this out.

This proposition shows that, when $\nu$ is small, for a long $\mathcal{O}(1/\nu)$ transient period, $B(t)$ is decaying at the rapid rate $\mathcal{O}(e^{-t/\nu})$, while $A(t)$ is not changing much. After this transient period, $B$ has become exponentially small and both $A$ and $B$ decay to zero at the background rate, $\mathcal{O}(e^{-2\nu t})$. Figures 3 and 4 below show the evolution of $A(t)$ and $B(t)$ on a logarithmic scale. These plots elucidate the results of Lemma 3.1. We see that following the initial period of growth, the curve L=$\log(A(t))$ decreases at a rate very close to $L=-2\nu t$, suggesting that the low modes indeed decay no faster than the global background decay rate. The particular line, $L=\max(A(t))-2\nu t$, that was graphed was an asthetic choice to clearly show the two curves become parallel. The trend can be observed to continue for longer times. Figure 4 illustrates that for $0<t<\frac{1}{\nu}=100$, $L=\log(B(t))$ initially decreases faster than the line $L=-\frac{2A_{0}}{5\nu e^{2}}t$, followed by less rapid decay. The point of intersection at $t=120$ was a choice made to most clearly show the change in decay rate of the higher modes once $t>\frac{1}{\nu}$. An analogous result will be shown for $\delta\neq 1$ in §4 in Lemma 4.1 and Proposition 4.2.

Remark 3.3**.**

Interestingly, one can check that if the general forms of $G_{(2,0)}$ and $G_{(0,2)}$ given by (2.5) were used in (3.1), $\dot{A}$ and $\dot{B}$ as defined in the statement of Lemma 3.1 would still satisfy (3.2). Therefore, the dynamics of $A$ and $B$ are exactly the same regardless of the choice of coefficients for $G_{(2,0)}$ and $G_{(0,2)}$.

Next, to precisely describe the behavior exhibited in Figure 2, we define the following variables (see also Remark 3.5):

[TABLE]

Equation (3.1) then implies that

[TABLE]

where P_{re}=\mathop{\mathrm{missing}}{Re}\nolimits(P), P_{im}=\mathop{\mathrm{missing}}{Im}\nolimits(P), Q_{re}=\mathop{\mathrm{missing}}{Re}\nolimits(Q), and Q_{im}=\mathop{\mathrm{missing}}{Im}\nolimits(Q).

Theorem 3.4**.**

There exists a family of local stable manifolds, $\mathcal{M}_{r}=\{R=f(A,w,z,P_{re},P_{im},Q_{re},Q_{im};r)\}$ for $r\geq 0$, corresponding to each of the fixed points of the system (3.4). Any initial condition on $\mathcal{M}_{r}$ will converge at an exponential rate to the fixed point $\{R=r,A=0,w=0,z=0,P=0,Q=0\}$ as $t\to\infty$. This family, up to and including quadradic terms, is given by

[TABLE]

Proof.

The existence of this family of stable manifolds is given by Theorem 4.1 in [4]. Before applying this theorem, it must be verified that the nonlinearity of (3.4), as well as its partial derivatives, vanish at each fixed point. To see that no singularities exist, note that the domain of (3.4) is determined by the way the variables are defined in (3.3). It is given by the set

[TABLE]

The vector field in (3.4) has no singular points in this domain since,

[TABLE]

Therefore, at any fixed point of the form $\vec{r}=(r,0,0,0,0,0,0,0)$, it can be checked that the nonlinearity and all of its partial derivatives vanish. Thus we obtain the existence of a local invariant (stable) manifold $\mathcal{M}_{r}$, associated to the flow toward the fixed point at $\vec{r}$. The rest of the proof is a long, but standard, computation of the family $\mathcal{M}_{r}$ and will be left to the appendix. ∎

The implication Theorem 3.4 has on the originial system (3.1) for $\omega_{1},$ $\omega_{3},$ $\omega_{5},$ and $\omega_{7}$ is that for a typical initial condition, system (3.1) will not evolve to an $x$-bar ($r=\infty$) or $y$-bar ($r=0$) state. For most $r\in(0,\infty)$, the asymptotic state will look like a dipole. We note, however, that for $r\ll 1$ or $r\gg 1$ the asymptotic state may appear to be more like a bar state than a dipole, even though it is not a pure bar state. This is again consistent with the numerical results of [3, 16].

Remark 3.5**.**

Theorem 3.4 establishes the existence of locally invariant manifolds only for initial conditions to (3.1) with $|\omega_{3}|\neq 0$. The same result can be proven for initial conditions to (3.1) with $|\omega_{1}|\neq 0$ by taking instead as the change of variables $\mathcal{R}=\frac{|\omega_{3}|^{2}}{|\omega_{1}|^{2}}$, $\mathcal{P}=\frac{\omega_{1}\bar{\omega}_{3}\bar{\omega}_{7}}{|\omega_{1}|^{2}}$, and $\mathcal{Q}=\frac{\bar{\omega}_{1}\bar{\omega}_{3}\omega_{5}}{|\omega_{1}|^{2}}$ with $A$, $w$, and $z$ unchanged. The set $\{\omega_{1}=\omega_{3}=0\}$ is invariant for (3.1), so this case is not relevant for the study of the selection of the bars versus dipoles.

Remark 3.6**.**

The invariant manifolds $\mathcal{M}_{r}$ in Theorem 3.4 are defined locally near each fixed point $\vec{r}$. However, Lemma 3.1 implies that, for any initial condition, the solution will eventually enter a neighborhood of this one-dimensional manifold of fixed points, and hence its dynamics will eventually be governed by the theorem.

Remark 3.7**.**

If we instead use the more general form of $\hat{\omega}(2,0)$ and $\hat{\omega}(0,2)$ on the center manifold given by (2.5), one can check that, although the explicit form of $\mathcal{M}_{r}$ would change, qualitatively the same structure will prevail. There remains a line of fixed points as the one-dimensional center manifold and a family of stable manifolds consisting of the initial conditions that evolve towards each of these fixed points.

4 Asymmetric Torus

In this section we study (2.6) with $\delta\neq 1$. In particular, we show in Proposition 4.2 that again the high modes decay much more rapidly than the low modes, and we show in Theorem 4.4 that $R(t)\to 0$ as $t\to\infty$ for $\delta<1$, thus indicating convergence to a $y$-bar state. A similar convergence result to an $x$-bar state exists when $\delta>1$. We begin by rewriting (2.6) for convenience:

[TABLE]

For the remainder of this section, we will assume that $\delta$ is fixed with $\delta\in(\sqrt{\frac{2}{3}},\sqrt{\frac{3}{2}})$, $\delta\neq 1$. It will occasionally be useful to fix a sufficiently small value of $\eta>0$ and require that $|\delta^{2}-1|<\eta$. For such a fixed $\delta$ and $\eta$, some of our results will hold for all $\nu>0$ sufficiently small.

First, we obtain a global decay rate via the energy

[TABLE]

Lemma 4.1**.**

For all $t\geq 0$, $E(t)\leq E(0)e^{-2K_{1}\nu t}$, with $K_{1}=\min({1,1/\delta^{2}})$.

Proof.

Using (4.1) we see that

[TABLE]

which proves the result. ∎

To obtain a faster decay rate for the higher modes, we again define $A=|\omega_{1}|^{2}+|\omega_{3}|^{2}$ and $B=|\omega_{5}|^{2}+|\omega_{7}|^{2}$, and denote the real and imaginary parts of $\omega_{i}$ as $\omega_{i}^{re}$ and $\omega_{i}^{im}$ respectively for $i=1,3,5,7.$ We find using (4.1) that

[TABLE]

.

Proposition 4.2**.**

Let $A_{0}=A(0)$, $B_{0}=B(0)$, and $K_{1}=\min(1,1/\delta^{2})$. For any fixed $\delta$ and $\eta$ with $0<|\delta^{2}-1|<\eta$, there exists a $\nu^{*}>0$ such that for all $0<\nu<\nu^{*}$ there exists a time $t_{*}=\mathcal{O}(\nu|\log(\nu\eta)|)$ and positive constants $M_{0}$, $M_{1}$ and $K_{2}$, independent of $\nu$, such that

[TABLE]

Moreover, for all $0\leq t\leq\min(1/\nu,1/\eta)$, $A(t)\geq A_{0}e^{-K_{2}}$, and for all $t\geq 0$, $A(t)+B(t)\leq(A_{0}+B_{0})e^{-2\nu K_{1}t}$.

Proof.

By Lemma 4.1, we have $A(t)+B(t)\leq(A_{0}+B_{0})e^{-2\nu K_{1}t}$. Next, notice that

[TABLE]

Using the fact that $|\omega_{1}\omega_{3}|\leq(|\omega_{1}|^{2}+|\omega_{3}|^{2})/2=A/2$, we have

[TABLE]

Thus, for all $t\geq 0$,

[TABLE]

If we now let $t_{1}=\min(1/\nu,1/\eta)$, then $A(t)\geq A_{0}e^{-K_{2}}$ for all $0\leq t\leq t_{1}$, where

[TABLE]

Here we used the assumption that $\delta\geq\sqrt{\frac{2}{3}}$. To obtain the estimates for $B$, notice that the above equation for $\dot{B}$ implies

[TABLE]

where

[TABLE]

Since $\delta\in(\sqrt{\frac{2}{3}},\sqrt{\frac{3}{2}})$, $D_{0}$ is bounded away from 0. Suppose first that

[TABLE]

(If the initial condition $B_{0}<B_{*}$, then it is already asymptotically small and we will consider this case afterwards.) Then we find that

[TABLE]

Thus, setting $M_{0}=(D_{0}A_{0})/(2e^{K_{2}})$ we find that $B(t)\leq B_{0}e^{-\frac{M_{0}}{\nu}t}$, at least until $B(t)=B^{*}$. The latter occurs at a time no larger than $t_{*}$, where $t_{*}$ is defined via

[TABLE]

Finally, consider times for which $B(t)\leq B_{*}$. By Lemma 4.1, $B(t)\leq B_{0}e^{-2\nu K_{1}t}$ for all time. Hence, if we restart the system at time $t=t_{*}$, when $B=\mathcal{O}(\eta^{2}\nu^{2})$, then $B(t)$ decays at the background rate for $t\geq t_{*}$, i.e. $B(t)\leq M_{1}\eta^{2}\nu^{2}e^{-2\nu K_{1}t}$. This completes the proof. ∎

Remark 4.3**.**

We comment on the relationship between Proposition 4.2 and Lemma 3.1. Lemma 3.1 is a result in the case $\delta=1$. In Proposition 4.2, one must begin by fixing $\delta$ and $\eta$ so that $0<|\delta^{2}-1|<\eta$. Although one can take $\delta$ arbitrarily close to $1$, but fixed, and hence $\eta$ arbitrarily close to zero, but fixed, one cannot take the limits $\delta\to 1$ and $\eta\to 0$ and obtain something meaningful from Proposition 4.2. Moreover, although $t_{*}\to 0$ as $\nu\to 0$, one must also be careful with this limit, because the bound on $|B(t)|$ also goes to zero in this limit: if we were to choose $t=t_{*}/2$ in this bound, we would have

[TABLE]

which goes to zero as $\nu\to 0$. Thus, one really must keep $\delta,\eta$, and $\nu$ fixed when considering Proposition 4.2.

We further note that the results of Proposition 4.2 do not contradict Lemma 3.1. Both Lemma 3.1 and Proposition 4.2 say that for an initial period of time, $B(t)$ is bounded above by a term of the order $O(e^{-\frac{\Lambda}{\nu}t})$, albeit with slightly different constants $\Lambda>0$ in each one of the two cases. After that initial period of time, both results say that $B(t)$ remains small, and to be precise in Lemma 3.1 we showed that for $t$ large enough $B(t)$ is bounded above by $e^{-\frac{2A_{0}}{5\nu^{2}e^{2}}}$, whereas in Proposition 4.2 we can prove an upper bound of the order $\eta^{2}\nu^{2}e^{-2\nu K_{1}t}$. In the case $\delta\neq 1$, we can only guarantee the upper bounds that appear in Proposition 4.2, which however are sufficient for our purposes, as Theorem 4.4 below demonstrates.

Now consider again the variable, $R=\frac{|\omega_{1}|^{2}}{|\omega_{3}|^{2}}$. Suppose $\delta<1$ and set $\gamma=2\nu(\frac{1}{\delta^{2}}-1)\leq C\nu\eta$. We wish to show that $R(t)\sim\mathcal{O}(e^{-\gamma t})$ for all $t\geq 0$, as long as $R(0)=:R_{0}$ is not too big. Using (4.1), we find

[TABLE]

Theorem 4.4**.**

For any fixed $\delta<1$ , sufficiently close to 1, and fixed $\eta$ sufficiently small, with $|\delta^{2}-1|<\eta$, there exist $R^{*},\nu^{*}>0$ such that, for all $0<\nu\leq\nu^{*}$ and $R_{0}\leq R^{*}$, $R(t)\leq M_{2}e^{-\gamma t}$ for some $M_{2}>0$ and all $t\geq 0$. $R_{*}$ is $\mathcal{O}(\min(\eta^{-1},\nu^{-1}))$ while $M_{2}$ is independent of $\nu$ and $\eta$.

Proof.

For ease of notation, define

[TABLE]

With this notation,

[TABLE]

Therefore, we have the bound

[TABLE]

Define

[TABLE]

where $T$ is defined to the largest time such that $e^{\gamma t}R(t)\leq R_{*}$, and $R_{*}$ will be determined below. If $|||R|||$ is finite, then $R(t)$ decays like $e^{-\gamma t}$. Multiplying the integral inequality above by the exponential weight $e^{\gamma t}$, we find that

[TABLE]

where

[TABLE]

We now estimate the three terms above using Proposition 4.2, by splitting the time interval into two pieces: $0\leq t\leq t_{*}$ and $t_{*}\leq t$. So Proposition 4.2 gives for term I,

[TABLE]

for some constant $C_{1}$ that is independent of $\nu$ and $\eta$. Similarly, for term II we have

[TABLE]

for some constant $C_{2}$ that is independent of $\nu$ and $\eta$. Finally, for term III we have

[TABLE]

for some constant $C_{3}$ that is independent of $\nu$ and $\eta$. Combining the above estimates, we find that

[TABLE]

We can rewrite this as

[TABLE]

Now, choose $R_{*}=1/(2\mbox{max}(\eta,\nu)C_{2})$. We then find

[TABLE]

as long as $\nu^{*}$ and $\eta$ are such that $\max(\eta,\nu^{*})\leq\frac{1}{4C_{1}}$. Using the above definition of $R^{*}$, the above right hand side will be less than $R^{*}$ if $R_{0}\leq R^{*}/8$ and $\max(\eta,\nu^{*})\leq 1/\sqrt{16C_{0}C_{2}}$. In this case,

[TABLE]

Because the bound on the right hand side is independent of $T$, the bound must hold for all time. Hence, there must exist an $M_{2}>0$ such that $R(t)\leq M_{2}e^{-\gamma t}$ for all $t\geq 0$. ∎

We emphasize again the importance of the transient period of rapid decay present in the dynamics of the higher Fourier modes established in Proposition 4.2. Without their rapid decay to a small enough order during this initial time period, the estimates in the proof of Theorem 4.4, particularly for term I, would not have gone through. Moreover, we see that given a small, fixed distance of $\delta$ from 1, a sufficiently small value for the viscosity can be selected to separate the decay rates into the two regimes established in Proposition 4.2 and used in the proof of Theorem 4.4, thus driving the system toward a bar state.

Let us now define $U(t)=R(t)^{-1}$. We see that $U(t)$ must satisfy

[TABLE]

The dynamics of $U$ for $\delta>1$ are analogous to the dynamics of R when $\delta<1$. With similar estimates to those in the proof of Theorem 4.4, one can show that $U(t)\to 0$ as $t\to\infty$, which indicates convergence to an $x$-bar state.

5 Perturbation Analysis

The purpose of this section is to provide further evidence of a selection mechanism through an alternate method. The system’s domain will be viewed as a perturbation of the symmetric torus. Using the approximations computed in this section, we will confirm the results of §3 and §4. To ultimately view (2.6) as a perturbed system, we define the perturbation parameter $\epsilon$ via $\delta=1+\epsilon_{0}\epsilon$, with $\epsilon_{0}=\pm 1$. Note that the sign of $\epsilon_{0}$ determines whether $\delta$ is greater or less than one. We begin by scaling $\nu$ as an appropriate power of $\epsilon$, which effectively relates $\nu$ and $\delta$. Subsequently, asymptotic expansions in $\epsilon$ are computed that are connected with the observed multiple time scales in the evolution of the vorticity. The properties of these expansions agree with the results of §3 and §4. Moreover, we are again able to observe the evolution to the appropriate bar state depending on if $\epsilon_{0}=\pm 1$. These expansions also have the property, at least among the $\mathcal{O}(1)$ and $\mathcal{O}(\epsilon)$ terms, that the emergence of a bar state will become faster as $|\delta-1|=\epsilon$ increases. This could be related to the observation in [3] that the bar states do not dominate the stochastic evolution (based on the stochastic forcing used there) unless $|\delta-1|$ is sufficiently large.

Motivated by geometric singular perturbation theory, we first scale (2.6) in a way that reveals a slow and a fast subsystem. We must choose appropriate values for the scaling parameters that accomplish two main objectives. First, to reveal a slow-fast system, we aim for the leading order terms in the scaled versions of $\dot{\omega}_{1}$ and $\dot{\omega}_{3}$ to be some order of magnitude in $\epsilon$ higher than those in the scaled versions of $\dot{\omega}_{5}$ and $\dot{\omega}_{7}$. Second, to ensure that the decay rates in the newly scaled system match those seen in the previous sections, we would like the leading order terms, once scaled, to match the terms in (5.1) below.

[TABLE]

If these correspond to the leading order terms in the scaled system, then the decay rates seen in the asymptotic expansions that are to be computed will match those observed in previous sections.

First scale the viscoscity and time by $\nu=\epsilon^{\alpha}\nu_{0}$ and $\tau=\epsilon^{\alpha}t$. Then scale the Fourier modes by $\omega_{1}=\epsilon^{\beta}\Omega_{1}$, $\omega_{3}=\epsilon^{\beta}\Omega_{3},$ $\omega_{5}=\epsilon^{\phi}\Omega_{5}$, and $\omega_{7}=\epsilon^{\phi}\Omega_{7}$. The acceptable range of values for the three scaling parameters $\alpha$, $\beta$, and $\phi$ that are relevant to this discussion will now be identified. The clearest way to do so will be to define $\rho$ and $\sigma$ as

[TABLE]

With these initial scalings, the system (5) below is obtained. Note here that the coefficient $\frac{\delta^{2}-1}{\delta}$, for $\delta=1+\epsilon_{0}\epsilon$ is $\mathcal{O}(\epsilon)$ as $\epsilon\to 0$.

[TABLE]

Looking closely at the powers of $\epsilon$ appearing in (5), one notices that if the following restrictions on $\sigma$ and $\rho$ hold, then the leading order terms are of the desired form given by (5.1):

[TABLE]

Remark 5.1**.**

As mentioned in the opening paragraph of this section, we aim to compute asymptotic expansions in $\epsilon$ for $\Omega_{1}$, $\Omega_{3}$, $\Omega_{5}$, and $\Omega_{7}$. If $\beta$ and $\phi$ satisfy the above inequalities, then the expansions computed have the expected properties for any value of $\alpha$. No restrictions on $\alpha$ frees us from being constrained to particular relative values of the viscocity and aspect ratio of the domain. Namely, after connecting the viscocity and aspect ratio through $\nu=\epsilon^{\alpha}\nu_{0}$, the value of $\alpha$ determines which of $\nu$ or $\epsilon$ is larger. Therefore the following results will hold for small values of $\nu$ and $\epsilon$ regardless of which is bigger relative to one another.

Proceeding with the computation of the asymptotic expansions, for simplicity, set $\phi=\alpha$ and $\beta=\alpha-\frac{1}{2}$ for any value of $\alpha>1/2$. With this choice of parameter values, the final scaled system that we work with from this point on can be obtained by substituting $\delta=1+\epsilon_{0}\epsilon$ into (5) and using the Taylor series expansions of the $\delta$ dependent coefficients. The resulting system is given below by (5.3).

[TABLE]

where

[TABLE]

[TABLE]

and

[TABLE]

With these scalings, it is evident that $\Omega_{1}$ and $\Omega_{3}$ are the slow variables, evolving with respect to $\tau$, to leading order, at an $\mathcal{O}(1)$ rate as $\epsilon\to 0$, while $\Omega_{5}$ and $\Omega_{7}$ evolve on the faster time scale $\mathcal{O}(\epsilon^{-1})$ as $\epsilon\to 0$. This matches the time scale separation we saw present in the previous sections, which can be seen below by reversing the scalings in the leading order terms.

[TABLE]

We now proceed using methods from geometric singular perturbation theory. Setting $\epsilon=0$ in (5.3) leads to the following leading order slow dynamics

[TABLE]

and so the leading order slow manifold is $M_{0}=\{\Omega_{5}=\Omega_{7}=0\}$. Observe that in the perturbed system (5.3), this manifold of fixed points for the $\epsilon=0$ reduced slow system is no longer invariant. This can be seen in the differential equations for $\Omega_{5}$ and $\Omega_{7}$. However, since $M_{0}$ is a normally hyperbolic manifold and the vector field in (5.3) satisfies the smoothness conditions of Fenichel’s theorems, a perturbed invariant manifold, $M_{\epsilon}$, exists for sufficiently small $\epsilon>0$ and is $\mathcal{O}(\epsilon)$ close to $M_{0}$. See Theorem 9.1 in [6]. Any trajectory in phase space will approach this manifold exponentially fast and then track the slow dynamics on $M_{\epsilon}$.

Defining the fast variable $s=\tau/\epsilon$ and setting $\epsilon=0$ in (5.3), we find the leading order fast dynamics to be given by

[TABLE]

Assuming expansions of $\Omega_{i}(s)$ for $i=1,3,5,7$ to be of the form $\Omega_{i}(s)=\Omega_{i0}(s)+\epsilon\Omega_{i1}(s)+\mathcal{O}(\epsilon^{2})$ away from the slow manifold, we find

[TABLE]

Here we see that away from the slow manifold, to leading order, the higher order modes are decaying at a rate $\mathcal{O}\left(e^{-\frac{|\Omega_{10}(0)|^{2}+|\Omega_{30}(0)|^{2}}{5\nu_{0}}s}\right)$, while the lowest modes are constant. This is consistent with the initial rapid decay rates among the higher modes seen in the previous sections, as seen in the calculation below:

[TABLE]

The dynamics on the slow manifold will determine whether solutions evolve towards a bar state or a dipole. To analyze this, we consider again system (5.3) and compute formal asymptotic expansions of the solutions in terms of the slow variable $\tau$. We begin by writing

[TABLE]

Next, we compute terms in these expansions up to and including $\mathcal{O}(\epsilon)$ terms. To do so, first consider the $\mathcal{O}(\epsilon^{-1})$ terms present in (5.3) and match these terms with the derivatives taken in (5.4). The only terms present are

[TABLE]

and so we find $|\Omega_{50}|^{2}=|\Omega_{70}|^{2}=0$. Next, matching the $\mathcal{O}(1)$ terms, we obtain the following system governing the dynamics of the $\mathcal{O}(1)$ terms of $\Omega_{1}$ and $\Omega_{3}$, as well as algebraic equations that determine the $\mathcal{O}(\epsilon)$ terms of $\Omega_{5}$ and $\Omega_{7}$:

[TABLE]

Solving these we obtain the following expressions:

[TABLE]

Yet to be computed are the $\mathcal{O}(\epsilon)$ terms for the lower modes, $\Omega_{1}$ and $\Omega_{3}$. The relevant equations are

[TABLE]

Using (5.5) to solve (5.6) we obtain

[TABLE]

Together, equations (5.5) and (5.7) make the approximations to $\Omega_{1}$ and $\Omega_{3}$ up to and including $\mathcal{O}(\epsilon)$, which will be denoted by $\bar{\Omega}_{1}$ and $\bar{\Omega}_{3}$:

[TABLE]

Of interest will be the magnitudes of the scaled low modes, which will be defined by $X:=|\Omega_{1}|^{2}$ and $Y:=|\Omega_{3}|^{2}$. In computing these, we obtain

[TABLE]

where we have used the notation $X(\tau)=X_{0}(\tau)+\epsilon X_{1}(\tau)+\mathcal{O}(\epsilon^{2})$ and $Y(\tau)=Y_{0}(\tau)+\epsilon Y_{1}(\tau)+\mathcal{O}(\epsilon^{2})$. For notational convenience, let $K_{x_{0},y_{0}}=\frac{20X_{0}(0)Y_{0}(0)}{X_{0}(0)+Y_{0}(0)}$.

Now take $\bar{X}:=X_{0}+\epsilon X_{1}$ and $\bar{Y}:=Y_{0}+\epsilon Y_{1}$ to be the approximations to $|\Omega_{1}|^{2}$ and $|\Omega_{3}|^{2}$. For each $0<\epsilon\ll 1$, there exists a finite interval of time on which these approximations are valid. Define $\tau_{+}$ and $\tau_{-}$ as follows

[TABLE]

Recall that, as approximations to the nonnegative quantities $|\Omega_{1}|^{2}$ and $|\Omega_{3}|^{2}$, $\bar{X}$ and $\bar{Y}$ must too be nonnegative. Observe that for $\epsilon_{0}=1$, $\bar{X}(\tau_{+})=0$ and $\bar{X}(\tau)<0$ on $\tau>\tau_{+}$. For $\epsilon_{0}=-1$, this property is shared by $\bar{Y}$ on $\tau\geq\tau_{-}$. These properties indicate that the approximations are certainly not valid for values of $\tau$ beyond $\tau_{\pm}.$

As a direct consequence of these observations, on a finite time interval $\epsilon_{0}$ will determine the bar state toward which the system evolves. This is summarized by Proposition 5.2 below and is consistent with the results in §4. The key property of the approximations $X(\tau)$ and $Y(\tau)$ in (5.8) is the linear growth exhibited by the $\mathcal{O}(\epsilon)$ terms. In particular, the opposite sign on the $\epsilon_{0}\nu_{0}\tau e^{-2\nu_{0}\tau}$ terms allows for the evolution toward the correct bar state. This linear growth originates from the resonant forcing terms in the differential equations for $\Omega_{11}$ and $\Omega_{31}$, which can be seen after substituting the expressions for $\Omega_{10}$ and $\Omega_{30}$ given in (5.5) into (5.6).

Proposition 5.2**.**

Let $0<\epsilon\ll 1$. Consider the approximations to $|\Omega_{1}|^{2}$ and $|\Omega_{3}|^{2}$ up to $\mathcal{O}(\epsilon)$ given by (5.8),

[TABLE]

There exists positive times $\tau_{+}$ and $\tau_{-}$, defined by (5.9), for which, when $\epsilon_{0}=1$, $\lim\limits_{\tau\rightarrow\tau_{+}}\frac{\bar{X}(\tau)}{\bar{Y}(\tau)}=\infty$, indicating evolution to an x-bar state, and, when $\epsilon_{0}=-1$, $\lim\limits_{\tau\rightarrow\tau_{-}}\frac{\bar{X}(\tau)}{\bar{Y}(\tau)}=0$, indicating evolution to a y-bar state. The critical times $\tau_{+}$ and $\tau_{-}$ are $\mathcal{O}(1/\epsilon)$ as $\epsilon\to 0$.

Proof.

Using the expressions given by (5.8), consider the ratio $\frac{\bar{X}}{\bar{Y}}$ as it is a measure of how close the system is to one bar state or another.

[TABLE]

Observe that, if $\epsilon_{0}=1$ the denominator of the ratio in (5.10) decreases monotonically and vanishes at $\tau_{+}=\frac{1}{\nu_{0}K_{x_{0},y_{0}}}\left(\frac{Y_{0}(0)}{\epsilon}+Y_{1}(0)\right)$, and if $\epsilon_{0}=-1$ then the numerator decreases monotonically and vanishes at $\tau_{-}=\frac{1}{\nu_{0}(K_{x_{0},y_{0}}+4X_{0}(0))}\left(\frac{X_{0}(0)}{\epsilon}+X_{1}(0)\right)$. Hence for $\epsilon_{0}=1$, we have $\lim\limits_{\tau\rightarrow\tau_{+}}\frac{\bar{X}(\tau)}{\bar{Y}(\tau)}=\infty$ (indicating an x-bar state) while for $\epsilon_{0}=-1$, we have $\lim\limits_{\tau\rightarrow\tau_{-}}\frac{\bar{X}(\tau)}{\bar{Y}(\tau)}=0$ (indicating a y-bar state). ∎

This proposition illustrates that as time $\tau$ increases in the finite interval on which the approximations in (5.8) are valid, the bar state toward which $\frac{\bar{X}}{\bar{Y}}$ tends is determined solely by $\epsilon_{0}$. The fact that $\tau^{\pm}=\mathcal{O}(1/\epsilon)$ implies that, for larger $\epsilon$, there will be more rapid convergence to the bar state. This could be related to the results of [3], which suggest that $|\delta-1|=\epsilon$ needs to be sufficiently large before a bar state will dominate the dynamics under appropriate stochastic perturbations.

6 Acknowledgements

The authors wish to thank the anonymous referees for their helpful comments. M.B was partially supported by National Science Foundation (NSF) grant DMS 1411460, and K.S. was partially supported by NSF DMS 1550918.

7 Appendix

*Completion of proof of Theorem 3.4

*We have already showed the existence of a family of local invariant manifolds for (3.4). Here we compute the form that this family takes. We begin by splitting (3.4) into real and imaginary parts to get the following 8 dimensional real system of ODEs

[TABLE]

Observe that the above system has a line of equilibrium points $\vec{r}=(r,0,0,0,0,0,0,0)$. The Jacobian of (7.1) at each of these equilibrium points has 7 negative eigenvalues and 1 zero eigenvalue. Thus for each point on the line $\vec{r}$, there is an associated 7 local dimensional stable manifold, $W^{s}(\vec{r})$, and 1 dimensional center manifold, $W^{c}(\vec{r})$, which is in fact globally defined and corresponds to the line $\vec{r}$ itself. Below we will proceed to explicitly compute this family of stable manifolds, denoted $\mathcal{M}_{r}$, up to and including the quadradic terms. This is a long and tedious calculation but comes from standard invariant manifold theorems from ODE, see Theorem 4.1 in [4]. Since the linearization of (7.1) has 7 negative eigenvalues and 1 zero eigenvalue, this theorem guarantees the existence of an 8 dimensional center-stable manifold associated to each fixed point. As we can identify the line of equilibrium points, $\vec{r}$, as the globally defined center manifold, we indeed know that the following computation will result in the codemension 1 stable manifold, i.e. the center-stable manifold is the union of $M_{r}$ and $\vec{r}$. We begin the derivation by first shifting coordinates to move the equilibrium point at $(r,0,0,0,0,0,0,0)$ to the origin.

[TABLE]

The resulting system with fixed point the fixed point at the origin corresponding to $\vec{r}$ is

[TABLE]

The Jacobian of (7.2) at the origin is

[TABLE]

Let us define

[TABLE]

[TABLE]

Then we can write the system as

[TABLE]

Before computing the stable manifold we must change variables and diagonalize the matrix $J$. The matrix of eigenvectors of $J$ and its inverse are given by

[TABLE]

So we have

[TABLE]

Defining

[TABLE]

the dynamics of $\mathbb{Y}$ are given by

[TABLE]

So now we have

[TABLE]

For notational convenience, let $\lambda=y_{1}+\frac{r+1}{4\nu}y_{3}-\frac{r+1}{4\nu}y_{4}-\frac{r+1}{2\nu}y_{5}+\frac{r+1}{2\nu}y_{7}$. Then, we have

$N(S\mathbb{Y})$=

[TABLE]

After computing $S^{-1}N(S\mathbb{Y})$ we see that the complete system in the $\mathbb{Y}$ variables is

[TABLE]

Note that in grouping $\Theta(2)$ terms in the $\dot{y}_{1}$ equation in (7.3) we have:

[TABLE]

Now we are ready to match coefficients on the stable manifold. We know that the stable manifold takes the form $y_{1}=h(y_{2},y_{3},y_{4},y_{5},y_{6},y_{7},y_{8})$ and goes through the origin with parabolic tangency. So let:

[TABLE]

Its derivative with respect to $t$ is:

[TABLE]

Matching coefficients of among the quadratic terms in $\dot{h}$ and $\dot{y}_{1}$, we obtain the nonzero coefficients below.

[TABLE]

[TABLE]

[TABLE]

The stable manifold for the origin in the $\mathbb{Y}$ variables becomes:

[TABLE]

In the original variables, the stable manifold can be written as:

[TABLE]

or, equivalently,

[TABLE]

concluding the proof of the theorem.