Exact coherent structures in fully developed two-dimensional turbulence

Abstract

This paper reports several new classes of unstable recurrent solutions of the two-dimensional Euler equation on a square domain with periodic boundary conditions. These solutions are in many ways analogous to recurrent solutions of the Navier–Stokes equation which are often referred to as exact coherent structures. In particular, we find that recurrent solutions of the Euler equation are dynamically relevant: they faithfully reproduce large-scale flows in simulations of turbulence at very high Reynolds numbers. On the other hand, these solutions have a number of properties which distinguish them from their Navier–Stokes counterparts. First of all, recurrent solutions of the Euler equation come in infinite-dimensional continuous families. Second, solutions of different types are connected, e.g. an equilibrium can be smoothly continued to a travelling wave or a time-periodic state. Third, and most important, they are only weakly unstable and, as a result, fully developed turbulence mimics some of these solutions remarkably frequently and over unexpectedly long temporal intervals.

1. Introduction

The important role of coherent structures in fluid turbulence has been both widely accepted and hotly debated for several decades. However, many open questions still remain. A quote from the seminal paper on this topic by Hussain (1986) sets the stage for the present study: ‘The interaction between coherent structures and incoherent turbulence is the most critical and least understood aspect of turbulent shear flows. This coupling appears to be rather different from the classical notion of cascade; even considering the large and fine scales, they are not decoupled as widely presumed. The coupling can be intricate and of different kinds…’. Indeed, this pretty much sums up the state of knowledge to the present day.

We will follow Hussain in distinguishing coherent structures as those of size comparable to the transverse length scale of the shear flow as opposed to coherent substructures whose characteristic size corresponds to the Taylor microscale. These scales are only clearly distinguished at high Reynolds numbers ($Re$) where coherent structures can be considered inviscid. Rather interestingly, the greatest progress in understanding coherent structures to date has been made in the context of weakly turbulent flows where this scale separation disappears. Advanced numerical methods such as Newton–Krylov solvers (Viswanath 2007) enabled computation of unstable, recurrent (e.g. steady, time-periodic or relative-periodic) solutions of the governing equations in a variety of canonical shear flows. Many of such numerical solutions were found to have spatiotemporal structure similar to that of familiar coherent structures such as streamwise vortices and velocity streaks and hairpin vortices near a wall (Waleffe 1998, 2001; Gibson, Halcrow & Cvitanović 2008; Itano & Generalis 2009; Shekar & Graham 2018); consequently, they have been termed exact coherent structures (ECSs). Recent numerical and experimental studies demonstrated that ECSs do not merely resemble turbulent flows, they also organize and guide the dynamics of weak turbulence in both two (Suri et al. 2018, 2020) and three spatial dimensions (Krygier, Pughe-Sanford & Grigoriev 2021; Crowley et al. 2022).

Exact coherent structures have already generated significant insight into weakly turbulent flows (Kawahara, Uhlmann & Van Veen 2012; Graham & Floryan 2021). Most notably, ECSs are found to capture self-sustaining processes that maintain wall-bounded turbulence (Waleffe 1997; Hall & Smith 1991). Exact coherent structures also elucidate the transition from laminar flow to turbulence explaining both the ‘bypass’ mechanism in linearly stable flows (Khapko et al. 2016) and the formation of chaotic sets underpinning sustained turbulence (Kreilos & Eckhardt 2012). Despite these successes, due to the lack of scale separation in weakly turbulent flows, it remains unclear how much of this understanding carries over to fully developed turbulence.

Extending the ECS framework to higher $Re$ proved challenging due to both conceptual and technical challenges. As $Re$ increases, the range of scales accessible to turbulence becomes larger, and the number of distinct ECSs grows very quickly. Furthermore, each ECS becomes more unstable; and, to complicate things even further, it becomes even more expensive to find ECSs through direct numerical simulations (DNS) which effectively become intractable at high $Re$ due to the high spatial and temporal resolution requirements. A key conceptual challenge is related to the limit $Re\to \infty$. The Euler equation is a singular limit of the Navier–Stokes equation, which makes it difficult to establish a relation between dynamically relevant recurrent solutions of the two equations beyond some relatively loose, although important, constraints (Batchelor 1956; Okamoto 1994; Gallet & Young 2013).

One of the most notable successes in continuing recurrent solutions to higher $Re$ is the computation of attached eddies (Eckhardt & Zammert 2018; Doohan, Willis & Hwang 2019; Yang, Willis & Hwang 2019; Azimi & Schneider 2020) and their bulk analogues (Deguchi 2015; Eckhardt & Zammert 2018) in a variety of canonical three-dimensional (3-D) shear flows. In wall units, these solutions become independent of $Re$, just like the wall-bounded fully developed turbulent flows. It is important to note that these solutions represent coherent substructures. To the best of our knowledge, no examples of ECS have been found in high-$Re$ 3-D turbulent flows at the scale comparable to either the system size (e.g. distance between the boundaries) or the scale of the forcing; such ECSs correspond to recurrent solutions of the Euler equation. It should be pointed out, however, that, some – not necessarily dynamically relevant – steady and time-periodic solutions representing large-scale flows at fairly high $Re$ have been computed through continuation in two spatial dimensions (Kim & Okamoto 2010, 2015; Kim, Miyaji & Okamoto 2017) and three spatial dimensions (Wang, Gibson & Waleffe 2007). In principle, recurrent large-scale flows could also be computed with the help of large eddy simulations (Rawat et al. 2015; Hwang, Willis & Cossu 2016), although such flows do not represent formal solutions of either the Navier–Stokes or Euler equation.

The objective of this paper is therefore to introduce such large-scale recurrent solutions of the incompressible Euler equation in two spatial dimensions, where fully resolved DNS of turbulent flow can be performed for relatively high $Re$. There is a rich history of computing solutions of the Euler equation using various analytical tools. Even leaving singular solutions involving point vortices aside, quite a few examples of absolute and relative equilibria have been found. They include isolated circular and elliptic vortices (Lamb 1924), translating pairs of counter-rotating vortices (Pierrehumbert 1980; Saffman & Tanveer 1982), rotating arrays of two or more vortices (Dritschel 1985; Carton, Flierl & Polvani 1989; Crowdy 2002b) and stationary multipolar vortex arrays (Crowdy 1999, 2002a; Tur & Yanovsky 2004; Xue, Johnson & McDonald 2017).

This paper reports several new classes of recurrent solutions – equilibria, travelling waves and time-periodic states – of the two-dimensional (2-D) Euler equation on a domain with periodic boundary conditions which can be thought of as generalizations of vortex crystals (Aref et al. 2003). A distinguishing feature of these solutions is their dynamical relevance for high-$Re$ statistically stationary turbulent flows which balance forcing and dissipation: just like in the case of weakly turbulent flows, we find recurrent solutions to organize and guide the dynamics of turbulent flow on large scales. As expected, the library of inviscid recurrent solutions is found to be far larger than that describing viscous flows: unlike their analogues for weakly turbulent flows which are isolated, solutions of the Euler equation are found to belong to continuous families spanned by an infinite number of parameters.

The rest of the paper is organized as follows. Section 2 describes the problem set-up investigated here, and § 3 discusses the relation between solutions of the Euler and Navier–Stokes equation in the high-$Re$ limit. Our results are presented in § 4, their implications for the problem of fully developed fluid turbulence are discussed in § 5, and conclusions are presented in § 6.

2. Problem description

We consider an incompressible Newtonian fluid in two spatial dimensions. In the presence of forcing and dissipation, its dynamics are governed by the Navier–Stokes equation,

