Non-equilibrium Ensembles for the three-dimensional Navier-Stokes equations

Here we consider the classical case of an incompressible fluid enclosed in a three-dimensional container with periodic boundary conditions, described by the NSE. Assuming incompressibility amounts to removing internal energy from the energy content of the fluid.

The NSE can be written as

where ${\bm{u}}({\bm{x}},t)$ is imposed to have zero divergence and zero spatial average, $p$ is the pressure field, $\nu$ is the viscosity, and ${\bm{f}}$ is a force field. The NSE is fundamentally difficult in 3D because it is not known whether solutions could be constructed with sufficient generality. For instance, if we suppose that ${{\bm{f}}}$ acts only at large scales (i.e., confined to a finite number of modes) and the initial data have only a finite number of modes, then it is not known whether, in such a generality and no matter how small $\nu>0$ is fixed, a smooth solution follows for all $t>0$. Note that this formulation of the problem is slightly weaker than the millennium problem B of the Clay Mathematics Institute (29).

On the other hand, equally hard problems arise in statistical mechanics systems; for instance, in 3D there is no existence theorem for an infinite system of hard spheres starting at an initial state in which the maximal speed and the minimal particles distance in any unit box are bounded away from $\infty$ and $0$. Notwithstanding, this has not been an obstacle for developing the statistical mechanical theory of phase transitions.

The obvious way forward is to introduce extra parameters which regularize the equations, turning them into equations which admit solutions and try to study only properties that can be shown to be independent of the regularization parameters. For instance, in the case of statistical mechanics the equations are typically regularized by confining the system to a finite box, say a cube of volume $V$.

In the present case we consider the truncated NSE, i.e., the regularized version of Eq. (1) in Fourier space obtained by requiring that ${\bm{u}}$ is periodic in the container and has only modes ${{\bm{k}}}=(k_{1},k_{2},k_{3})$ with $k\leq N$, and $k{\buildrel def\over{=}}||{{\bm{k}}}||_{2}=\sqrt{k_{1}^{2}+k_{2}^{2}+k_{3}^{2}}$. We focus on properties of the solutions which hold uniformly in the cut-off $N$, and the container is supposed to be of size $[0,2\pi]^{3}$. Therefore, one can introduce the complex scalars $u_{\beta,{{\bm{k}}}}={\overline{u}}_{\beta,-{{\bm{k}}}}$, where the overline notation denotes the complex conjugate, and two unit vectors ${\bf e}_{\beta}({{\bm{k}}})=-{\bf e}_{\beta}(-{{\bm{k}}})$ mutually orthogonal and orthogonal to ${{\bm{k}}}$, so that the velocity field can be written as

where the sum is restricted to $k\leq N$. Furthermore, by defining the kernel $T_{{{\bm{p}}},{{\bm{q}}},{{\bm{k}}}}^{\beta_{1},\beta_{2},\beta}=-({\bm{e}}_{\beta_{1}}({{\bm{p}}})\cdot{{\bm{q}}})({\bm{e}}_{\beta_{2}}({{\bm{q}}})\cdot{\bm{e}}_{\beta}({{\bm{k}}})),$ and making the time notation implicit unless strictly necessary, the NSE can be expressed as

with $k,p,q\leq N$. Due to the symmetry between ${{\bm{q}}},{{\bm{k}}}$, the sum over $\beta_{1},\beta_{2}$ keeping ${{\bm{p}}}+{{\bm{q}}}+{{\bm{k}}}={\bm{0}}$ yields the identity

reflecting the conservation of both the total energy, $\int d{{\bm{x}}}\,{\bm{u}}^{2}$, and the total helicity, $\int d{{\bm{x}}}\,{\bm{u}}\cdot({\mbox{\boldmath$\partial$}}\wedge{\bm{u}})$ in the unforced and inviscid ($\nu=0,{{\bm{f}}}={\bm{0}}$) limits.

Historically, discussions on the origin of dissipative macroscopic properties, out of purely reversible Newtonian dynamics at molecular level, go back at least to Maxwell [see Eq. (128) in (30)]. Therefore, it is natural to expect that fluid motion can also be described by microscopic reversible equations. One can thus inquire whether even at the macroscopic level, although molecular motion is no longer explicitly playing a role, the same phenomena can be described by macroscopic reversible equations.

