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 (Reference Hussain1986) 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 Reference Viswanath2007) 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 Reference Waleffe1998, Reference Waleffe2001; Gibson, Halcrow & Cvitanović Reference Gibson, Halcrow and Cvitanović2008; Itano & Generalis Reference Itano and Generalis2009; Shekar & Graham Reference Shekar and Graham2018); 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. Reference Suri, Tithof, Grigoriev and Schatz2018, Reference Suri, Kageorge, Grigoriev and Schatz2020) and three spatial dimensions (Krygier, Pughe-Sanford & Grigoriev Reference Krygier, Pughe-Sanford and Grigoriev2021; Crowley et al. Reference Crowley, Pughe-Sanford, Toler, Krygier, Grigoriev and Schatz2022).
Exact coherent structures have already generated significant insight into weakly turbulent flows (Kawahara, Uhlmann & Van Veen Reference Kawahara, Uhlmann and Van Veen2012; Graham & Floryan Reference Graham and Floryan2021). Most notably, ECSs are found to capture self-sustaining processes that maintain wall-bounded turbulence (Waleffe Reference Waleffe1997; Hall & Smith Reference Hall and Smith1991). Exact coherent structures also elucidate the transition from laminar flow to turbulence explaining both the ‘bypass’ mechanism in linearly stable flows (Khapko et al. Reference Khapko, Kreilos, Schlatter, Duguet, Eckhardt and Henningson2016) and the formation of chaotic sets underpinning sustained turbulence (Kreilos & Eckhardt Reference Kreilos and Eckhardt2012). 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 Reference Batchelor1956; Okamoto Reference Okamoto1994; Gallet & Young Reference Gallet and Young2013).
One of the most notable successes in continuing recurrent solutions to higher 
$Re$ is the computation of attached eddies (Eckhardt & Zammert Reference Eckhardt and Zammert2018; Doohan, Willis & Hwang Reference Doohan, Willis and Hwang2019; Yang, Willis & Hwang Reference Yang, Willis and Hwang2019; Azimi & Schneider Reference Azimi and Schneider2020) and their bulk analogues (Deguchi Reference Deguchi2015; Eckhardt & Zammert Reference Eckhardt and Zammert2018) 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 Reference Kim and Okamoto2010, Reference Kim and Okamoto2015; Kim, Miyaji & Okamoto Reference Kim, Miyaji and Okamoto2017) and three spatial dimensions (Wang, Gibson & Waleffe Reference Wang, Gibson and Waleffe2007). In principle, recurrent large-scale flows could also be computed with the help of large eddy simulations (Rawat et al. Reference Rawat, Cossu, Hwang and Rincon2015; Hwang, Willis & Cossu Reference Hwang, Willis and Cossu2016), 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 Reference Lamb1924), translating pairs of counter-rotating vortices (Pierrehumbert Reference Pierrehumbert1980; Saffman & Tanveer Reference Saffman and Tanveer1982), rotating arrays of two or more vortices (Dritschel Reference Dritschel1985; Carton, Flierl & Polvani Reference Carton, Flierl and Polvani1989; Crowdy Reference Crowdy2002b) and stationary multipolar vortex arrays (Crowdy Reference Crowdy1999, Reference Crowdy2002a; Tur & Yanovsky Reference Tur and Yanovsky2004; Xue, Johnson & McDonald Reference Xue, Johnson and McDonald2017).
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. Reference Aref, Newton, Stremler, Tokieda and Vainchtein2003). 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,
\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$,
\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
\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 Reference Gallet and Young2013; Kim et al. Reference Kim, Miyaji and Okamoto2017). Note that, unlike, say, a Kolmogorov forcing profile (Tithof et al. Reference Tithof, Suri, Pallantla, Grigoriev and Schatz2017), a checkerboard profile breaks continuous translational symmetry in both spatial directions, constraining the types of ECSs that can be found at moderate 
$Re$ (Suri Reference Suri2021). 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 (Reference Kraichnan1967), Leith (Reference Leith1968) and Batchelor (Reference Batchelor1969). This discrepancy is well known (Boffetta & Ecke Reference Boffetta and Ecke2012) 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 Reference Kraichnan1967, Reference Kraichnan1971). 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.
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. Reference Chen, Ecke, Eyink, Wang and Xiao2003).
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
\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 Reference Smith and Yakhot1993; Bracco et al. Reference Bracco, McWilliams, Murante, Provenzale and Weiss2000), 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
\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
\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 Reference Eyink1996; Clercx & Van Heijst Reference Clercx and Van Heijst2009), 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 (Reference Batchelor1956) who derived an integral condition on the stream function of the Euler flow. Okamoto (Reference Okamoto1994) 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
\begin{equation} \hat{N}[\boldsymbol{u}_1,p_1]=[\boldsymbol{H},0], \end{equation}where we have defined
\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
\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
\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
\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}$,
\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
\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}$,
\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
\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
\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.
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 Reference Cvitanović and Gibson2010; Lucas & Kerswell Reference Lucas and Kerswell2015):
\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. Reference Budanur, Cvitanović, Davidchack and Siminos2015). 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 (Reference Marcotte and Grigoriev2015). 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.
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.
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
\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 Reference Liu, Li and Zhao2019). For instance, on an unbounded domain, a single circular vortex
\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
\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 (Reference Moffatt1987). For instance, the Taylor–Green vortex with
\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.
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).
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.
(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.
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 Reference Chandler and Kerswell2013; Lucas & Kerswell Reference Lucas and Kerswell2014, Reference Lucas and Kerswell2015). 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
\begin{equation} d(\boldsymbol{u})=\frac{\|\boldsymbol{u}-\boldsymbol{u}_0\|}{\|\boldsymbol\|{u_0}} \end{equation}
from the equilibrium 
$\boldsymbol {u}_0$.
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 Reference Legras, Santangelo and Benzi1988) and experiments in rotating tanks (van Heijst, Kloosterziel & Williams Reference van Heijst, Kloosterziel and Williams1991). Our results show that such coherent structures correspond to numerically exact solutions of the Euler equation.
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
\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$.
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$.
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.
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.
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.
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 Reference Cartwright, Feingold and Piro1996; Grigoriev Reference Grigoriev2005). 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 (Reference Kraichnan1975). 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 Reference Ott1999).
The presence of large-scale coherent structures is well known (Clercx & Van Heijst Reference Clercx and Van Heijst2009; Boffetta & Ecke Reference Boffetta and Ecke2012) to lead to substantial deviations from the predictions of KLB theory which assumes implicitly or explicitly (Kraichnan Reference Kraichnan1975) 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 (Reference Ottino1989) 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 Reference Kim and Okamoto2010; Gallet & Young Reference Gallet and Young2013; Kim & Okamoto Reference Kim and Okamoto2015; Kim et al. Reference Kim, Miyaji and Okamoto2017). 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.
Supplementary movies
Supplementary movies are available at https://doi.org/10.1017/jfm.2023.584.
Funding
This material is based upon work supported by the National Science Foundation under grant no. 2032657.
Declaration of interests
The authors report no conflict of interest.
Data availability statement
Exact coherent structures and their key properties as well as MATLAB codes for solving the governing equations and Newton-GMRES solvers are available on GitHub at https://github.com/cdggt/euler2D/tree/main.
Appendix A. Numerical methods
A.1. Solving the Euler and Navier–Stokes equations
Both the Euler equation and the Navier–Stokes equation can be written entirely in terms of the stream function 
$\psi =-\nabla ^{-2}\omega$,
\begin{equation} \partial_t\psi=g(\psi), \end{equation}
with an appropriately defined function 
$g(\psi )$. In the case of the Navier–Stokes equation,
\begin{equation} g(\psi)=\nabla^{{-}2}\left(\partial{\psi}{x}{}\partial{}{y}{}\nabla^2\psi-\partial{\psi}{y}{}\partial{}{x}{}\nabla^2\psi -\varphi(\boldsymbol{x},\boldsymbol{y})\right)+\frac{1}{Re}\nabla^2\psi, \end{equation}
where 
$\varphi$ is the forcing profile introduced in (2.2).
To handle a wide range of time scales present in turbulence, time stepping was performed using variable time step Runge–Kutta–Fehlberg scheme
\begin{equation} \psi(t+\Delta t) = \psi(t) + \Delta t\left(\frac{47}{450}k_1+\frac{12}{25}k_3+\frac{32}{225}k_4+\frac{1}{30}k_5+\frac{6}{25}k_6\right), \end{equation}where
\begin{equation} \left.\begin{gathered} k_1 = g\left(\psi(t)\right) , \\ k_2 = g\left(\psi(t)+\frac{2\Delta t}{9}k_1\right) , \\ k_3 = g\left(\psi(t)+\frac{1\Delta t}{12}k_1+\frac{1\Delta t}{4}k_2\right) ,\\ k_4 = g\left(\psi(t)+\frac{69\Delta t}{128}k_1-\frac{243\Delta t}{128}k_2+\frac{135\Delta t}{64}k_3\right) , \\ k_5 = g\left(\psi(t)-\frac{17\Delta t}{12}k_1+\frac{27\Delta t}{4}k_2-\frac{27\Delta t}{5}k_3+\frac{16\Delta t}{15}k_4\right), \\ k_6 = g\left(\psi(t)+\frac{65\Delta t}{432}k_1-\frac{5\Delta t}{16}k_2+\frac{13\Delta t}{16}k_3+\frac{4\Delta t}{27}k_4+\frac{5\Delta t}{144}k_4\right). \end{gathered}\right\} \end{equation}
The step size 
$\Delta t$ was adjusted according to
\begin{equation} \Delta t_{new} =\epsilon \mathcal{E}^{{-}1/5}\Delta t_{old} \end{equation}
with fixed tolerance 
$\epsilon$ (typically 
$\epsilon =1.5$) and
\begin{equation} \mathcal{E} = \Delta t \left\|\frac{1}{360}k_1-\frac{128}{4275}k_3-\frac{2197}{75240}k_4+\frac{1}{50}k_5+\frac{2}{55}k_6 \right\|. \end{equation}In the case of the Euler equation, where
\begin{equation} g(\psi) = \nabla^{{-}2}\left(\partial{\psi}{x}{}\partial{}{y}{}\nabla^2\psi-\partial{\psi}{y}{}\partial{}{x}{}\nabla^2\psi\right), \end{equation}time-dependent solutions were computed using a fixed-step-size fourth-order Runge–Kutta method,
\begin{equation} \psi(t+\Delta t) = \psi(t)+\Delta t\left(\frac16k_1+\frac13k_2+\frac13k_3+\frac16k_4\right), \end{equation}where
\begin{equation} \left.\begin{gathered} k_1 = g(\psi(t)), \\ k_2 = g\left(\psi(t)+\frac{\Delta t}{2}k_1\right), \\ k_3 = g\left(\psi(t)+\frac{\Delta t}{2}k_2\right), \\ k_4 = g(\psi(t)+\Delta t k_3). \end{gathered}\right\} \end{equation}
To prevent contamination of smooth solutions by high-frequency noise associated with discretization and truncation errors, we added a term 
$\hat {D}\psi$ representing hyperviscosity to 
$g(\psi )$, where
\begin{equation} \hat D\, {\rm e}^{{\rm i}\boldsymbol{k}\,\boldsymbol{\cdot}\,\boldsymbol{x}} =\begin{cases} \alpha(k_v-k)\exp\left(\dfrac{\beta}{k_v-k} \right) {\rm e}^{{\rm i}\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{x}}, & \text{if } k>k_v,\\ 0, & \text{otherwise}. \end{cases} \end{equation}
We used the following choices of parameters: 
$k_v=40$; 
$\alpha =0.5$; and 
$\beta =1/13$. We also checked that converged solutions change negligibly if 
$k_v$ is increased to 64. The choice of the functional form of the hyperviscous term is arbitrary, so long as it has no effect on wavenumbers below the threshold 
$k_v$, while the high wavenumbers are strongly suppressed. For our choice of 
$k_v$, hyperviscosity only affects the high frequencies that are essentially absent in large-scale flows.
Spatial derivatives were all computed spectrally for both the Euler and Navier–Stokes equation. Additionally, a 
$2/3$ dealiasing scheme was implemented for numerical stability. Typically, we computed solutions to the Navier–Stokes equation on a 
$512\times 512$ grid, although we also checked the key results for consistency using grids up 
$2048\times 2048$. A time-averaged energy spectrum shown in figure 17 indicates that, for a 
$2048\times 2048$ grid, the energy decreases by 10 orders of magnitude (corresponding to an 
$O(10^{-5})$ relative error for the velocity field) between modes with 
$k=1$ and modes with 
$k=170$ (indicated by the dashed line in the figure), the latter corresponding to the dealiasing cutoff on a 
$512\times 512$ grid. This is consistent with the result shown in figure 1(a) and suggests that our simulations of turbulent flow are fairly well-resolved even on 
$512\times 512$ grids.
Time-averaged energy spectra of (a) the periodic orbits in the family presented in figure 12 computed on a 
$256\times 256$ grid and (b) turbulent trajectory computed on a 
$2048\times 2048$ grid. The dashed vertical line in (a) represents the hyperviscosity threshold 
$k_v$, and in (b) it represents the two-thirds dealiasing cutoff corresponding to a 
$512\times 512$ grid. The spectra at 
$k=1$ have been normalized to unity in both cases.
Solutions of the Euler equation were computed on a 
$256\times 256$ grid, which provides sufficient resolution for the smooth solutions we are interested in. Indeed, as seen in figure 17, for a family of time-periodic ECSs, we find the time-averaged energy to decrease by 10 orders of magnitude between modes with 
$k=1$ and modes with 
$k=k_v=40$ (indicated by the dashed line in the figure), suggesting that hyperviscosity has a negligible effect on the solution. Note that dealiasing cutoff corresponds to a far higher frequency 
$k=85$ in this case.
A.2. Newton-GMRES solver
Exact coherent structures of all types were computed using a Krylov subspace Newton method (see, e.g. Viswanath Reference Viswanath2007), which uses GMRES to solve the equation
\begin{equation} J(\omega_j) \Delta\omega_j=F(\omega_j) \end{equation}
for the vorticity field, where 
$F$ represents the relevant conditions for either equilibria or periodic orbits (specified below), and 
$J$ is the Jacobian of 
$F$. The initial condition 
$\omega _0$ is either taken to be a properly smoothed snapshot of the turbulent flow field or constructed using homotopy from previously converged solutions. Vorticity is then updated iteratively as 
$\omega _{j+1}=\omega _j-\Delta \omega _j$ until the relative error 
$\zeta =\|F(\omega )\|/\|\omega \|$ becomes less than 
$5\times 10^{-5}$. Note that, in terms of the velocity field, the accuracy of our converged solutions corresponds to an 
$O(10^{-5})$ relative residual, which is the same as the relative error of our numerical simulations. Hence, converged solutions can be considered numerically exact, with their accuracy limited mainly by the spatial and temporal resolution of the simulations.
Since ECSs of the Euler equation have marginal directions, 
$J$ becomes non-invertible near a solution of 
$F(\omega )=0$. Although, due to the presence of marginal directions, equation (A11) has infinitely many solutions in the full state space, its projection
\begin{equation} \hat{P}_KJ(\omega_j) \Delta\omega_j=\hat{P}_KF(\omega_j) \end{equation}onto the Krylov subspace
\begin{equation} K = {\rm span}\{Jx, J^2x, J^3 x,\ldots\}, \end{equation}
constructed using Arnoldi iteration (where vector 
$x$ does not lie in the kernel of 
$J$) has a well-defined unique solution so long as 
$\dim K\leq \dim \textrm {image}(J)$. In practice, the latter condition is always satisfied.
A.3. Equilibria
For equilibria, we used the function
\begin{equation} F(\omega,\boldsymbol{v}) = ((\boldsymbol{u-v})\boldsymbol{\cdot}\boldsymbol{\nabla}\omega,\,E - E_0,\, \phi_x ,\, \phi_y), \end{equation}
where the additional constraints fix the energy 
$E$ and the spatial position 
$(\phi _x,\phi _y)$ of the solution. The energy 
$E_0$ was typically set equal to that of the initial condition 
$\omega _0$. The velocity 
$\boldsymbol {v}$ of the reference frame was set to be zero for absolute equilibria and treated as part of our solution space for relative equilibria.
A.4. Periodic orbits
For time-periodic solutions, we used the objective function
\begin{equation} F(\omega) = \omega-\hat L\hat\varPhi_T\omega, \end{equation}
where 
$\hat \varPhi _T$ denotes the time-
$T$ map, i.e. 
$\hat {\varPhi }\omega (0)=\omega (T)$, and 
$\hat L$ is a linear operator that includes translations (to allow relative solutions) and discrete rotations (to allow preperiodic solutions such as those presented in figure 6). For time-periodic solutions, we fixed the period 
$T$ instead of the energy, since the former constraint was more straightforward to implement numerically.
The introduction of the hyperviscous term (A10) was crucial for converging time-periodic orbits, as it prevented Newton iterations from introducing high frequencies into the solution. As the solution become better converged, the threshold frequency 
$k_v$ can be increased to make sure the converged solution satisfies the Euler equation to a higher accuracy.
A.5. Homotopy
Solution families reported here have been constructed using homotopy. In its standard form, homotopy represents a linear interpolation
\begin{equation} \boldsymbol{u}(\sigma)=\sigma \boldsymbol{u}(0)+(1-\sigma)\boldsymbol{u}(1) \end{equation}
between a pair of states 
$\boldsymbol {u}(0)$ and 
$\boldsymbol {u}(1)$ characterized by a continuous parameter 
$\sigma \in [0,1]$. Note that the flow fields defined by (A16) generally do not satisfy the Euler equation even when both 
$\boldsymbol {u}(0)$ and 
$\boldsymbol {u}(1)$ do. Hence, one could use these flow fields as initial conditions to the Newton-GMRES solver. This naive approach works sufficiently well for solutions that are close to each other.
For solutions that are not particularly close, simple interpolation generates initial conditions that are not close to solutions either unless 
$\sigma$ is close to 0 or 1. In such cases, a sequential interpolation scheme can be used which employs a homotopy between previously converged solutions. Namely, for some discrete set 
$\{\sigma _j\}\in (0,1)$ and previously converged solutions 
$\boldsymbol {u}(\sigma _{j-1})$ and 
$\boldsymbol {u}(1)$ we define the initial condition for 
$\boldsymbol {u}(\sigma _{j-1})$ as
\begin{equation} \boldsymbol{u}^*(\sigma_j) = \frac{(1-\sigma_j)\boldsymbol{u}(\sigma_{j-1})+(\sigma_j-\sigma_{j-1})\boldsymbol{u}(1)}{1-\sigma_{j-1}}. \end{equation}
This method yields better (smoother) approximations to continuous families connecting distant states. Indeed, while the homotopy may yield a smooth family of initial conditions 
$\boldsymbol {u}^*(\sigma _j)$, converged states 
$\boldsymbol {u}(\sigma _j)$ may not form a smooth family since the difference 
$\boldsymbol {u}^*(\sigma _j)-\boldsymbol {u}(\sigma _j)$ is not guaranteed to be small unless the initial condition is close to a solution. Most of the families of solutions described in the paper were constructed via this method.
A.6. Smoothing of initial conditions
In order to speed up convergence of Newton-GMRES iterations and improve the success rate of the solver, initial conditions found from recurrence analysis of turbulent flow should be properly smoothed. The rate of success depends rather sensitively on the choice of both the smoothing method and its various parameters. While quantifying the rate of success and convergence speed as a function of all these choices is outside the scope of our study, we note that sufficiently aggressive smoothing yields an effectively 100 % success rate, as is the case for the solutions shown in figure 15. Below we describe three different smoothing methods that have been used to prepare initial conditions for the Newton-GMRES solver, with the results compared in figure 18 for an initial condition that corresponds to a nearly time-periodic large-scale flow.
(i) Spectral smoothing. The flow field can be smoothed by applying an axially symmetric mask in Fourier space. The most common and effective choice was the Gaussian blur which corresponds to the mask
(A18)with a characteristic wavenumber\begin{equation} M(\boldsymbol{k};k_s) =\exp\left(-\frac{k^2}{k_s^2}\right) \end{equation}$k_s$ controlling the degree of smoothing. Another mask that was occasionally used is
(A19)In a few instances we also used the Fourier filtering operator\begin{equation} M(\boldsymbol{k};k_s) = \min\left(1,\frac{k_s}{k} \right). \end{equation}$\hat {L}_{k_s}$ used to split the flow into the large- and small-scale component, which corresponds to the mask
(A20)where we set\begin{equation} M(\boldsymbol{k};k_s) = \frac12\left[1-\tanh\left(\frac{k-k_s}\ell\right)\right], \end{equation}$\ell =5$.Spectral smoothing is a very standard and simple approach which is effective at removing small-scale structure while preserving some large-scale structure of the flow. It can noticeably change the dynamics, however.
(ii) Hyperviscous smoothing. The fine structures in the flow field can also be effectively eliminated via the action of hyperviscosity over a suitable long period of time. For time-periodic states, we evolve the flow over one period using the time-integrator for the Euler equation with
$k_v$ set to a relatively low value (compared with that used in the Newton-GMRES solver). This approach allows small-scale structure to be removed without significantly affecting the spatial or temporal structure of the large-scale flow.(iii) Stream function smoothing. Solutions of the Euler equation can vary quickly transverse to, but not along, the streamlines of the flow. Hence, it is helpful to average the vorticity field
$\omega$ along the instantaneous streamlines. The smoothed vorticity field is then given by
(A21)Unlike the other two smoothing approaches, this method is better at preserving the edges of the vortex cores where vorticity changes rapidly, while damping out asymmetrical structures around the vortices. For instance, the flow shown in figures 18(a), 18(b) and 18(d), features vortices with pronounced azimuthal structure; this structure has been largely removed by the stream function smoothing as shown in figure 18(c). Hence, this smoothing method is unsuitable for computing solutions featuring structured vortices, such as the tripolar vortex seen in figure 11.\begin{equation} \bar\omega(x,y) = \iint \exp\left(\frac{[\psi(x,y)-\psi(x',y')]^2}{2\sigma^2}\right)\omega(x',y')\,{\rm d} x'\,{{\rm d}\kern0.05em y}'. \end{equation}
A snapshot of turbulent flow field (a) and initial conditions for the Newton-GMRES solver obtained by applying spectral smoothing with 
$k_s=7$ (b), stream function smoothing with 
$\sigma =0.5$ (c), and hyperviscous smoothing with 
$k_v=7$, 
$\alpha =1$, and 
$\beta =1/13$. The relative error 
$\zeta$ for the corresponding initial conditions is 
$0.226$ in (a), 
$0.0932$ in (b), 
$0.105$ in (c) and 
$0.163$ in (d).
In particular, in computing the solutions reported in figure 15, we created initial conditions by sequentially performing stream function smoothing followed by spectral smoothing and, ultimately, hyperviscous smoothing, with the combination of three methods reducing the relative error 
$\zeta$ from 
$O(0.1)$ to 
$O(10^{-3})$. A few tens of Newton iterations were then sufficient to achieve convergence with 
$\zeta <5\times 10^{-5}$ for all initial conditions. This aggressive smoothing is the reason why some properties of converged solutions differ substantially from those of the corresponding turbulent flow snapshots, especially at later times.
Open access