(2.1)\begin{equation} \partial_t\boldsymbol{u}+ (\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla})\boldsymbol{u} ={-}\boldsymbol{\nabla} p+Re^{{-}1}\nabla^2\boldsymbol{u}+\boldsymbol{f}, \end{equation}

where $p$ is the pressure, $\boldsymbol {f}$ represents the external forcing and $Re$ is the Reynolds number. In two spatial dimensions, it can be conveniently rewritten in terms of vorticity $\omega =\partial _xu_y-\partial _yu_x=-\nabla ^2\psi$,

(2.2)\begin{equation} \partial_t\omega+ (\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla})\omega = Re^{{-}1}\nabla^2\omega+\varphi, \end{equation}

where $u_x=\partial _y\psi$ and $u_y=-\partial _x\psi$ are the components of the velocity field, $\psi$ is the stream function, and $\varphi =\partial _xf_y-\partial _yf_x$. In this study, we confine our attention to spatially periodic domains $0\le x,y < 2{\rm \pi}$ with stationary checkerboard forcing

(2.3)\begin{equation} \varphi=\sin(Nx)\sin(Ny), \end{equation}

where $N=4$. With this choice of non-dimensionalization, the length, time, velocity and vorticity scales are all $O(1)$.

The frequency of the spatial forcing, $k_f=\sqrt {2}N\approx 6$ is chosen to be reasonably well separated from both the dominant wavenumber $k_0=1/\sqrt {2}$ describing the large-scale flow on the low wavenumber side and the inertial range on the higher wavenumber side. This choice of the forcing frequency ensures that the direct (enstrophy) cascade takes place over most of the wavenumber range resolved in the simulations while the inverse (energy) cascade is constrained to the narrow range $[k_0,k_f]$. As far as the structure of the large-scale flows at high $Re$ is concerned, the particular choice of the forcing profile does not appear to play a noticeable role according to both our simulations and the results of previous systematic studies (Gallet & Young 2013; Kim et al. 2017). Note that, unlike, say, a Kolmogorov forcing profile (Tithof et al. 2017), a checkerboard profile breaks continuous translational symmetry in both spatial directions, constraining the types of ECSs that can be found at moderate $Re$ (Suri 2021). This symmetry breaking effectively disappears at higher $Re$, however, as discussed below.

To generate turbulent flows, we computed solutions of (2.2) numerically using a pseudospectral method. We used a variable time step Runge–Kutta–Fehlberg scheme where spatial derivatives and all linear terms were computed in Fourier space and the nonlinear term was computed in physical space. Additionally, we used a two-thirds dealiasing scheme for numerical stability. Most of the results reported here used the grid resolution of $512\times 512$.

In this study we set the Reynolds number to a relatively high value of $Re=10^5$ in order for structure on a broad range of scales to develop, as illustrated by figure 1(a). The energy spectrum is found to exhibit a clear power-law scaling $E(k)\propto k^\alpha$ over at least a decade in the wavenumbers ($16\le k\le 170$). This scaling indicates the presence of an inertial range characteristic of fully developed turbulence and clear separation between the $O(1)$ length scale of the forcing and the Taylor microscale $k_t^{-1}$. The exponent $\alpha \approx -4.5$ of the power law is found to differ substantially from the value of $-$3 predicted by the classical theory of turbulent cascades in 2-D turbulence developed by Kraichnan (1967), Leith (1968) and Batchelor (1969). This discrepancy is well known (Boffetta & Ecke 2012) and demonstrates the limitations of the Kraichnan–Leith–Batchelor (KLB) theory in properly accounting for the effect of coherent structures on the direct cascade, despite its acknowledgement of the non-local nature of interactions (Kraichnan 1967, 1971). Indeed, figure 1(b) shows that, although the energy spectrum still retains the power-law scaling when averaged over a short time interval, the exponent $\alpha$ exhibits substantial fluctuations in time reflecting changes in the large-scale structure of the flow.

Figure 1. The energy of turbulent flow. (a) The energy spectrum averaged over a long time interval ($10^3$ non-dimensional units) in the asymptotic regime exhibits a clear power-law scaling $E(k)\propto k^\alpha$ (shown as dashed line) in the inertial range. (b) The exponent $\alpha$ of the power law computed for energy spectra averaged over a characteristic time scale $T_c$. (c) Variation in the energy over a long time interval. Note that $t=0$ in different plots corresponds to arbitrary and different times.

The general structure of coherent components of turbulent flow at large and small scales can be easily recognized by applying a Fourier filter. For illustration, we use here a low-pass Fourier filter, denoted by the linear operator $\hat {L}_k$, which corresponds to a smoothed circular window of radius $k$ in Fourier space. We arbitrarily consider wavenumbers $k\le 16$ as representing large scales and $k\ge 24$ as representing small scales. For reference, taking dealiasing into account, our numerical simulations resolve spatial frequencies from the lowest, $k_{min}=k_0$, to the highest, $k_{max}=\lfloor 512/3\rfloor = 170$. The wavenumber associated with the Taylor microscale at which viscous effects become important is somewhat higher, $k_t\sim Re^{1/2}=O(300)$. As figure 1(a) illustrates, we find power-law scaling over the entire wavenumber range corresponding to small scales while, for large scales, significant deviations from a power law are observed.

A typical snapshot of turbulent flow is shown in figure 2(a). Figures 2(b) and 2(c) show, respectively, the large-scale flow $\hat {L}_{16}\omega$ and the small-scale flow $(1-\hat {L}_{24})\omega$. The large-scale component, for sufficiently low-frequency forcing such as the one considered here, features a pair of counter-rotating vortices. This is an expected result of the inverse cascade that accumulates the energy in the largest scales accessible to the flow. Small scales represent filaments of vorticity which mainly occupy the space between the two vortices. Stretching (or thinning) of these vorticity filaments in the hyperbolic regions of the large-scale flow is believed to be a key physical mechanism behind the direct cascade (Chen et al. 2003).

Figure 2. A typical snapshot of turbulent vorticity field $\omega$. The corresponding large-scale flow $\hat {L}_{16}\omega$ (b) and small-scale flow $(1-\hat {L}_{24})\omega$ (c). Panel (d) shows a snapshot of the converged time-periodic ECS. This and all subsequent plots use the same colour bar. For solutions of the Euler equation, such as that shown in panel (d), vorticity scale is arbitrary due to scaling invariance. Hence, the colour bar shown in panel (a) can be used to interpret vorticity fields shown in all subsequent figures.

For higher-frequency forcing, the vortices at the scale closer to that of the forcing become more prominent, and the spectrum has two separate scaling regions, one between the domain scale and the forcing scale controlled by the inverse cascade and another between the forcing scale and the Taylor microscale controlled by the direct cascade. This more complicated situation is outside the scope of the present study.

The spatial resolution used in our simulations, which may be considered low by modern standards, was motivated by the need to evolve the flow over extremely long time scales. As illustrated in figure 1(c), for the high $Re$ considered here, it takes on the order of $10^5$ non-dimensional time units for the flow to come to a statistical equilibrium and for the energy $E=\|\boldsymbol {u}\|^2/2$ to approach its long-term average $\bar {E}\approx 0.44$. Here and below, we use the bar to denote the temporal mean and the $L_2$-norm defined by the spatial mean

(2.4)\begin{equation} \|\boldsymbol{u}\|^2=\frac{1}{4{\rm \pi}^2}\int_0^{2{\rm \pi}}\int_0^{2{\rm \pi}}\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{u}\,{\rm d}x\,{\rm d}y. \end{equation}