Here we shall study stationary states of the truncated NSE. The basic idea is that viscosity controls the regularity of the flow by forbidding uncontrolled growth of energy, and dissipates the input work by the forcing: $\int{{\bm{f}}}\cdot{\bm{u}}\,d{{\bm{x}}}$, proportionally to the product of viscosity times enstrophy

and where at steady state there is a statistical balance between input and output of energy. One can envisage that by replacing the viscous force $\nu\Delta{\bm{u}}$ with a reversible force that keeps the enstrophy constant many statistical properties of stationary states will remain correctly described, independently of the cut-off $N$ (31).

The analogy with statistical mechanics is helpful; there, the cut-off is the volume of the container, and the equilibrium states can be described by different microscopic dynamics, like the energy conserving microcanonical distribution or the isokinetic evolution (14). Both choices of states and evolution laws give the right statistics for local observables, i.e., observables depending on the configurations of particles located in a volume small compared to the total one, i.e., to the cut-off $V$. Most importantly, the statistical properties of such observables depend on $V$ less and less as $V\to\infty$.

In the case of NS fluids subject to large scale forcing, locality is defined with respect to the Fourier space. The analogs of local observables are functions of the velocity field ${\bm{u}}$ that only depend on components ${\bm{u}}_{{\bm{k}}}$ with $k$ small compared to the cut-off $N$. More specifically, $O({\bm{u}})$ is a local observable, when it depends only on a finite number of modes with $k<K$ for some $K<N$. This leads to the proposal that replacing the standard viscous force with a dissipative reversible term that keeps the enstrophy constant will lead to stationary states completely equivalent to the original dynamics up to an arbitrary $K$, as far as the statistical properties of local observables are concerned, provided $N$ is large (ideally $N\to\infty$, in the same sense in which equilibrium thermodynamics becomes ensemble independent as the volume tends to infinity). Here $K$ can in principle depend on the viscosity and it is important to understand its scaling properties in the $\nu\to 0$ limit (see below).

Specifically, we consider a different evolution equation where (i) in Eq. (3) the external force ${{\bm{f}}}$ is fixed, with $||{{\bm{f}}}||_{2}=\text{const}\in\mathbb{R}$, and acts on “large scale”, i.e., it has finitely many modes: ${{\bm{f}}}_{{\bm{k}}}=0$ if $k\geq k_{f}$ for some $k_{f}$; (ii) the viscous force $-\nu k^{2}{\bm{u}}_{{\bm{k}}}$ is replaced by $-\alpha({\bm{u}})k^{2}{\bm{u}}_{{\bm{k}}}$, with $\alpha({\bm{u}})$ defined such that the enstrophy

is a constant of motion. Note that a similar choice was followed by the authors in (28), while in (27) the kinetic energy was instead kept fixed. Such assumptions lead to the equations

where

with

obtained by multiplying both sides of Eq. (6) by $k^{2}{\overline{u}}_{\beta,{{\bm{k}}}}$, then summing over ${{\bm{k}}}$, and finally imposing that the right-hand side equals 0.

We call Eq. (3) irreversible NSE (abbreviated as ${\cal I}$) and Eq. (6) reversible NSE (abbreviated as ${\cal R}$). The name reversible refers to the property that if $S_{t}{\bm{u}}={\bm{u}}(t)$ is a solution of ${\cal R}$, with initial data ${\bm{u}}_{0}$, and if ${\cal P}{\bm{u}}=-{\bm{u}}$, then ${\cal P}S_{t}{\bm{u}}=S_{-t}{\cal P}{\bm{u}}$ is also a solution, as a consequence of $\alpha({\cal P}{\bm{u}})=-\alpha({\bm{u}})$. Instead, such an identity does not hold for the solutions of ${\cal I}$.