The results presented here describe the statistically stationary turbulence that is found after the long transients have died down we initialize the simulation using a small random perturbation about $\boldsymbol {u}=0$. A sample movie of turbulent flow in the asymptotic regime is provided as supplementary movies available at https://doi.org/10.1017/jfm.2023.584.

To fully resolve the flow at the target $Re$, one formally needs to use an $n\times n$ grid with $n\gtrsim 3k_t\approx 10^3$. Hence, simulations on a $2048\times 2048$ grid can be considered fully resolved. While a number of previous studies of 2-D turbulence have used simulations with this resolution (Smith & Yakhot 1993; Bracco et al. 2000), they computed flows over time intervals many orders of magnitude shorter than those considered here. To make sure that our numerical results are representative of fully resolved simulations of turbulence, we compared them with those obtained on a $2048\times 2048$ grid and found the differences to be negligible on the time scale $T_c=10$ characteristic of the recurrent solutions discussed below. We also found the flows, and the recurrence diagrams discussed below, computed on $512\times 512$ and $2048\times 2048$ grids, to remain qualitatively similar on time scales of $O(10^3)$.

3. The $Re\to \infty$ limit

It is well known that the viscous term in the Navier–Stokes equation represents a singular perturbation in the limit of $Re\to \infty$. As a result, the solutions of the Euler equation

(3.1)\begin{equation} \partial_t\omega+ (\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla})\omega =0, \end{equation}

which describes the inviscid limit, may have properties that are dramatically different from those of the Navier–Stokes equation. In the inviscid limit, vorticity is materially conserved, giving rise to an infinite number of conserved quantities such as

(3.2)\begin{equation} I_n=\frac{1}{4{\rm \pi}^2}\int_0^{2{\rm \pi}}\int_0^{2{\rm \pi}}\omega^n\,{\rm d}x\,{\rm d}y \end{equation}

with integer $n\ge 1$ (Eyink 1996; Clercx & Van Heijst 2009), in addition to the conserved energy $E$. Hence, by Noether's theorem, the Euler equation has an infinite number of continuous symmetries and inviscid flows belong to infinite-dimensional solution families spanned by continuous parameters $\sigma _n$, $n=1,2,\ldots$, associated with the corresponding symmetries. On the other hand, the Navier–Stokes equation generally breaks all of these continuous symmetries (except for the temporal translation, so long as the forcing is time-independent), retaining only the discrete translational symmetries of the forcing, if there are any.

As a result of their different symmetries, the Euler equation has a much broader variety of solutions than the Navier–Stokes equation. Consequently, it is natural to ask which solutions of the Euler equation are physically relevant, i.e. correspond to a solution of the Navier–Stokes equation at a high, but finite $Re$. For steady flows with closed streamlines, this limit was originally considered by Batchelor (1956) who derived an integral condition on the stream function of the Euler flow. Okamoto (1994) has subsequently shown that Batchelor's criterion is equivalent to a solvability condition for this singularly perturbed problem.

This idea can be developed further to obtain several specific, interpretable constraints on the dynamically relevant solutions of the Euler equation, not only steady, but also time-periodic. Let $[\boldsymbol {u}_0,p_0]$ be a solution of the Euler equation and $[\boldsymbol {u},p]=[\boldsymbol {u}_0,p_0]+Re^{-1}[\boldsymbol {u}_1,p_1]$ be a solution of the Navier–Stokes equation (2.1) with $Re\gg 1$. The perturbation $[\boldsymbol {u}_1,p_1]$ satisfies the equation

(3.3)\begin{equation} \hat{N}[\boldsymbol{u}_1,p_1]=[\boldsymbol{H},0], \end{equation}

where we have defined

(3.4)\begin{equation} \left.\begin{gathered} \hat{N}[\boldsymbol{u}_1,p_1]=[\partial_t \boldsymbol{u}_1 + (\boldsymbol{u}_0\boldsymbol{\cdot}\boldsymbol{\nabla})\boldsymbol{u}_1+(u_1\boldsymbol{\cdot}\boldsymbol{\nabla})\boldsymbol{u}_0+\boldsymbol{\nabla} p_1,\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}_1],\\ \boldsymbol{H} = \nabla^2\boldsymbol{u}_0+\boldsymbol{h}, \end{gathered}\right\} \end{equation}

and assumed that $\boldsymbol {h}= Re\,\boldsymbol {f}$ is $O(1)$.

Subject to the proper boundary conditions in space (and, for time-periodic solutions, in time), (3.3) defines a boundary value problem for the perturbation $[\boldsymbol {u}_1,p_1]$. Since the linear operator $\hat {N}$ has null eigenvalues associated with every continuous symmetry of the Euler equation, where $\boldsymbol {e}_n=[\partial \boldsymbol {u}_0/\partial \sigma _n,\partial p_0/\partial \sigma _n]$ are the corresponding eigenfunctions, the boundary value problem is ill-posed and only has solutions provided the solvability condition

(3.5)\begin{equation} \langle\partial \boldsymbol{u}_0/\partial \sigma_n,\boldsymbol{H}\rangle=0 \end{equation}

is satisfied for each of the group parameters $\sigma _n$, where

(3.6)\begin{equation} \langle\boldsymbol{u},\boldsymbol{v}\rangle= \int_\varOmega(\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{v})\,{\rm d}\varOmega \end{equation}

and the integral is taken over the entire domain in space (and, for time-periodic solutions, time) on which the inviscid solution $[\boldsymbol {u}_0,p_0]$ is defined.

In particular, the scaling symmetry of the Euler equation implies that $\boldsymbol {u}_0|_{\sigma _t}=\textrm {e}^{\sigma _t}\boldsymbol {u}_0|_{\sigma _t=0}$ and $p_0|_{\sigma _t}=\textrm {e}^{2\sigma _t}p_0|_{\sigma _t=0}$ is a solution for any real $\sigma _t$. Hence, $\boldsymbol {e}_t=[\boldsymbol {u}_0,2p_0]_{\sigma _t=0}$, and the corresponding solvability condition reduces to

(3.7)\begin{equation} \langle\boldsymbol{u}_0,\boldsymbol{h}\rangle ={-}\langle\boldsymbol{u}_0,\nabla^2\boldsymbol{u}_0\rangle. \end{equation}

For $Re\gg 1$, we can replace the inviscid solution with the viscous one yielding, to leading order in $Re^{-1}$,

(3.8)\begin{equation} \langle\boldsymbol{u},\boldsymbol{f}\rangle={-}Re^{{-}1}\langle\boldsymbol{u},\nabla^2\boldsymbol{u}\rangle, \end{equation}

so the energy injection by the forcing has to balance the energy dissipation by viscosity at every instant for steady-state solutions and over one period for time-periodic solutions of the Navier–Stokes equation. While energy conservation clearly has to be satisfied for steady and time-periodic viscous flows, our numerical simulations show that, in the asymptotic regime, turbulent flow also satisfies this condition, with the energy becoming essentially time-independent or slowly varying, as figure 1(c) illustrates. Moreover, for smooth solutions of the Euler equation describing large-scale flows, $\langle \boldsymbol {u}_0,\nabla ^2\boldsymbol {u}_0\rangle =O(1)$. Hence, so long as the turbulent flow field is close to such an inviscid solution, we have $\langle \boldsymbol {u},\nabla ^2\boldsymbol {u}\rangle \approx \langle \boldsymbol {u}_0,\nabla ^2\boldsymbol {u}_0\rangle$, hence $\boldsymbol {u}$ remains essentially orthogonal to the forcing field.