Let us introduce the collection ${\mathcal{E}}^{{\cal I},N}$ of the stationary distributions $\mu^{{\cal I},N}_{\nu}$ for the ${\cal I}$  evolutions, with cut-off $N$; for a given choice of ${{\bm{f}}}$ the distributions are parametrized by the Reynolds number ${\rm Re}\propto\nu^{-1}$. In principle, several distributions may correspond to the same ${\rm Re}$ although it is plausible that for ${\rm Re}$ large enough there is only one stationary distribution. Likewise, let ${\mathcal{E}}^{{\cal R},N}$ be the collections of the stationary distributions $\mu^{{\cal R},N}_{D}$ for the ${\cal R}$ evolutions, with the same cut-off $N$. These are parametrized by the value $D$ of the (constant) enstrophy. Again, in general, several distributions may correspond to the same $D$; similarly to the irreversible case, it is plausible that for generic initial data and $D$ large, $\mu_{D}^{R,N}$ is unique.

It is natural to associate to each distribution in ${\mathcal{E}}^{{\cal I},N}$ and ${\mathcal{E}}^{{\cal R},N}$ the average enstrophy $\langle{\cal D}({\bm{u}})\rangle_{\nu}^{{\cal I},N}$ and, respectively, the enstrophy $D$, as well as the work per unit time dissipated:

where $\langle O\rangle^{{\cal I},N}_{\nu}$, $\langle O\rangle^{{\cal R},N}_{D}$ denote the averages of an observable $O$ over the distributions $\mu^{{\cal I},N}_{\nu}$ and $\mu^{{\cal R},N}_{D}$. Parameters $\nu,D$ at given $N$, will be said to be in correspondence if

Equation (10) defines an implicit relationship between $D$ in ${\cal R}$ and $\nu$ in ${\cal I}$.

In this work, we perform a few preliminary checks, based on a numerical analysis of the truncated NSE, of the following two conjectures on the distributions of local observables, originally presented in (17; 18; 1), respectively.

Let us introduce the Kolmogorov scale $k_{\nu}=\left(\frac{\varepsilon}{\nu^{3}}\right)^{1/4}$, with $\varepsilon=\nu\langle{\cal D}({\bm{u}})\rangle$ being the energy dissipation rate, defined as the typical length where inertial and irreversible dissipative effects balance. We can then state a different conjecture in the limit of $N\to\infty$.

Remarks:

1. 1.

   Conjecture 1 relies on the chaoticity of the evolution at small $\nu$ and fixed $N$. Instead, conjecture 2 concerns only the limit $N\to\infty$ and relies on the chaoticity of the microscopic motions leading to the NSE. Conjecture 2 is formulated for all $\nu$, including the case where the asymptotic motion consists of periodic attracting sets. It just provides a quantitative version of conjecture 1: given a local observable living on scale $K$ the question of how small should $\nu$ be to achieve equivalence receives the quantitative answer that, in the limit $N\to\infty$, $\nu$ has to be such that $K<c_{0}k_{\nu}$ which becomes eventually true, as lower and lower values of $\nu$ are considered, because $k_{\nu}$ diverges as $\nu\to 0$.

2. 2.

   Conjecture 2 could be extended to more general macroscopic equations derived via scaling limits from microscopic reversible evolution.

3. 3.

   In cases with more than one distribution for a given $\nu$ or $D$ the interpretation should be that extra labels should be added to distinguish the various distributions and follow the remark in conjecture 1. Conjecture 2 first appeared in (24) and was formally proposed in (25).

4. 4.

   In (2; 26) a stronger conjecture is directly formulated with $c_{\nu}=\infty$ for all $\nu$.

The two conjectures differ in the order of consideration of the two limits, $\nu\to 0$ and $N\to\infty$. In particular, it is important to stress that the so-called fully developed turbulent limit is achieved in nature by $N\to\infty$ first and then $\nu\to 0$, i.e., it is captured by conjecture 2. On the other hand, the limit $\nu\to 0$ for fixed $N$ leads to a quasi-thermalization (in the sense of (32)) of the high wave numbers range cut-offed by the maximum $N$. An additional feature of this regime is that the dimension of the attractor approaches the number of degrees of freedom of the system (20).