Continuous translational symmetry in the $x$-direction corresponds to the Goldstone mode $\boldsymbol {e}_x=[\partial _x\boldsymbol {u}_0,\partial _xp_0]$ and yields the solvability condition

(3.9)\begin{equation} \langle\partial_x\boldsymbol{u}_0,\boldsymbol{h}\rangle ={-}\langle\boldsymbol{u}_0,\partial_x\boldsymbol{h}\rangle={-}\langle\partial_x\boldsymbol{u}_0,\nabla^2\boldsymbol{u}_0\rangle. \end{equation}

For a corresponding viscous solution $\boldsymbol {u}$, we find, to leading order in $Re^{-1}$,

(3.10)\begin{equation} \langle\boldsymbol{u},\partial_x\boldsymbol{f}\rangle ={-}Re^{{-}1}\langle\partial_x\boldsymbol{u},\nabla^2\boldsymbol{u}\rangle. \end{equation}

Again, since $\langle \partial _x\boldsymbol {u},\nabla ^2\boldsymbol {u}\rangle \approx \langle \partial _x\boldsymbol {u}_0,\nabla ^2\boldsymbol {u}_0\rangle =O(1)$ for viscous flows that are close to smooth solutions of the Euler equation, $\boldsymbol {u}$ should be (nearly) orthogonal to $\partial _x\boldsymbol {f}$. Coupled with the condition (3.8), this implies local continuous translational symmetry of the viscous flow field. Similar statements can be made regarding translational symmetry in the $y$-direction. Therefore, although forcing generally breaks continuous translational symmetry at finite $Re$, this symmetry should be effectively restored, at least for infinitesimal shifts, in the limit $Re\to \infty$ for the forcing profiles that do not significantly affect the structure of the large-scale flow. The same conclusions should apply to a turbulent flow field at high $Re$, so long as it stays in the neighbourhood of smooth inviscid solutions.

The implications of the solvability conditions for turbulent flow can be quantified by computing the rate of energy injection

(3.11)\begin{equation} P=\frac{1}{4{\rm \pi}^2}\int_0^{2{\rm \pi}}\int_0^{2{\rm \pi}}\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{f}\,{\rm d}x\,{\rm d}y \end{equation}

and the measure of local translational symmetry breaking

(3.12)\begin{equation} Q=\frac{1}{4{\rm \pi}^2N}\left|\int_0^{2{\rm \pi}}\int_0^{2{\rm \pi}}\sum_iu_i\boldsymbol{\cdot}\boldsymbol{\nabla} f_i\,{\rm d}x\,{\rm d}y\right|. \end{equation}

Figure 3 shows the results for a portion of turbulent trajectory, $300< t<400$, for which the large-scale flow is nearly time-periodic with the apparent period of $T\approx 10$, as indicated by the instantaneous value of $P$ in figure 3(a). Assuming there is a nearby smooth, time-periodic solution of the Euler equation, we should expect the average of $P$ computed over that period to be $O(Re^{-1})$ for all times, and this is what we find. The temporal average becomes almost constant, with the magnitude representing the balance between the energy injection and viscous dissipation. Moreover, the temporal averages of both $P$ and $Q$ shown in figures 3(b) and 3(c) are found to be $O(Re^{-1})$ not only for the turbulent flow $\boldsymbol {u}(\boldsymbol {x},t)$ but also for its shifted version $\boldsymbol {u}(\boldsymbol {x}-\boldsymbol {a},t)$, where the shift $\boldsymbol {a}$ is allowed to vary continuously over one period of the forcing in both directions, i.e. $-{\rm \pi} /4< a_x,a_y<{\rm \pi} /4$. This suggests that, in the limit $Re\to \infty$, the large-scale flow recovers continuous translational symmetry with respect to finite, not just infinitesimal, shifts.

Figure 3. The energy injection $P$ as a function of time (a). The instantaneous value is shown in blue and the running average computed over the apparent period of the large-scale flow is shown in orange. The temporal average of $P$ (b) and $Q$ (c) computed over the entire time interval $300< t<400$ as a function of the shift $\boldsymbol {a}$.

4. Numerical results

Temporal evolution of the large-scale flow can be conveniently analysed by inspecting a suitably defined recurrence function (Cvitanović & Gibson 2010; Lucas & Kerswell 2015):

(4.1)\begin{equation} R(t,\tau) = \min_{\boldsymbol{a}}\|\hat L_{16} [\boldsymbol{u}(\boldsymbol{x},t)- \boldsymbol{u}(\boldsymbol{x}+\boldsymbol{a},t-\tau)]\|. \end{equation}

To avoid explicit minimization over the shift $\boldsymbol {a}$, effective continuous translational symmetry of the large-scale flow was reduced using Fourier-mode slicing (Budanur et al. 2015). Specifically, we shifted $\boldsymbol {u}(\boldsymbol {x},t)$ in both spatial directions so that the phase $\phi _x$ ($\phi _y$) of the first Fourier mode $\boldsymbol {k}=(1,0)$ ($\boldsymbol {k}=(0,1)$) becomes zero before the filtering is applied and the norm is computed in (4.1). The components of the shift $\boldsymbol {a}$ were then reconstructed from the original phases, e.g. $a_x=\phi _x(t)-\phi _x(t-\tau )$. Deep (relative to the mean value) minima of $R(t,\tau )$ which correspond to $|\boldsymbol {a}|\ll 1$ represent segments of turbulent flow that come close to equilibria or time-periodic states, while deep minima with $|\boldsymbol {a}|=O(1)$ correspond to relative equilibria (travelling waves) or relative periodic orbits.

A representative example of such recurrence function is shown in figure 4. Note that qualitatively similar recurrence functions are produced by DNS on finer meshes (e.g. $2048\times 2048$). Deep minima of $R(t,\tau )$ at integer multiples of $\tau \approx 10$ which, for this time interval, correspond to $\boldsymbol {a}\approx 0$, suggest that the large-scale component of the turbulent flow closely follows (or shadows, in the language of dynamical systems) a time-periodic solution of the unforced Euler equation (3.1) with a period $T\approx 10$. The corresponding numerically exact solutions can be found using the tuples $\{\boldsymbol {u}(\boldsymbol {x},\boldsymbol {t}),\tau,\boldsymbol {a}\}$ representing the minima as initial conditions for a Newton-generalized minimal residual (Newton-GMRES) solver accounting for translational symmetry of the Euler equation. A detailed description of such a solver is provided, for instance, by Marcotte & Grigoriev (2015). In this study, our focus is entirely on smooth solutions describing large-scale flows. Hence, initial conditions should be prepared by smoothing the flow field $\boldsymbol {u}(\boldsymbol {x},t)$ to eliminate the fine structure, as discussed in the Appendix. Smoothing also helps speed up convergence and improve the success rate of the solver. For illustration, figures 2(b) and 2(d) show, respectively, snapshots of the vorticity field representing an initial condition (here the large-scale flow without additional smoothing) and the almost indistinguishable converged solution of the Euler equation.

Figure 4. Recurrence function $R(t,\tau )$ for a Fourier-filtered turbulent flow. The minima at multiples of $\tau \approx 10$ show that large-scale flow exhibits nearly periodic dynamics with a period $T\approx 10$ over extremely long intervals.