We have performed direct numerical simulations of both ${\cal I}$ and ${\cal R}$ equations for incompressible fluid in a triply periodic domain of size $L=2\pi$. We used a dealiased (following the $\frac{2}{3}$ rule (33)) parallel 3D pseudospectral code, based on the P3DFFT implementation (34), on cubic grids of size $N_{0}=64$, $128$, $256$, and $512$ collocation points in each direction, effectively corresponding, via $N=N_{0}/3$, to $N=21$, $42$, $85$, and $170$ in Fourier space, respectively. $N$ is the same as the cut-off scale $k_{max}$ and will be used interchangeably. The time integration has been implemented with a second-order Adams-Bashforth scheme and a small timestep equal to $\Delta t=10^{-4}$ was chosen for all the simulations; tests with different time steps have been done confirming that the one chosen was sufficient to ensure the robustness of the presented results. The grid size is $\Delta x=2\pi/N_{0}$. The zero mode is not forced, ensuring that the mean velocity field remains zero at all times, and the external forcing is deterministic and only acts on large scales, specifically in the $k_{f}\in[1,\sqrt{2}]$ domain. The phases were chosen as quenched random variables and kept equal in all simulations with $\sum_{{\bm{k}}}k^{2}||{\bm{f}}||_{2}^{2}=1$. We refer to Table 1 in Appendix C for a summary of the numerical simulations we performed, which are labeled as R#, where R stands for run in this context, followed by the number of the run, and a symbol, indicating the statistical regime each run belongs to; see discussion in Sec. III.2.

We examine conjectures 1 and 2 by testing a set of viscosities: $\nu=5\times 10^{-2},10^{-2},10^{-3},10^{-4},10^{-5}$ for different values of $N$. We proceed as follows. First we run a reasonably long ${\cal I}$ simulation at fixed $\nu$, $N$ and measure the average enstrophy $\langle D\rangle^{{\cal I},N}_{\nu}$ in the statistically steady state. Next, we start the ${\cal R}$ run from a steady ${\cal I}$ state for which the instantaneous enstrophy is $D=\langle D\rangle^{{\cal I},N}_{\nu}$. This step reduces the transient time required to reach the attracting set, but it is not strictly necessary, and in principle one can start from any initial condition with the appropriate enstrophy. We run the ${\cal R}$ for the same duration as the ${\cal I}$, and at each timestep we make tiny corrections by rescaling the velocity field as ${\bm{u}}\rightarrow{\bm{u}}\sqrt{\frac{\langle{D}\rangle^{{\cal I},N}_{\nu}}{D(t)}}$. Although, given the very small timestep we use, the deviations from condition (10) were tested to be small, the latter correction ensures that the total enstrophy stays constant within machine precision.

Two distinct statistically steady regimes can be identified for the ${\cal I}$ and ${\cal R}$ systems, depending on the cut-off $N$ and Re (or, equivalently, on $k_{\nu}$). One is what we call the hydrodynamic regime, achievable when $N$ is large enough for any $\nu$, such that $k_{\nu}\lesssim N$. It is characterized by a well developed inertial range (for small $\nu$) with a $E(k){\sim}k^{-5/3}$ scaling for $k_{f}\ll k\ll k_{\nu}$ and a well resolved exponential or superexponential decay for large wave numbers, where

are the averaged spectral properties, or simply the energy spectrum. A scaling close to $E(k){\sim}k^{-5/3}$ is observed in nature because of the regularizing properties of viscosity at small scales (3; 4; 6). The other, -quasithermalized- regime (32), features a $E(k){\sim}k^{2}$ scaling for large values of $k$ (32) and it is obtained at sufficiently small $\nu$ for any fixed $N$, such that $k_{\nu}\gg N$. We call this regime quasithermalized because the average energy content of the individual modes depends only weakly on $k$, given the geometrical degeneracy of the energy shells. A crossover region can be located between the two in the ($N,\nu$) parameter space (27).

In Fig. 1 we present a first comparison between ${\cal I}$ and ${\cal R}$ by looking at $E(k)$ at changing $\nu$ and for fixed spatial resolution, $N_{0}=128$. Here, we have omitted in $\langle\cdot\rangle$ the label ${\cal I}$ and ${\cal R}$ for simplicity and the same will be done in what follows whenever it does not lead to ambiguities. All regimes can be accessed as the viscosity is varied, and in terms of the average energy spectrum, the ${\cal I}$ and ${\cal R}$ agree remarkably well. This will be further checked in the following with respect to the two conjectures. We anticipate that deviations will appear at scales $k$ beyond the Kolmogorov scale $k_{\nu}$, which will be further explored later. Next we will deal with the two regimes separately.

A rigorous (yet non trivial, see below) consequence of both conjectures is the relation