Another example of near-recurrences of turbulent flow is shown in figure 5, which corresponds to the zoomed-in version of figure 4. In particular, a dark blue triangular region at $t\approx 603$ and $\tau \lesssim 2$ represents a short interval when the large-scale flow is nearly stationary. The minima of the recurrence function at $605\lesssim t\lesssim 615$ and $\tau \approx 1$ represent a longer interval when the large-scale flow is nearly time-periodic. The corresponding solutions of the Euler equation and their properties are discussed in detail in the subsequent sections.

Figure 5. A zoomed-in version of figure 4 showing a shorter interval of nearly time-periodic dynamics with a period $T\approx 1$.

4.1. Equilibria

Equilibria are the simplest type of recurrent solutions. Perhaps the best-known example is Taylor–Green vortex equilibria which, in two spatial dimensions, are described by

(4.2)\begin{equation} \left.\begin{gathered} \psi=A\cos(\gamma(x+a_x))\cos(\gamma(y+a_y)),\\ \omega={-}A\gamma^2\cos(\gamma(x+a_x))\cos(\gamma(y+a_y)). \end{gathered}\right\} \end{equation}

Due to the translational symmetry of the Euler equation, the two constant phases $a_x$ and $a_y$ are arbitrary. The constant amplitude $A$ (the analogue of the energy $E$) as well as wavenumber $\gamma$ are also arbitrary due to the scaling symmetry of the Euler equation which, on an unbounded domain, implies that if $\boldsymbol {u}(\boldsymbol {x},t)$ (and $\omega (\boldsymbol {x},t)$) is a solution of the Euler equation, then so is $\gamma ^{-1}\lambda \boldsymbol {u}(\gamma \boldsymbol {x},\lambda t)$ (and $\gamma ^{-1}\lambda ^2\omega (\gamma \boldsymbol {x},t)$) for any $\lambda$ and $\gamma$. Spatial periodicity restricts $\gamma$ to integer multiples of 1 or $1/\sqrt {2}$, but $\lambda$ can take a continuum of values.

As discussed in § 3, each of the Taylor–Green vortices should belong to a family of equilibria spanned by an infinite number of continuous parameters including $a_x$, $a_y$ and $\lambda$. Other continuous parameters describe the shape of the vortices. These shape parameters correspond to Lie point symmetries of the Euler equation (Liu, Li & Zhao 2019). For instance, on an unbounded domain, a single circular vortex

(4.3)\begin{equation} \omega=F(r),\quad r=\sqrt{x^2+y^2} \end{equation}

is a solution to the steady-state version of the Euler equation

(4.4)\begin{equation} (\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla})\omega = 0 \end{equation}

for any choice of the shape function $F(r)$. Indeed, we can expand $F(r)$ in a basis of, say, Bessel functions, with the infinite set of Fourier–Bessel coefficients playing the role of continuous shape parameters.

For our purposes, it is convenient to classify different solution families using topology of their streamlines following Moffatt (1987). For instance, the Taylor–Green vortex with

(4.5)\begin{equation} \omega=\cos x+\cos y \end{equation}

is obtained by setting $\gamma =1/\sqrt {2}$ and $\boldsymbol {a}=0$ in (4.2) and rotating the coordinate system by ${\rm \pi} /4$. It belongs to a family of solutions featuring two counter-rotating vortices which also includes an infinite number of other equilibria unrelated by either translational or scaling symmetry.

These solutions can be found either by applying Newton-GMRES solver to flows time-integrated over a short time interval or by solving (4.4) directly using GMRES iterations for different initial conditions. For instance, using the state $\omega =(\cos x+\cos y)^4\,{\mathrm {sign}}(\cos x+\cos y)$ as a non-equilibrium initial condition, we found the equilibrium shown in figure 6(a). We then used homotopy (see the Appendix) to construct a continuous family of equilibria connecting that numerical solution with the analytical solution (4.5). Note that homotopy uses known smooth solutions to generate new smooth initial conditions, so no additional smoothing is needed prior to Newton-GMRES iterations. Representative intermediate equilibria belonging to this family are shown in figure 6(b,c). The members of this family feature two symmetric, counter-rotating vortices and differ mainly in the width of the vortices.

Figure 6. A family of symmetric equilibria, computed via homotopy between the state shown in panel (a) and a Taylor–Green vortex (4.5) (not shown). Intermediate states are shown in panels (b) and (c).

A different continuous family of equilibria can be obtained by breaking the reflection symmetry between the two vortices. For instance, adding a Gaussian patch of vorticity to the reflection-symmetric equilibrium state shown in figure 7(a) and then subtracting the mean amplifies one of the vortices and weakens and broadens the other. Reconverging the resulting (non-equilibrium) state using Newton-GMRES solver, we found the asymmetric equilibrium state shown in figure 7(c). We then used homotopy to construct another continuous family of solutions, with a representative intermediate equilibrium belonging to this family shown in figure 7(b).

Figure 7. A family of asymmetric equilibria, computed via homotopy between the states shown in panels (a) and (c). An intermediate state is shown in panel (b).

Which members of such continuous families, if any, are dynamically relevant for (e.g. are visited by) high-$Re$ turbulent flow needs to be established, however. As discussed previously, such dynamically relevant equilibria can be found using the Newton-GMRES solver seeded with initial conditions identified through recurrence analysis. Good initial conditions correspond to deep minima (compared with the mean value) of $R(t,\tau )$ with a fixed non-vanishing $\tau \ll 1$. One example of such an equilibrium is shown in figure 8(a); the corresponding initial condition was obtained by applying spectral and stream function smoothing to the snapshot of turbulent flow at $t\approx 603$ corresponding to figure 4. This solution is qualitatively quite similar to the equilibria shown in figures 6 and 7 which also speaks in favour of all these equilibria belonging to one family parameterized by multiple continuous parameters.

Figure 8. (a) A dynamically relevant equilibrium computed using a turbulent flow snapshot and (b) the corresponding spectrum of stability eigenvalues $\lambda _i$ computed on a $100\times 100$ grid. The horizontal (vertical) axis corresponds to the real (imaginary) part of $\lambda$. (c) The number $M$ of marginal stability eigenvalues found on an $n\times n$ grid.

The number of continuous parameters for any family will be reflected in the stability spectrum of any member of the family. Since generators of the Lie group describing different continuous parameters commute with the evolution operator, each continuous parameter would correspond to a marginal eigenvalue associated with a corresponding marginal mode. Consider, for instance, the dynamically relevant equilibrium shown in figure 8(a). Its stability spectrum computed in Fourier space using the MATLAB function eig() is shown in figure 8(b). Due to memory constraints, this computation has been performed on grids up to $100\times 100$. Since the total number and accuracy of the eigenvalues depends on both the spatial resolution of the flow field and the accuracy of eig(), it is difficult to estimate the number of continuous parameters precisely. Let us define, rather arbitrarily, an eigenvalue $\lambda _i$ to be marginal so long as $|\lambda _i|<10^{-5}\ll \max _i(Re(\lambda _i))\approx 0.06$. The number $M$ of such marginal eigenvalues, computed on an $n\times n$ grid, grows roughly linearly with $n$, as illustrated by figure 8(c), suggesting that $M\to \infty$ as $n\to \infty$. Note that the limited accuracy of eig() is also responsible for the deviations from the expected inversion symmetry $\lambda \to -\lambda$ in figure 8(b) which reflects the time-reversal symmetry of the Euler equation.