which holds when Eq. (10) holds. In Fig. 2(a) we validate Eq. (14). The errors are calculated by $\delta\mathcal{O}=\sigma(\mathcal{O})\sqrt{1+2\tau_{\mathcal{O}}}/\sqrt{j},$ where $\sigma(\mathcal{O})$ is the standard deviation, $\tau_{\mathcal{O}}$ the non-dimensional autocorrelation time, and $j$ the ensemble size.

Figure 2(a) provides a consistency check, which looks at first not entirely obvious because $\alpha$ is not a local observable; see the definition given in Eq. (7). Yet, if the conjectures and condition (10) hold, it must be true that the ensemble average of $\alpha$ is equal to the viscosity. This follows from the balance between energy input and output in steady state conditions that hold both for ${\cal I}$ and ${\cal R}$:

Under the equivalence condition $\langle{\cal D}({\bm{u}})\rangle^{{\cal I},N}_{\nu}=D$, the value $\lim_{N\to\infty}\langle{{\bm{f}}}\cdot{\bm{u}}\rangle$ must be the same in ${\cal I}$ and ${\cal R}$, because ${{\bm{f}}}\cdot{\bm{u}}$ is a local observable. Hence, the relation $\langle{\cal D}({\bm{u}})\rangle^{{\cal I},N}_{\nu}=D$ imposes the equality of the averages $\lim_{N\to\infty}\langle{{\bm{f}}}\cdot{\bm{u}}\rangle$ for the ${\cal I}$ and ${\cal R}$ systems, which, in turn, implies the non trivial relation (14). We remark that, as $\alpha({\bm{u}})$ can be measured also in ${\cal I}$, we can also check its statistics. We observe that $\langle\alpha({\bm{u}})\rangle^{{\cal I},N}_{\nu}\approx\nu$ in all statistical regimes, which is nontrivial because it is not a consequence of the equivalence conjecture because $\alpha$ is not a local observable. Although the distributions of $\alpha({\bm{u}})$ in both ${\cal R}$ and ${\cal I}$ are in all cases peaked around $\nu$, as well as the mean values of both $\langle\alpha({\bm{u}})\rangle^{{\cal R},N}_{D}$ and $\langle\alpha({\bm{u}})\rangle^{{\cal I},N}_{\nu}$ are approximately equal to $\nu$, their tails are different in the hydrodynamic regime, but become almost identical in the quasithermalized regime (see Sec. IV.2).

Since $\alpha$ is a sum of the forcing contribution $W({\bm{u}})/\Gamma({\bm{u}})$ and of the internal nonlinear contribution $\Lambda({\bm{u}})/\Gamma({\bm{u}})$ (8), it is interesting to remark that, on average, the non linear term dominates the forcing one as shown in Fig. 2(b), where the ratio $\langle\Lambda\rangle/\langle W\rangle$ versus $\nu$ is plotted. The ratio increases linearly as $\nu$ decreases in the hydrodynamic regime, before it approaches a (large) constant depending on $N$ in the quasithermalized regime. Thus, on average the internal nonlinear exchanges dominate over the forcing for the ${\cal R}$ dynamics, introducing a non-local (in scale) coupling among all modes and making the validity of conjecture 2 even less trivial.

In this section we address the validity of conjecture 1 for 3D NSE. We remark that it has been confirmed in simpler dynamical systems in (18; 20; 21; 22). conjecture 1 applies when we consider decreasing values of $\nu$ at fixed $N$, which leads to the quasithermalized regime. In terms of observables, we consider the single mode kinetic energy, $U({{\bm{k}}})=||{\bm{u}}(k_{1},k_{2},k_{3})||_{2}^{2}$ and the energy spectrum, $E(k)$ (13) which will be presented separately.

We now present our numerical tests of conjecture 2, which applies at fixed $\nu$ when $N$ is large, corresponding to the hydrodynamic regime. We follow a systematic approach by first keeping the value of $\nu$ fixed, and then progressively increase the value of $N$. We then repeat this protocol for different values of $\nu$. Again, we consider the single mode kinetic energy and the energy spectrum.