The stability spectrum also shows that the equilibrium is only moderately unstable. Unstable perturbations grow on a time scale of $\max (Re(\lambda ))^{-1}\approx 17$ comparable to the characteristic time scale $T_c\approx 10$ associated with turbulent flow. At first sight, this result appears very surprising. Indeed, this equilibrium describes inviscid flow and therefore corresponds to the limit $Re\to \infty$, where ECSs are expected to be very strongly unstable. While that may indeed be the case for recurrent solutions describing small-scale flows (i.e. coherent substructures), the presence of the inverse cascade in two spatial dimensions implies that recurrent solutions describing large-scale flows (i.e. coherent structures) can only be unstable with respect to perturbations with low spatial frequencies. These low-frequency modes are largely aligned with the infinite-dimensional subspace spanned by the marginal modes for any large-scale Euler equilibrium. This near-alignment is responsible for the relatively small real parts of the eigenvalues associated with unstable eigenmodes.

Overall, we find that although the large-scale turbulent flow does occasionally visit the neighbourhoods of equilibrium solutions of Euler, recurrence analysis shows that these instances are quite rare. In order to understand why that is the case, we computed the leading eigenmodes associated with a pair of representative equilibrium states. The most unstable eigenmode corresponding to the symmetric equilibrium from figure 7(a) is shown in figure 9(a). It is almost indistinguishable from the marginal mode associated with spatial translation in the direction indicated by the black arrow. Indeed, numerically evolving the equilibrium, we find that it transitions to a travelling wave moving in the corresponding direction. The similarity of the leading eigenmode to the translational marginal mode implies that the relative position of the vortices for the travelling wave does not change in time and hence neither does the flow generated by the vortex array in the comoving reference frame. As a result, the direction in which the vortices travel never changes either.

Figure 9. The leading mode of instability for equilibria calculated using power iteration. (a) The leading eigenfunction for the equilibrium shown in figure 7(a) represents an infinitesimal shift associated with translation of the flow pattern in the direction shown by the arrow. (b) The leading Floquet mode for the equilibrium shown in figure 7(b) represents an infinitesimal shift of the vortex cores associated with circular motion of the flow pattern indicated by the arrows. The snapshot shown corresponds to a particular phase of the oscillation.

For the asymmetric equilibrium from figure 7(b), the situation is quite different. The leading eigenvalues in this case form a complex conjugate pair, so the corresponding eigenmode is time-periodic; its snapshot is shown in figure 9(b). The spatial structure of this mode is quite interesting and sheds some light on the dynamics of both large- and small-scale flow structures in turbulence. In the regions associated with the vortex cores, the eigenfunction is similar to the translational marginal mode, and represents a shift of the vortices. As time progresses, the corresponding dipolar structures rotate in the direction indicated by the arrows, which suggests that this equilibrium evolves into a time-periodic flow, with the vortex centres executing a circular motion, as confirmed by numerical simulation.

The structure of the eigenmode in the hyperbolic regions does not correspond to a rigid spatial shift, as can be seen by comparison with figure 9(a), suggesting that the flow pattern not only translates but also deforms. Instead, the eigenmode features thin vorticity filaments aligned with the streamlines of the base flow. These vorticity filaments are qualitatively similar to those found in turbulent flow (cf. figure 2c). This observation suggests that it is the time-dependence of the large-scale flow that is responsible for generation of small scales and the direct cascade overall.

Since symmetric equilibria are rare compared with asymmetric ones (they represent a set of measure zero in the corresponding continuous parameter space), we should expect equilibria to transition predominantly to time-periodic flows. Indeed, the recurrence function (cf. figure 4) illustrates that the large-scale flow is nearly time-periodic over a significant fraction of the time, but only rarely resembles travelling waves. Both types of solutions are discussed in more detail below.

4.2. Travelling waves

While travelling waves appear to represent time-periodic flows on a domain with periodic boundary conditions, they become equilibria in a comoving reference frame, i.e. they correspond to relative equilibria. Absolute equilibria discussed in the previous section correspond to the special case of the comoving frame having a zero velocity. Since the velocity of a generic comoving frame is quite unlikely to vanish, equilibria are expected to be far less common than travelling waves, which explains why time-independent large-scale flows are so rarely observed in the DNS of turbulence.

Travelling waves and equilibria of the Navier–Stokes equation are distinct and the velocity of the comoving frame for travelling waves can only take discrete values (Chandler & Kerswell 2013; Lucas & Kerswell 2014, 2015). In contrast, travelling waves and equilibria of the Euler equation are connected, and the velocity of the comoving frame can take a continuum of values. This is illustrated in figure 10, which describes a continuous family of travelling wave solutions obtained with the help of homotopy between an equilibrium shown in figure 10(a), which belongs to the family shown in figure 6, and a travelling wave shown in figure 10(b). The initial condition for the travelling wave was generated by shifting one of the vortices. Specifically, the vorticity field was updated according to $\omega \to \omega +\varepsilon \boldsymbol {n}\boldsymbol {\cdot }\boldsymbol {\nabla }\omega$ in the regions with $\omega >0$, where $\varepsilon \ll 2{\rm \pi}$ and $\boldsymbol {n}=(1,-1)$, followed by spectral smoothing, and then integrating the flow in time until visible transients disappeared. While the corresponding vorticity fields for these two solutions appear extremely similar, they are slightly different. Their difference is described by the marginal mode associated with this continuous family which is shown in figure 10(c) and represents a relative displacement of the two vortices. Just like arrays of point vortices, equilibrium arrangements of extended vortices require high symmetry balancing advection flows generated by actual or virtual neighbours. When the symmetry of a vortex array is broken, advection causes an overall drift in the pattern. As shown in figure 10(d), the drift velocity $v$ increases with the degree of asymmetry quantified by the normalized state space distance

(4.6)\begin{equation} d(\boldsymbol{u})=\frac{\|\boldsymbol{u}-\boldsymbol{u}_0\|}{\|\boldsymbol\|{u_0}} \end{equation}

from the equilibrium $\boldsymbol {u}_0$.

Figure 10. A continuous family of solutions connecting an equilibrium (a) and a travelling wave solution (b). The marginal mode representing this solution family (c). The speed of the comoving frame as a function of the state space distance to the equilibrium (d).

In the example considered here, the two vortices have the same shape and strength and drift with the same velocity, yielding a travelling wave. When the shapes and/or strengths of the two vortices are different, their drift speeds will not be the same, yielding a gradual change in the relative position of the vortices and therefore a change in the direction of the drift. Hence, asymmetric vortices will generally have dynamics that are more complicated than travelling waves; one example is time-periodic states discussed in the next section. We can again invoke the argument that symmetric vortex arrangements are rare, so ECSs in the form of travelling waves should also be relatively uncommon in turbulence. This prediction is indeed born out by recurrence analysis.

We will conclude our discussion of travelling waves by noting that they, just like absolute equilibria, are expected to belong to families with an infinite number of continuous parameters, some of which correspond to translational and scaling symmetry while the rest describe the shape and relative position of the vortices. The non-trivial example shown in figure 10 represents continuous variation in the relative position for vortices with a particular shape. However, given the arbitrary choice of the vortex shape in this example, vortices with any other shape can also be continued to corresponding travelling waves by varying their relative position.

4.3. Periodic orbits

As illustrated by the recurrence function shown in figure 4, turbulent flow exhibits long intervals during which large-scale flow is almost time-periodic in some reference frame. The shift $\boldsymbol {a}$ during these intervals tends to be quite small, suggesting that the comoving reference frame is stationary and, therefore, the large-scale flow is well-described by an unstable periodic orbit solution (UPO) of Euler. The recurrence function shows that there are two distinct types of UPOs characterized by substantially different periods: $T\approx 10$ or $T\approx 1$. A representative example of UPO with period 10.017 is shown in figure 2(d), and a corresponding movie is provided as supplementary movies. This type of UPO features a pair of axially symmetric vortices, just like the equilibria and travelling waves discussed previously. The shapes of the vortices do not change noticeably over the period, while their centres undergo a nearly circular motion.

A representative example of UPO with the shorter period is shown in figure 11(b); a corresponding movie is provided as supplementary movies. The initial condition was obtained by applying a combination of hyperviscous smoothing and spectral smoothing to a snapshot of a turbulent flow. This type of UPO features a pair of vortices only one of which is axially symmetric. The other vortex has a strongly asymmetric shape, e.g. elliptical or, as in the present example, tripolar. We find that the centres of both vortices are essentially stationary, so the time-dependence is associated entirely with the rigid rotation of the asymmetric vortex. Coherent structures featuring tripolar vortices in particular have been frequently observed in both numerical simulations of 2-D turbulence (Legras, Santangelo & Benzi 1988) and experiments in rotating tanks (van Heijst, Kloosterziel & Williams 1991). Our results show that such coherent structures correspond to numerically exact solutions of the Euler equation.

Figure 11. Snapshot of turbulent flow $\omega$ (a), the corresponding large-scale flow $\hat {L}_{16}\omega$ (b) and the converged UPO with the period $T\approx 1.05$ (c). The arrow indicates the direction of rotation for the tripolar vortex.

This is illustrated by figure 11 which compares a snapshot of DNS of 2-D turbulence featuring a tripolar vortex with the corresponding large-scale flow and numerically converged UPO of Euler. It is important to note that smoothing affects the smaller-scale structures of the flow, and Newton-GMRES solver is not guaranteed to converge to the ECS that is the closest to the initial condition, so some difference in the vorticity fields shown in figures 11(b) and 11(c) is expected. It is worth noting that solutions featuring a single tripolar vortex on an unbounded spatial domain correspond to relative equilibria of Euler, e.g. equilibria in a corotating reference frame. In the present case, time-periodicity arises from breaking of the rotational symmetry of Euler by both the boundary conditions and the presence of another (stationary and axially symmetric) vortex.

Just like equilibrium and travelling wave solutions of the Euler equation, UPOs are found to come in families parameterized by numerous continuous parameters, some of which correspond to translational and scaling symmetry of Euler while others describe the shape and relative position of the vortices. Figure 12 shows one such family of UPOs; a corresponding movie is provided as supplementary movies. A reference time-periodic solution shown in figure 12(b) was obtained by further smoothing the flow shown in figure 2(b). The solution family was constructed by converging initial conditions obtained by modifying the vorticity field of the reference solution according to $\omega \to |\omega |^\kappa {\mathrm {sign}}(\omega )$ over a continuous range of parameter $\kappa$. Solutions corresponding to $\kappa =0.5$ and $\kappa =1.8$ are shown in figure 12(a,c), respectively. This procedure is an alternative to the homotopy that allows constructing a continuous family of solutions using one reference state rather than two. The continuous nature of this solution family is seen in figure 12(d) which shows the period $T$ and the amplitude

(4.7)\begin{equation} A=\frac{1}{T}\int_0^T\sqrt{(\phi_x(t)-\overline{\phi_x})^2 +(\phi_y(t)-\overline{\phi_y})^2}\,{\rm d}t, \end{equation}

of the UPO, which quantifies the degree of unsteadiness of the flow, as a function of the enstrophy $H=I_2$.

Figure 12. A continuous family of UPOs which reflects variation in the shape of the vortices. Panels (a–c) show snapshots of representative states from the same family. (d) Dependence of the amplitude and period of UPOs on the normalized enstrophy for this solution family.

Just like travelling wave solutions, for the Euler equation, UPOs are not isolated from the other types of solutions. Figure 13 illustrates this for a continuous family of solutions connecting an equilibrium shown in figure 13(a) and a UPO shown in figure 13(b); a corresponding movie is provided as supplementary movies. The equilibrium featuring two asymmetric vortices was obtained through the same approach as that used to generate the state shown in figure 7(c) and the UPO was obtained from the equilibrium through the same procedure as that used to generate the state shown figure 10(c). The marginal mode describing this family is shown in figure 13(c) and represents a relative shift of the vortices. The equilibrium here, unlike that shown in figure 10, features vortices with different shapes. As discussed previously, the shape asymmetry leads to the drift of the vortices in a direction that changes with time. The vortex pattern moves in a circular path, as indicated by an arrow in figure 13(b), repeating after period $T$.

Figure 13. A continuous family of UPOs connecting the equilibrium shown in panel (a) and a UPO shown in panel (b). The arrow represents the direction in which the pattern moves. (c) The marginal mode associated with this solution family. (d) The period and amplitude of UPOs as a function of the state space distance to the equilibrium.

All UPOs characterized by the overall drift of the pattern that we found are formally preperiodic: the state at time $t=T/4$ is related to the state at time $t=0$ by a ${\rm \pi} /2$ rotation about some point in the domain. As a result, the pattern moves in an almost circular trajectory (an example is shown in figure 14 in red). The amplitude $A$ effectively describes the radius of the circle. For the family of UPOs shown in figure 13, the amplitude is shown in figure 13(d) and is found to vanish at the corresponding equilibrium and grow monotonically as the distance $d$ from that equilibrium increases. As the supplementary movies show, both vortices move along nearly circular trajectories with somewhat different radii, but the vortex pattern remains qualitatively the same at all times.

Figure 14. The position of the vortex pattern describing turbulent large-scale flow (in black) corresponding to figure 4 over the time interval $77\lesssim t\lesssim 400$. Subintervals of recurrent (non-recurrent) dynamics are shown as solid (dashed) lines. The trajectory of a representative UPO, computed from an initial condition along the recurrent portion of the flow, is shown in red.

To illustrate that families of UPOs characterized by such circular motion do indeed describe the dynamics of turbulent flow, we have plotted the position of the vortex pattern describing the large-scale flow, as characterized by the first-Fourier-mode phases $\phi _x$ and $\phi _y$, in figure 14. The trajectory of the turbulent flow in this low-dimensional projection of the state space is nearly time-periodic over the time intervals $77\lesssim t\lesssim 200$ and $270\lesssim t\lesssim 400$. These intervals are separated by an interval $200\lesssim t\lesssim 270$ during which the dynamics are aperiodic. Indeed, during the time-periodic intervals, the pattern moves in a (nearly) circular trajectory (shown as a solid black line) which closely resembles the trajectory of particular UPOs (shown as red circles). During the aperiodic interval, the trajectory (shown as a dashed black line) becomes much more complicated and is not described by any UPO. The observation that the phases $\phi _x$ and $\phi _y$ vary continuously in time instead of taking discrete values is consistent with the analysis of § 3 which suggests that, at high $Re$, turbulent flow recovers continuous translational symmetry.

As the recurrence function shown in figure 4 further illustrates, the temporal period of the large-scale flow is not constant but varies slowly, and so does the amplitude $A$. Therefore, turbulent flow follows different UPOs. This variation represents a slow drift in the infinite-dimensional parameter space describing a family of the UPOs. Due to the scaling invariance of the Euler equation, one of these continuous parameters is the energy $E$; the energy of dynamically relevant UPOs of the Euler equation should be equal to the energy of the large-scale turbulent flow. Recall that the Euler equation conserves energy, so $E$ is constant for any solution, including UPOs. Since the period of any UPO scales with its energy, $T\propto E^{-1/2}$, a slow variation in $E$ should cause a slow variation in the temporal period. However, the energy is not the only parameter that varies and affects the period.