In this section we inspect the statistics of the fluctuations of $\alpha({\bm{u}})$ for ${\cal R}$, in order to assess the probability to observe negative values at varying the control parameters. In Fig. 15(a) we show in light colors the temporal evolution of $\alpha/\nu$, and in solid colors its running average, performed as $\frac{1}{t-t_{0}}\sum_{t^{\prime}=t_{0}}^{t}\alpha({\bm{u}}(t^{\prime}))/\nu$ for ${\cal R}$ at $\nu=10^{-3}$ and $N_{0}=64,\,512$ (green - R3$\triangle$ and red - R16$\bigcirc$, respectively). The presented time window, which is scaled by the Kolmogorov time $\tau_{K}=\sqrt{\nu/\varepsilon}$, is the maximum achieved for the run with $N_{0}=512$ (R16$\bigcirc$), while for R3$\triangle$ this is only about 1% of its total temporal extend. Then, in Fig. 15(b) the PDFs of $\alpha/\nu$ are presented for $N_{0}=64$, 128, 256, and 512 at $\nu=10^{-3}$. Notably, all PDFs agree, even though at $\nu=10^{-3}$ only $N_{0}=512$ is fully resolved, belonging therefore to the hydrodynamic regime, while the other simulations belong to the crossover regime. The inset plots (i) and (ii) are respectively the energy spectrum and the PDF in semilogarithmic scale.

In Fig. 15(c) we show the same of Fig. 15(b) but at lower viscosity, $\nu=10^{-4}$, where the data set at $N_{0}=64$ (R4$\square$) is in the quasithermalized regime, while the other runs are in the crossover regime. As one can see, except for the case of the quasithermalized regime, in all other simulations we do not observe negative values of $\alpha$ within our statistical sample.

Additionally, in Fig. 15(d) we perform a detailed study at fixed $\nu=10^{-3}$ and low resolution, in order to observe when, and how the transition to the occurrence of frequent negative values of $\alpha$ takes place. Starting from $N_{0}=64$ and gradually decreasing $N_{0}$, we observe that the PDF of $\alpha$ first turn from non-Gaussian to quasi-Gaussian, while becoming narrower (see $N_{0}=48$ and $N_{0}=44$), before starting to widen (see $N_{0}=40$), until occurrence of $\alpha<0$ events appears (see $N_{0}\leq 36$). The rapid change in the left tail of the PDF at changing $N_{0}$ suggests that in order to further substantiate the asymptotic probability to observe negative $\alpha$ values, when $N_{0}$ is large and in the presence of intermediate or hydrodynamical regimes, one would need statistical samples significantly (exponentially) larger than those we could generate here, which is by far out of the scope of this paper and probably out of reach given the available computational power. This result is in agreement with some previous observations made in the context of a reduced model of turbulence (21).

In (28) negative values of $\alpha$ are observed in an reversible ensemble at $N_{0}=1024$ and Taylor-based Reynolds number Re${}_{\lambda}=300$. Indeed, such results are somewhat puzzling and are not expected on the basis of our data. We do not expect to be able to observe $\alpha<0$ at such large $N_{0}$ and such Reynolds number by extrapolating from our results for the fully resolved hydrodynamical regime (see Sec. IV.2 below).

Interestingly enough, $\alpha({\bm{u}})$ is an observable that can be measured also in the ${\cal I}$ ensemble, by computing Eq. (7). In Fig. 16 we present the cases for two viscosities, $\nu=10^{-2}$ [(a) and (b)] and $\nu=10^{-3}$ [(c) and (d)]. In Figs. 16(a), 16(c) we show the time series of $\alpha({\bm{u}})/\nu$ for a time window scaled by the Kolmogorov time scale $\tau_{K}$, and in Figs. 16(b), 16(d) the corresponding PDF for different different resolutions $N_{0}$.

First, from Figs. 16(a)–16(c) we note that also in ${\cal I}$ the non-trivial result $\langle\alpha({\bm{u}})\rangle^{{\cal I},N}_{\nu}=\nu$ holds (see also Sec. III.3). On the other hand, from panels Figs. 16(b)–16(d) we observe some different behavior from the corresponding PDF of the ${\cal R}$ ensemble in particular for the case of $\nu=10^{-3}$. In the latter case the high wave number range starts to depart from the hydrodynamic regimes and probably has a direct impact on the fluctuations of $\alpha$.