In particular, consider the variation in the enstrophy $H$ of the turbulent flow shown as a continuous blue line in figure 15, computed after the low-pass filter has been applied, during the time interval where the dynamics are nearly time-periodic. We used appropriately smoothed snapshots of the turbulent flow as initial conditions for the Newton-GMRES solver to compute the corresponding UPOs. The enstrophy describing these UPOs (shown as black circles) is found to track that of the turbulent flow, providing further evidence that the evolution of large-scale flow at high but finite $Re$ can be well-described as a slow drift in the infinite-dimensional space of parameters describing continuous families of unstable time-periodic solutions of Euler.

Figure 15. The normalized enstrophy $H$ of turbulent large-scale flow (in blue) corresponding to the interval $t\in [300,400]$ in figure 4. The symbols correspond to UPOs converged from turbulent flow at the corresponding points in time.

In conclusion of this section, let us address the issue of stability of dynamically relevant UPOs which explains why large-scale flow can remain nearly time-periodic over time intervals far exceeding the characteristic time scale $T_c$. Computing stability spectra of UPOs of the Euler equation with meaningful precision is extremely challenging due to the very large number of marginal (unit) Floquet multipliers associated with various continuous parameters and the absence of a spectral gap with unstable multipliers. As a result, standard approaches such as Arnoldi iterations become intractable. To assess the stability of a representative UPO $\boldsymbol {u}(\boldsymbol {x},t)$ with period $T=12.23$, (whose snapshot is shown in figure 16a), we computed the evolution of a random perturbation by time-integrating the Euler equation linearized about this UPO. The magnitude of an infinitesimal perturbation $\delta \boldsymbol {u}$ at integral multiples of the period is shown in figure 16(b). The growth is found to be algebraic rather than exponential over 30+ periods, with the growth rate steadily decreasing over time. This suggests that the leading Floquet multiplier is extremely close to unity and we simply observe transient growth associated with the non-normality of the evolution operator. Additional validation of this conclusion is provided by the spatial profile of the perturbation which, at long times, becomes virtually indistinguishable from that of the marginal mode $\partial _t\boldsymbol {u}$ associated with time translation (cf. figure 16c,d). The corresponding perturbation essentially represents a slow drift in the period of the UPO.

Figure 16. Stability of a dynamically relevant periodic orbit shown in panel (a). Transient growth in the amplitude of a random infinitesimal perturbation (b). The shape of the perturbation at long times (c) is nearly identical to the marginal mode $\partial _t \boldsymbol {u}$ (d).

5. Discussion

We found that sustained high-$Re$ turbulence in a spatially periodic 2-D domain driven by stationary forcing, for extended intervals of time, exhibits nearly time-periodic dynamics on large scales. This observation has important implications for our understanding of the direct cascade. During these intervals, the large-scale flow is well-described by particular weakly unstable time-periodic solutions of the Euler equation. These UPOs all feature a pair of counter-rotating vortices separated by hyperbolic regions of extensional flow which contain most, if not all, of the small-scale structure. The stretching and folding of thin vorticity filaments due to, respectively, the advection by the large-scale flow and its time-dependent nature are a familiar mechanism which generates small-scale structure in chaotic advection of passive tracers by time-periodic flows (Cartwright, Feingold & Piro 1996; Grigoriev 2005). Indeed, the analogy between small-scale vorticity and passive tracers whose evolution is described by formally the same equation, has been noted a long time ago by Kraichnan (1975). However, the realization that flow at a single (e.g. large) scale could be solely responsible for generating a fractal structure characterized by power-law scaling with a non-integral exponent is much more recent (Ott 1999).

The presence of large-scale coherent structures is well known (Clercx & Van Heijst 2009; Boffetta & Ecke 2012) to lead to substantial deviations from the predictions of KLB theory which assumes implicitly or explicitly (Kraichnan 1975) that small scales add up incoherently. Indeed the presence of large-scale structure implies strong coherence at small scales as well. As our results illustrate, the relative position of the vortices and hence the orientation of the extensional flow remains nearly constant, implying strong correlations between the orientation of vorticity filaments at all scales. These strong correlations explain why the observed power law exponent $\alpha$ deviates so strongly from the KLB prediction for the turbulent flow considered here.

Rather counter-intuitively, while the properties of the large-scale flow, which control the dynamics of small-scale vorticity and hence the direct cascade, tend to change slowly and smoothly, the scaling exponent changes rather abruptly and unpredictably, as illustrated by figure 1(b). This observation echoes a similar one made by Ottino (1989) who noticed that transport properties of time-periodic 2-D flows can change in a rather irregular fashion when their parameters are varied smoothly. Transport and mixing of passive scalars is a result of Lagrangian chaos, and chaotic sets are known to exhibit irregular dependencies on parameters. Perhaps the best-known example of such irregular dependence is provided by the windows of periodic dynamics of the logistic map (characterized by poor mixing) interspersed within the chaotic parameter interval (characterized by good mixing).

6. Conclusions

The present study identified several classes of numerically exact unstable recurrent solutions of the Euler equation on a 2-D square domain with periodic boundary conditions: equilibria; travelling waves; and periodic orbits. All three types of solutions feature a pair of extended counter-rotating vortices. These solutions are the infinite-$Re$ analogues of the so-called ‘unimodal’ stationary and time-periodic solutions of the Navier–Stokes equation at high $Re$ (Kim & Okamoto 2010; Gallet & Young 2013; Kim & Okamoto 2015; Kim et al. 2017). Unlike the vast majority of exact solutions of 2-D Euler reported previously, many of the solutions reported here are dynamically relevant: they provide a reasonably accurate description of the spatial structure and dynamics of the large-scale component of high-Reynolds number turbulent flow. Hence, they should be thought of as the proper exact coherent structures in the limit of infinite Reynolds number. The time-periodic solutions were found to be particularly important: large-scale turbulent flow was found to mirror this type of recurrent solutions of the Euler equation surprisingly frequently.

Another important result of this study is that recurrent solutions of Euler belong to families parameterized by an infinite number of continuous parameters. Some of these parameters correspond to translational and scaling symmetries of Euler, while the rest can be thought of as ‘shape’ parameters describing the spatial profile and relative position of the vortices. Clearly, only a fraction of the members of those solution families are dynamically relevant. What distinguishes the region of the huge parameter space that is dynamically relevant from the rest is yet to be determined, although solvability conditions can provide useful insight.

Rather unexpectedly, recurrent solutions of Euler featuring a pair of counter-rotating vortices were found to be only weakly unstable. This is in contrast with the intuition generated by previous numerical studies of weakly turbulent flows which suggests that ECSs become progressively more unstable as the Reynolds number increases. The likely reason this intuition breaks down is that the recurrent solutions of Euler investigated here are analogues of classical coherent structures describing large-scale flows. On the other hand, continuation of recurrent solutions describing weakly turbulent flows to higher Reynolds numbers typically yields coherent substructures describing small-scale flows. It is these substructures that become increasingly unstable.

Another interesting fundamental question that remains unresolved is the relation between the spatial and temporal structure of the prevailing time-periodic ECSs and the enstrophy flux towards small scales (direct cascade). The KLB theory completely fails to describe the turbulent flow studied here, yet we still find a power-law scaling of the energy spectrum, albeit with a non-integer exponent. Subsequent work should determine why a power law emerges and what the exponent is.