2.1. Kolmogorov flow

We consider two-dimensional turbulence in a square domain with periodic boundary conditions, driven by a monochromatic body force in the streamwise direction. The out-of-plane vorticity $\omega = \partial _x v - \partial _y u$ , where the velocity $\boldsymbol{u} = (u,v)$ , evolves according to

(2.1) \begin{equation} \partial _t \omega + \boldsymbol{u} \boldsymbol{\cdot }\boldsymbol{\nabla }{\omega } = \frac {1}{\textit{Re}} {\nabla} ^2 \omega - n \cos {ny}. \end{equation}

In this non-dimensionalisation, we have chosen a length scale $1/k^{*}$ , which is the inverse of the fundamental wavenumber of the domain $k^* = 2\pi / L^*$ , and a time scale $1/\sqrt {k^*\chi ^*}$ , where $\chi ^*$ is the amplitude of the forcing in the momentum equation. These length and time scales lead to the definition of the Reynolds number in this flow ${\textit{Re}} := \sqrt {\chi ^*/k^{*3}}/\nu$ , where $\nu$ is the kinematic viscosity. Throughout this work, we set the forcing wavenumber $n = 4$ , as has been common in previous studies (Chandler & Kerswell Reference Chandler and Kerswell2013; Page, Brenner & Kerswell Reference Page, Brenner and Kerswell2021, Reference Page, Norgaard, Brenner and Kerswell2024b ).

Equation (2.1) is equivariant under continuous shifts in the streamwise ( $x$ ) direction, $\mathscr{T}_s : \omega (x,y) \to \omega (x+s, y)$ , under discrete shift-reflects by a half-wavelength in $y$ , $\mathscr{S}: \omega (x,y) \to -\omega (-x, y+\pi /4)$ , and under discrete rotations by $\pi$ , $\mathscr{R} : \omega (x,y) \to \omega (-x,-y)$ . The existence of the continuous symmetry means that we generically expect exact solutions to be relative periodic orbits that shift a finite distance in $x$ over one period. This is the case for the vast majority of solutions in our library.

Some key integral observables reported throughout this work are the total kinetic energy,

(2.2) \begin{equation} E(t) := \frac {1}{2} \langle \boldsymbol{u}^2 \rangle _V, \end{equation}

the total dissipation rate,

(2.3) \begin{equation} D := \frac {1}{{\textit{Re}}} \langle | \boldsymbol{\nabla }\boldsymbol{u} | ^2 \rangle _V = \frac {1}{{\textit{Re}}} \langle \omega ^2 \rangle _V, \end{equation}

and the total production rate,

(2.4) \begin{equation} I := \langle u \sin (ny) \rangle _V, \end{equation}

where the average over the volume $V$ is defined as

(2.5) \begin{equation} \langle \bullet \rangle _V := \frac {1}{(2\pi)^2} \iint \bullet \; \textrm{d}^2\boldsymbol{x}. \end{equation}

We solve (2.1) using the spectral version of the JAX-CFD solver (Kochkov et al. Reference Kochkov, Smith, Alieva, Wang, Brenner and Hoyer2021; Dresdner et al. Reference Dresdner, Kochkov, Norgaard, Zepeda-Núñez, Smith, Brenner and Hoyer2022). At each time step in the solver, the velocity field is computed by solving the Poisson equation ${\nabla} ^2 \psi = -\omega$ , where the streamfunction $\psi$ is related to the induced velocity components via $u = \partial _y \psi, v = -\partial _x \psi$ . The resolution is varied depending on the value of ${\textit{Re}}$ : $N_x\times N_y = 128 \times 128$ when ${\textit{Re}} \leqslant 150$ , rising to $256\times 256$ when $150 \lt {\textit{Re}} \leqslant 300$ and finally $512 \times 512$ when ${\textit{Re}} \gt 300$ .

Our analysis begins with a large set of RPOs converged at ${\textit{Re}}=40$ and ${\textit{Re}}=100$ , which are documented by Page et al. (Reference Page, Holey, Brenner and Kerswell2024a ) and Page et al. (Reference Page, Norgaard, Brenner and Kerswell2024b ) (see those papers for full details of the solutions, including periods, shifts, Floquet exponents etc.). Altogether, we begin with 174 RPOs at ${\textit{Re}} = 40$ and 151 RPOs at ${\textit{Re}} = 100$ . To investigate how the dynamical importance of these solutions varies with ${\textit{Re}}$ , we perform a continuation of the solutions at ${\textit{Re}} = 40$ and ${\textit{Re}} = 100$ to higher Reynolds number. The continuation results in a large set of solutions in the region ${\textit{Re}} \in (30, 1100)$ , with each RPO belonging to a solution branch initialised at either ${\textit{Re}} = 40$ or ${\textit{Re}} = 100$ . Branches can go through multiple fold bifurcations throughout the continuation, so that ${\textit{Re}}$ is not necessarily monotonic along the curves.

2.2. Vortex identification

Motivated by a possible connection of exact solutions in viscous Kolmogorov flow to solutions of the Euler equation as ${\textit{Re}} \to \infty$ (Kim & Okamoto Reference Kim and Okamoto2015; Kim et al. Reference Kim, Miyaji and Okamoto2017; Zhigunov & Grigoriev Reference Zhigunov and Grigoriev2023), part of this work explores the ability of point vortex solutions to capture the dynamics of the large-scale vortices in the Navier–Stokes RPOs. Note that while a Kolmogorov forcing profile is used throughout this work, previous work (Gallet & Young Reference Gallet and Young2013; Kim et al. Reference Kim, Miyaji and Okamoto2017; Zhigunov & Grigoriev Reference Zhigunov and Grigoriev2023) indicates that the large-scale structure of the flow may detach from the forcing in (2.1) at very high ${\textit{Re}}$ .

To model the vortex cores of the Kolmogorov RPOs with point vortices, we must first identify and extract coherent vortices from the turbulent flows. We follow the methodology of Page et al. (Reference Page, Holey, Brenner and Kerswell2024a ) and first compute the root-mean-square (r.m.s.) vorticity fluctuations $\omega _{rms} := \sqrt {\langle (\omega (\boldsymbol{x},t) - \overline {\omega }(y))^2 \rangle }_V$ , where the overline denotes an average over time and $x$ direction. Spatially localised vortices are then extracted as any connected regions $\mathcal{V}_i$ , where

(2.6) \begin{equation} | \omega (\boldsymbol{x},t) - \omega _{rms} | \geqslant 2 \omega _{rms}. \end{equation}

Some examples of these connected regions extracted from a vorticity snapshot are shown in figure 1 at both ${\textit{Re}}=100$ and ${\textit{Re}}=400$ . Note that this method requires that a reference vorticity magnitude must be chosen ( $2 \omega _{rms}$ in (2.6)), which influences the boundaries of the extracted vortex cores. This method is agnostic to the shape of the connected region, so that strong vorticity filaments are also occasionally extracted as vortex ‘cores’ – note the thin structure stretching across the periodic boundary in figure 1(b).

Figure 1. Application of the vortex extraction criterion discussed in the text to two example snapshots at (a,b) ${\textit{Re}}=100$ and (c,d) ${\textit{Re}}=400$ . Contours are of the out-of-plane vorticity, white lines identify regions identified by the threshold (2.6), purple text indicates the region area and red lines (which here enclose areas so small they appear as points in the visualisation) highlight which of these regions were discarded due to the bound imposed on the region area.

Following vortex boundary identification, the circulation $\varGamma _i := \int _{\boldsymbol{x} \in \mathcal{V}_i} \omega (\boldsymbol{x})\, \textrm{d}^2\boldsymbol{x}$ , vortex area $A_i := \int _{\boldsymbol{x} \in \mathcal{V}_i} \mathbb{I}(\boldsymbol{x}) \,\textrm{d}^2\boldsymbol{x}$ and the centre of vorticity $\tilde {\boldsymbol{x}}_{i} = (\tilde {x}_i, \tilde {y}_i) = ({1}/{\varGamma _i}) \int _{\boldsymbol{x} \in \mathcal{V}_i} \omega (\boldsymbol{x}) ( x \,\textrm{d}x, y \,\textrm{d}y)$ of each vortex core can all be computed. Similar to the approach by Benzi, Patarnello & Santangelo (Reference Benzi, Patarnello and Santangelo1987), we consider only ‘large-scale’ vortices, ignoring any structures with area less than 10 % of the largest vortex $A_i \lt 0.1 \max _j(A_{\!j})$ (see discarded vortices in red in figure 1). Finally, to satisfy the Gauss constraint, $\int \omega (\boldsymbol{x})\,\textrm{d}^2\boldsymbol{x} = 0$ , we must ensure that $\sum _j \varGamma _j = 0$ , which is enforced by equally modifying the extracted circulation of each vortex core, $\varGamma _j \to \varGamma _j - \langle \varGamma \rangle$ , prior to initialising a point vortex computation.

Our attempt to fit point vortex solutions is performed at ${\textit{Re}}=100$ . At this value, we find an average $\langle |\varGamma |\rangle = 2.84$ with standard deviation $\sigma _{\varGamma } = 2.68$ for the extracted cores in the RPO library (the distribution has a long right tail, see Page et al. Reference Page, Holey, Brenner and Kerswell2024a ). As a result, the induced velocities are $O(1)$ . This should be contrasted to the much weaker time average velocity, which is approximately an order of magnitude smaller (statistics are presented in § 3.5).

2.3. Doubly periodic point vortex model

We now outline the formulation for the point vortex problem in a doubly periodic domain. The following equations of motion in this system were first derived by Weiss & McWilliams (Reference Weiss and McWilliams1991). An alternative derivation is presented here following the approach adopted by Aref et al. (Reference Aref, Newton, Stremler, Tokieda and Vainchtein2003) for a periodic strip. We begin by considering $N_v$ point vortices in an infinite two-dimensional domain without a background velocity, where the $j\textrm{th}$ vortex is located at the complex position $z_{j} = x_{j} + iy_{j}$ . Each point vortex moves under the induced velocity from the $N_v-1$ other vortices in the domain. The equations of motion for these point vortices are

(2.7) \begin{equation} \dot {\overline {z}}_j = \dot {x}_j - i\dot {y}_j = \frac {1}{2\pi i} \sum _{k = 1}^{N_v} \phantom {}^{\prime } \frac {\varGamma _k}{ z_j - z_k}, \end{equation}

where the overbar represents the complex conjugate, the dot denotes a time derivative, $\varGamma _k$ is the circulation of the $k\textrm{th}$ point vortex and the prime indicates the omission of any singular terms (i.e. no self-induced velocity from $j = k$ ). These equations form a Hamiltonian system (e.g. see Batchelor Reference Batchelor1967),

(2.8) \begin{equation} \varGamma _j \dot {\overline {z}}_j = -\frac {1}{i}\frac {\partial H_u}{\partial z_j}, \end{equation}

where the unbounded Hamiltonian is

(2.9) \begin{equation} H_u = -\frac {1}{4\pi }\sum _{j=1}^{N_v}\sum _{k=1}^{N_v}\phantom {}^{\prime } \varGamma _{j} \varGamma _{k} \log {|z_{j} - z_{k}|}. \end{equation}

This result can be generalised to the bounded, doubly periodic $L \times L$ domain by introducing ghost vortices at $z_j + nL + imL$ with $m,n \in \mathbb{Z} \backslash \{ 0 \}, j \in \{1, \ldots, N_v \}$ . The equations of motion then become

(2.10) \begin{equation} \dot {\overline {z}}_j = \frac {1}{2\pi i} \sum _{n,m = -\infty }^{\infty } \sum _{k = 1}^{N_v} \phantom {}^{\prime } \frac {\varGamma _k}{z_j - z_k - nL - imL}. \end{equation}

We can make use of the identity

(2.11) \begin{equation} \sum _{n=-\infty }^{\infty } \frac {1}{x-n} = \pi \cot (\pi x) \end{equation}

to simplify (2.10) by evaluating either the summation over $n$

(2.12) \begin{equation} \dot {\overline {z}}_j = \frac {1}{2 Li} \sum _{k} \varGamma _k \sum _{m=-\infty }^{\infty } \cot \left (\pi \left (\frac {z_j-z_k}{L} - im\right)\right)\!, \end{equation}

or over $m$

(2.13) \begin{equation} \dot {\overline {z}}_j = \frac {1}{2L} \sum _{k} \varGamma _k \sum _{n=-\infty }^{\infty } \cot \left (\pi i\left (\frac {z_j - z_k}{L} - n\right)\right)\!. \end{equation}

Equations (2.12) and (2.13) correspond physically to an infinite number of strips, periodic in the horizontal and vertical directions, respectively (Aref et al. Reference Aref, Newton, Stremler, Tokieda and Vainchtein2003). We can then use the trigonometric identity

(2.14) \begin{equation} \cot (a + ib) = \frac {\sin (2a) - i\sinh (2b)}{\cosh (2b) - \cos (2a)} \end{equation}

to extract $\dot {x}_j$ from (2.13) and $\dot {y}_j$ from (2.12). This results in the final equations of motion for the doubly periodic domain

(2.15) \begin{equation} \begin{pmatrix} \dot {x}_j \\ \dot {y}_j \end{pmatrix} = \frac {1}{2L} \sum _{k=1}^{N_v} \phantom {}^{\prime } \varGamma _k \sum _{n=-\infty }^{\infty } \begin{pmatrix} -S_n\left (\frac {2\pi }{L}(y_j - y_k), \frac {2\pi }{L}(x_j - x_k)\right) \\ S_n\left (\frac {2\pi }{L}(x_j - x_k), \frac {2\pi }{L}(y_j - y_k)\right) \end{pmatrix}\!, \end{equation}

where

(2.16) \begin{equation} S_n(a,b) = \frac {\sin ( a)}{\cosh (b - 2\pi n) - \cos (a)}. \end{equation}

The reason for extracting $\dot {x}_j$ from (2.13) rather than (2.12) is due to the need to numerically truncate the summations. Consider, for example, the equilibrium configuration of $N_v = 4$ point vortices with equal circulation magnitude, equispaced in the periodic domain. Truncating the summation in (2.13) at finite $n$ results in two unmatched ghost vortices in the final vertically periodic strip in the summation, which causes a very small, but non-zero induced velocity in both the $x$ and $y$ directions. However, this induced numerical ‘drift’ in $x$ is substantially weaker than in $y$ . The same argument can be made for extracting $\dot {y}_j$ from (2.12), where the number of horizontal strips is truncated.

These equations of motion with $L = 2\pi$ match those derived via the method of Laplace transforms by Weiss & McWilliams (Reference Weiss and McWilliams1991) (up to a rescaling of time to match their conventions). Weiss & McWilliams (Reference Weiss and McWilliams1991) also derived the Hamiltonian for the doubly periodic system of point vortices

(2.17) \begin{equation} H = -\frac {1}{4\pi } \sum _{j,k = 1}^{N_v} \phantom {}^{\prime } \frac {\varGamma _j\varGamma _k}{2} h\left ( \frac {2\pi }{L}(x_j - x_k), \frac {2\pi }{L}(y_j - y_k) \right)\!, \end{equation}

where

(2.18) \begin{equation} h(a,b) = \sum _{m=-\infty }^{\infty } \log {\left ( \frac { \cosh (a - 2\pi m) - \cos (b) }{ \cosh (2\pi m) } \right)} - \frac {a^2}{2\pi }. \end{equation}

This Hamiltonian is clearly even in both $x$ and $y$ . It is clearly periodic in $y$ – periodicity in $x$ can be shown by manipulation of the sum truncated at large $M$ , before taking the $M\to \infty$ limit. The periodic boundary conditions break the rotational symmetry of the unbounded domain, leaving two additional constants of motion,

(2.19) \begin{align} P_x = \sum _{j=1}^{N_v} \varGamma _j x_j, \quad P_y = \sum _{j=1}^{N_v} \varGamma _j y_j, \end{align}

arising from the translational symmetry in $x$ and $y$ , respectively.

We solve (2.15) using a differentiable point vortex solver JAX-PV, developed in our previous work (Cleary & Page Reference Cleary and Page2023). This numerical solver is built on the JAX library (Bradbury et al. Reference Bradbury2018), to allow for efficient computation of the gradient of the time-forward map $\boldsymbol{f}^t(\textbf{x})$ of (2.15), where $\textbf{x} := (x_1, y_1, \ldots, x_{N_v}, y_{N_v})$ , with respect to $\mathbf{x}$ . The JAX framework also brings efficiency benefits such as just-in-time compilation and auto-vectorisation. Time integration is performed with the symplectic second order Runge–Kutta scheme at a fixed time step of $\delta t = 10^{-3}$ .

3.1. Continuation in Re

In this section, an arclength continuation effort of the library of RPOs in two-dimensional, doubly periodic Kolmogorov flow is presented. As discussed in § 2.1, the starting library of RPOs for the continuation comprises 174 solutions at ${\textit{Re}} = 40$ and 151 solutions at ${\textit{Re}} = 100$ . Each of these solutions is defined by the state vector $\boldsymbol{X} = [\omega, s, T]^T$ , where $\omega$ is a vorticity snapshot along the RPO, $s$ is the required translational shift and $T$ is the period of the RPO. To begin the continuation of a branch, a starting solution $\boldsymbol{X}({\textit{Re}}_0)$ at ${\textit{Re}}_0$ is naively perturbed in ${\textit{Re}}$ by $\delta {\textit{Re}}$ , ${\textit{Re}}_1 = {\textit{Re}}_0 + \delta {\textit{Re}}$ . Using $\boldsymbol{X}({\textit{Re}}_0)$ as the initialisation at ${\textit{Re}}_1$ , this solution is converged via a standard Newton-GMRES-Hookstep solver to yield $\boldsymbol{X}({\textit{Re}}_1)$ . This is robust as long as $\delta {\textit{Re}}$ is small, but cannot work when there are turning points in the solution branch. For this reason, the continuation proceeds via arclength continuation after the first successful convergence. The branch is then instead parametrised by its arclength $r$ , which increases monotonically along the branch. The state vector is extended to include ${\textit{Re}}$ , i.e. $\boldsymbol{X}(r) = [\omega, s, T, {\textit{Re}}]^T$ , along with an additional constraint

(3.1) \begin{equation} \frac {\partial \boldsymbol{X}}{\partial r} \boldsymbol{\cdot }\frac {\partial \boldsymbol{X}}{\partial r} = 1, \end{equation}

to match the additional unknown. The next step size in arclength $\delta r$ is controlled by

(3.2) \begin{equation} \delta r_{i-1} = \sqrt {(\boldsymbol{X}(r_{i-1}) - \boldsymbol{X}(r_{i-2}))^2}, \end{equation}

so that $\boldsymbol{X}(r_i = r_{i-1} + \delta r_{i-1})$ is computed via the modified Newton-GMRES-Hookstep solver with the extra constraint (3.1). The initialisation for $\boldsymbol{X}(r_i)$ was set by linearly extrapolating along $r$ from the previous two states along the branch $\boldsymbol{X}(r_{i-2})$ and $\boldsymbol{X}(r_{i-i})$ . The methodology matches that presented by Chandler & Kerswell (Reference Chandler and Kerswell2013), and full details can be found in the appendix of that paper.

The continuation algorithm was automated in the following ways. For each branch, the initial perturbation in ${\textit{Re}}$ was set to $\delta {\textit{Re}} = 2.5$ . If the perturbed solution was not converged by the standard Newton-GMRES-Hookstep solver, $\delta {\textit{Re}}$ was halved repeatedly until a convergence was found, up to a maximum of four halvings. After this first successful step in ${\textit{Re}}$ , the step size in arclength $\delta r$ was set according to (3.2), and $\boldsymbol{X}(r_i)$ was initialised using the two existing states along the branch. If the initialised state was not converged by the modified Newton-GMRES-Hookstep solver, then $\delta r$ was halved repeatedly until a convergence was found, again up to a maximum of four halvings. The number of halvings of $\delta r$ required for convergence was recorded along the branch, as a proxy for the ease of convergence along the branch. If both previous solutions along the branch did not require any halvings, then $\delta r$ was doubled, to speed up the continuation. Convergence of a solution was deemed to have failed if either (1) more than 100 Newtons iterations were required or (2) more than 10 Hooksteps were required more than three times. Continuation of a branch was also terminated if either (1) more than 50 solutions were converged along that branch, a restriction added due to many curves repeatedly folding, or (2) $\delta r$ was halved more than four times for the current solution. The starting ${\textit{Re}} = 40$ RPOs were continued upwards in ${\textit{Re}}$ , while the starting ${\textit{Re}} = 100$ RPOs were continued in separate computations both upwards and downwards in ${\textit{Re}}$ .

Long trajectories of Kolmogorov flow were also simulated at various ${\textit{Re}} \in [30, 1000]$ to compute reference turbulence statistics. We use the probability density function (p.d.f.) of the dissipation in these long computations as a simple measure of the dynamical significance of the RPOs to the turbulence, labelling a solution as ‘dynamically relevant’ if its mean dissipation rate $D$ lies between the $1\textrm{st}$ and $99\textrm{th}$ percentile of the reference p.d.f. This method is chosen for its simplicity and ease of evaluation, but will lead to some mislabelling. Better observables should be considered in future work (see e.g. the strategies of Krygier et al. (Reference Krygier, Pughe-Sanford and Grigoriev2021) and Page et al. (Reference Page, Norgaard, Brenner and Kerswell2024b )). Each long trajectory was simulated for $10^5$ advective time units, and sampled at intervals of 1 advective time unit. Contour plots of the p.d.f.s of the dissipation of these trajectories are shown in various figures throughout the paper to aid identification of the turbulent attractor, with boundaries of the dynamically important region of the turbulent attractor indicated by dashed black lines.

The time-averaged dissipation rate for the RPOs along each branch was computed and is plotted in figure 2 over the background turbulent p.d.f.s. This quantity was chosen to reveal the scaling of the volume-averaged velocity gradients with ${\textit{Re}}$ for each branch. The scaling of the time-and-volume-averaged dissipation, scaled by ${\textit{Re}}$ , $\overline {{\textit{Re}} \, D}$ , of the fully turbulent flow is also shown in the bottom panel of figure 2 and matches the ‘asymptotic’ scaling of ${\textit{Re}}^{1/2}$ identified by Chandler & Kerswell (Reference Chandler and Kerswell2013) (their observations held for ${\textit{Re}}\leqslant 200$ ). The majority of branches seeded from RPOs at ${\textit{Re}} = 100$ do not leave the region ${\textit{Re}} \in [60, 200]$ , with the branches repeatedly turning (this was the reason for terminating many continuations at 50 steps).

Figure 2. Continuation of periodic orbits and their overlap with the turbulent dissipation p.d.f. (Top) Histograms of the number of monotonic-in- ${\textit{Re}}$ subsections $N_{upo}$ of all solution branches (grey) and monotonic-in- ${\textit{Re}}$ subsections fully within the dynamically important region (red) as a function of ${\textit{Re}}$ . Both histograms are computed using 50 bins, spaced logarithmically over the range ${\textit{Re}} \in [20, 1000]$ . (Middle) The time-averaged dissipation rate ${\textit{Re}} \, D$ against ${\textit{Re}}$ of the arclength continuation of the initial library of RPOs starting at both ${\textit{Re}} = 40$ and ${\textit{Re}} = 100$ . The contour plot shows the reference p.d.f.s of dissipation rate, with the $1\textrm{st}$ and $99\textrm{th}$ percentiles indicated by the dashed black lines. The red/blue/green branches indicate those with terminal RPOs above/within/below this dynamically important region. Circles denote the terminal solution along each branch. Filled circles denote the branches which were terminated as 50 states were converged, while empty circles denote the branches which could not be continued further. (Bottom) The time-averaged, scaled dissipation rate $\overline { {\textit{Re}} \, D }$ of a long-time simulation is shown as a function of ${\textit{Re}}$ (black), as well as the scaling law ${\textit{Re}}^{1/2}$ (blue).

In the top panel of figure 2, the number of monotonic-in- ${\textit{Re}}$ subsections of solution branches is tracked – this is a proxy for the number of unique RPOs in a particular ${\textit{Re}}$ interval. The number of dynamically important solutions decreases more sharply with ${\textit{Re}}$ than the total number of monotonic subsections. This departure of the RPOs from ‘dynamically relevant’ dissipation values could be attributed to a number of causes. One is the (unlikely) exit from the attractor of RPOs via boundary crises. Alternatively, the RPOs may remain in the attractor but with an increasingly rare probability of visits from a turbulent orbit, hence, it is deemed dynamically irrelevant under our scalar-observable based criterion. Finally, there is the possibility that the solutions were not in the attractor when they converged, but close by. This is not something that is obvious from low-dimensional projections which are typically used to assess the closeness of an RPO to the turbulence, but does seem to be borne out by the following analysis.

To facilitate a discussion about the ${\textit{Re}} \to \infty$ behaviour of the solutions, the branches are labelled as belonging to one of three classes according to the dynamical importance of the final RPO on each branch, as assessed by comparison to the turbulent PDF described previously. The final RPO of branches in class 1 have a larger dissipation rate than the dynamically important region, class 2 solutions lie in the dynamically important region (many remain confined to a small range of ${\textit{Re}}$ ) and class 3 solutions have a smaller dissipation rate than the dynamically important region. These three classes are coloured in red, blue and green, respectively, in figure 2. At high ${\textit{Re}}$ , the average dissipation rate of RPOs in class 2 scales similarly to the turbulent attractor, while those in class 1 and class 3 scale more strongly (most with a scaling $D\sim {\textit{Re}}$ ) or weakly ( $D \sim 1/{\textit{Re}}$ – though we are lacking large numbers of solutions in this regime), respectively. Note that the vast majority of solutions in classes 1 and 3 appear to continue monotonically with ${\textit{Re}}$ after they leave the dynamically relevant region.

3.2. Class 1: connections to exact unforced Euler solutions

Some example branches from solution class 1 are shown in more detail in figure 3. The RPOs in this class have stronger dissipation rates than the turbulence, or equivalently larger enstrophy $\langle \omega ^2 \rangle _V$ . As a result, class 1 solutions have very defined vortical structures – e.g. see the horizontal bars and clear vortex cores in figure 3. The strength of the local vorticity associated with these features increases dramatically as ${\textit{Re}}$ is raised. For example, the maximum vorticity nearly triples for the last two points highlighted in the second row of figure 3, while ${\textit{Re}}$ increases from approximately 100 to 220.

Figure 3. Time-averaged, scaled dissipation rate ${\textit{Re}} \, D$ against ${\textit{Re}}$ for three representative branches in class 1. A single vorticity snapshot for each highlighted RPO (circled numbers 1–3) at points indicated by crosses on the branch are shown in a row on the right of the figure, with arclength increasing from left to right. Note the different colour bars for each panel, highlighting the strengthening vortical structures.

The increasingly strong vorticity fields in class 1 meant that the baseline resolution was often inadequate as the continuation proceeded upwards in ${\textit{Re}}$ . As a result, some of the initial computations terminated due to a failure to converge before reaching 50 states. These branches are identified with filled diamonds in figure 3, and were continued further at this point by doubling the number of Fourier modes in both spatial directions. For solutions which began at ${\textit{Re}}=40$ , it was sometimes necessary to double the resolution again (to a maximum of $512\times 512$ ) to follow the branch further upwards in ${\textit{Re}}$ .

The dissipation scaling $D \sim {\textit{Re}}$ is consistent with an enstrophy $\langle \omega ^2 \rangle _V \sim {\textit{Re}}^2$ , which implies $\omega = O({\textit{Re}})$ if the vortex cores do not shrink. This scaling is confirmed for a number of branches in class 1 in figure 4, and the velocity field (largely induced by the high amplitude vortex cores) also follows this scaling. In addition, the ${\textit{Re}}$ -dependence of the period of the RPOs is also shown in figure 4, and scales like $T_{\textit{RPO}} \propto 1/{\textit{Re}}$ . The horizontal shift associated with the class 1 RPOs (figure 4 d) does not exhibit a consistent scaling with ${\textit{Re}}$ – though the majority solutions appear to exhibit $s_{\textit{RPO}} \sim \text{constant}$ as ${\textit{Re}}$ increases. These scalings motivate the introduction of new variables $\varOmega := \omega /{\textit{Re}}$ , $\boldsymbol U := \boldsymbol u /{\textit{Re}}$ and $\tau := {\textit{Re}} \, t$ , from which (2.1) becomes

(3.3) \begin{equation} \partial _{\tau }\varOmega + \boldsymbol U \boldsymbol{\cdot }\boldsymbol{\nabla }\varOmega = \frac {1}{{\textit{Re}}^2}\left (\Delta \varOmega - n\cos ny\right)\!. \end{equation}

An asymptotic solution would then be sought in the form of a regular perturbation series:

(3.4) \begin{align} \varOmega (\boldsymbol x, \tau) &= \varOmega _0(\boldsymbol x, \tau) + \frac {1}{{\textit{Re}}^2}\varOmega _1(\boldsymbol x, \tau) + {\cdots}, \end{align}

(3.5) \begin{align} \boldsymbol U(\boldsymbol x, \tau) &= \boldsymbol U_0(\boldsymbol x, \tau) + \frac {1}{{\textit{Re}}^2}\boldsymbol U_1(\boldsymbol x, \tau) + {\cdots} . \end{align}

Figure 4. Scalings with ${\textit{Re}}$ of (a) the periods $T_{\textit{RPO}}$ , (b) maximum vorticity values $\max |\omega |$ , (c) maximum velocity values $\max |\boldsymbol u|$ and (d) shift for a subset of branches in class 1.

At each order, $\Delta \varPsi _i = -\varOmega _i$ , where $\varPsi _i$ is the streamfunction associated with $\boldsymbol U_i$ . The dynamics at leading order are therefore governed by the inviscid vorticity equation (the curl of the Euler equation)

(3.6) \begin{equation} \partial _{\tau } \varOmega _0 + \boldsymbol U_0\boldsymbol{\cdot }\boldsymbol{\nabla }\varOmega _0 = 0, \end{equation}

while the first-order correction is

(3.7) \begin{equation} \partial _{\tau } \varOmega _1 + \boldsymbol U_0\boldsymbol{\cdot }\boldsymbol{\nabla }\varOmega _1 + \boldsymbol U_1\boldsymbol{\cdot }\boldsymbol{\nabla }\varOmega _0 = \Delta \varOmega _0 - n\cos ny. \end{equation}

The existence of a solution to the first-order problem requires that a series of solvability conditions are satisfied,

(3.8) \begin{equation} \int _0^T \langle \zeta _j \left (\Delta \varOmega _0 - n\cos ny \right)\rangle _V \,\textrm{d}\tau = 0, \end{equation}

where – as shown by Zhigunov & Grigoriev (Reference Zhigunov and Grigoriev2023) – the $\{\zeta _j\}$ are eigenfunctions associated with continuous symmetries of the $T-$ (relative) periodic Euler solution $\varOmega _0(\boldsymbol x, \tau)$ . This is a generalisation of the solvability condition of Okamoto (Reference Okamoto1994) for steady solutions, which is identical to the physical constraint derived by Batchelor (Reference Batchelor1956).

The analysis here suggests a direct connection of the finite- ${\textit{Re}}$ solutions to solutions of the unforced Euler equation, where it is well known that solutions exist as continuous families due to the conservation of Casimir functions (Stuart Reference Stuart1967; Zhigunov & Grigoriev Reference Zhigunov and Grigoriev2023). However, we expect that each viscous solution should connect to a single ‘member’ of an Euler family, which is selected through the viscous solvability conditions (3.8) (Okamoto Reference Okamoto1994) – no further bifurcations associated with other ‘members’ of a solution family at finite ${\textit{Re}}$ are necessary as there is no direct connection.

The scalings of class 1 solutions with ${\textit{Re}}$ reported in figure 4 confirm the expected asymptotic behaviour. These results indicate that many Kolmogorov flow solutions found as low as ${\textit{Re}}=40$ – which have a structure at low ${\textit{Re}}$ that depends intimately on the forcing – connect directly to solutions of the unforced Euler equation. The various ‘turbulent’ RPOs found here complement the ‘unimodal’ states of Kim & Okamoto (Reference Kim and Okamoto2015), which also appear to connect to solutions of the unforced Euler equation as ${\textit{Re}}\to \infty$ . Those authors used a framework in which the body force was assumed to be $O(1/{\textit{Re}})$ in the momentum equation, though this is equivalent to the $O(1)$ force considered here after the rescaling used to produce (3.3); see also the discussion of Kim & Okamoto (Reference Kim and Okamoto2010). Kim & Okamoto (Reference Kim and Okamoto2015) observed unimodal Euler solutions to emerge for $500 \lesssim {\textit{Re}}_{\textit{KO}} \lesssim 100\,000$ , where their Reynolds number, ${\textit{Re}}_{\textit{KO}}$ , is related to ours by ${\textit{Re}} \approx \sqrt {{\textit{Re}}_{\textit{KO}}}$ (approximate because of differences in the forcing profile). Their vorticity profiles, which have $\omega _{\textit{KO}} = O(1)$ as ${\textit{Re}}_{\textit{KO}}\to \infty$ , would also scale like $\omega = O({\textit{Re}})$ using our conventions. The ${\textit{Re}}$ values where unimodal states are found by Kim & Okamoto (Reference Kim and Okamoto2015) correspond in our problem to a range $22 \lesssim {\textit{Re}} \lesssim 300$ , which agrees well with the ${\textit{Re}}$ values for which the class 1 scaling becomes apparent.

Furthermore, some of our class 1 solutions share qualitative similarities with the Euler solutions reported by Zhigunov & Grigoriev (Reference Zhigunov and Grigoriev2023). Of the 184 solutions labelled as ‘class 1’, 18 consist of a pair of opposite-signed vortices – similar profiles were reported by those authors (e.g. see example (3) in figure 3). Ratios of the time-averaged enstrophy and energy for these class 1 RPOs are typically found to be between approximately 5 and 10 (not shown), somewhat higher than those reported by Zhigunov & Grigoriev (Reference Zhigunov and Grigoriev2023), which were between 2 and 4, while the appropriately scaled periods for our RPOs are also much larger, $T\, {\textit{Re}} = O(100)$ . However, it is still possible that some of our states are connected to those shown by Zhigunov & Grigoriev (Reference Zhigunov and Grigoriev2023), which could presumably be checked by converging the leading (Euler) component of the solution and performing homotopy through the continuous solution families.

The class 1 solutions discussed here, some of which were highlighted in figure 3, indicate that a direct viscous-to-inviscid connection is possible for a wide variety of flow states beyond the unimodal states discussed here (Kim & Okamoto Reference Kim and Okamoto2015; Kim et al. Reference Kim, Miyaji and Okamoto2017; Zhigunov & Grigoriev Reference Zhigunov and Grigoriev2023). The scaling of this solution class relative to the turbulence would seem to suggest that these states, while ‘dynamically relevant’ by our definition at the ${\textit{Re}}$ values where they were converged and potentially similar to turbulent trajectories, were outside of the turbulent attractor.

3.3. Classes 2 and 3

The RPOs in classes 2 and 3 tend to be dominated by two oppositely signed vortices, arranged similarly to the vortex condensates typically seen in high- ${\textit{Re}}$ and Euler solutions (Zhigunov & Grigoriev Reference Zhigunov and Grigoriev2023; Page et al. Reference Page, Holey, Brenner and Kerswell2024a ). The large-scale opposite-signed vortex pair in class 3 branches is less distinct from the background vorticity than in class 2 branches – class 3 solutions still appear to show evidence of the $n=4$ forcing wave. Multipolar structures (Carton et al. Reference Carton, Flierl and Polvani1989; Morel & Carton Reference Morel and Carton1994) are observed in the class 2 solutions. For example, note the co-rotating same-signed vortex ‘bound state’ in the top row of snapshots in figure 5, and the ‘tripole’ structures (central core with two opposite-signed satellites) in the second and third rows of figure 5. In the middle row (solution ‘2’) in figure 5, the tripole structure weakens with increasing ${\textit{Re}}$ , and the vorticity field becomes dominated by the opposite-signed large-scale vortices. In contrast to the class 2 solutions, the vorticity strength of solutions in class 3 remains practically constant, even as ${\textit{Re}}$ is increased by an order of magnitude from ${\textit{Re}} = 100$ to ${\textit{Re}} \approx 1000$ . The dominant structures visible in the solutions in class 3 also vary remarkably little with these large changes in ${\textit{Re}}$ .

Figure 5. Time-averaged, scaled dissipation rate ${\textit{Re}} \, D$ against ${\textit{Re}}$ for three representative branches in class 2. Three RPOs along each branch are sampled (indicated by crosses) and visualised on the right of the figure. A single vorticity snapshot for each highlighted RPO (circled numbers 1–3) at points indicated by crosses on the branch are shown in a row on the right of the figure, with arclength increasing from left to right.

Figure 6. Time-averaged, scaled dissipation rate ${\textit{Re}} \, D$ against ${\textit{Re}}$ for three representative branches in class 3. Three RPOs along each branch are sampled (indicated by crosses) and visualised on the right of the figure. A single vorticity snapshot for each highlighted RPO (circled numbers 1–3) at points indicated by crosses on the branch are shown in a row on the right of the figure, with arclength increasing from left to right.

Figure 7. Full range of production $I$ at each ${\textit{Re}}$ for the three furthest continued class 3 solutions. A logarithmic colourmap is used, and the numbering corresponds to the numbered branches in figure 6.

The apparent asymptotic behaviour of the class 3 solutions, with $D\sim 1/{\textit{Re}}$ , is consistent with vorticity fields $\omega = O(1)$ (as shown in figure 6). We also observe that $T_{\textit{RPO}} \approx \text{constant}$ as ${\textit{Re}}$ increases (not shown). This behaviour is consistent with a connection to solutions of a forced Euler equation. To see this, consider a regular perturbation expansion $(\omega, \boldsymbol u) = (\omega _0, \boldsymbol u_0) + {\textit{Re}}^{-1}(\omega _1, \boldsymbol u_1) + {\cdots}$ , from which the leading order contribution satisfies

(3.9) \begin{equation} \partial _t \omega _0 + \boldsymbol u_0 \boldsymbol{\cdot }\boldsymbol{\nabla } \omega _0 = -n\cos ny. \end{equation}

There are no dissipation mechanisms, hence, for $T$ -periodic solutions, the forcing term must extract the same amount of kinetic energy it injects over a cycle,

(3.10) \begin{equation} \int _0^T I\, \textrm{d}t = \int _0^T\langle u_0 \sin ny\rangle _V\, \textrm{d}t = 0. \end{equation}

The production over a complete cycle of three class 3 RPOs is reported in figure 7 along the solution branches. For two of the branches (labelled ‘1’ and ‘3’), negative values of the production are realised at higher ${\textit{Re}}$ , which is consistent with a forced-Euler connection in the limit ${\textit{Re}}\to \infty$ . Further continuation would be necessary to verify this trend.

3.4. New RPOs at ${\textit{Re}} = 100$

Due to the non-monotonicity of the solution branches, the original library at ${\textit{Re}}= 100$ can be supplemented with new RPOs by slicing branches anywhere they cross ${\textit{Re}} = 100$ in the arclength continuation. As a result, 101 new RPOs at ${\textit{Re}} = 100$ were converged. The initialisations in the Newton convergence for these new RPOs were generated by linearly interpolating between the two RPOs which bracket the point where the branch crosses ${\textit{Re}} = 100$ . Some branches crossed ${\textit{Re}} = 100$ more than once to yield more than one new unique RPO.

These new RPOs are reported in figure 8 in terms of their production, $I$ , and dissipation $D$ , normalised by the laminar value. Due to the localised nature in phase space of these RPOs, full coverage of the turbulent p.d.f. at this ${\textit{Re}}$ is difficult, and requires a large number of solutions. These new RPOs go some way to increasing the coverage of the turbulent p.d.f., though gaps remain at higher $D$ and the tails of the $I$ marginal p.d.f. are not covered.

Figure 8. Energy dissipation rate $D$ against energy production rate $I$ at ${\textit{Re}} = 100$ , both normalised by the laminar dissipation $D_{lam} = {\textit{Re}} / (2n^2)$ . The 101 new RPOs at ${\textit{Re}} = 100$ from the continuation are shown in red, the starting library of 151 RPOs are shown in blue. The grey background is the p.d.f. computed from a trajectory of $10^5$ samples, separated by 1 advective time unit. The contour levels of the p.d.f. are spaced logarithmically. (a) New very high dissipation RPOs, far from the turbulent attractor, which resulted from the continuation of the very high dissipation RPOs at ${\textit{Re}} = 40$ . (b) Zoom-in on the turbulent attractor.

3.5. Reconstruction of turbulent statistics

A variety of statistics, including full dissipation and energy p.d.f.s, were recently successfully reconstructed via an expansion in the statistics of individual RPOs at ${\textit{Re}} = 40$ by Page et al. (Reference Page, Norgaard, Brenner and Kerswell2024b ). In that paper, the weights in the expansion were determined from the invariant measure of a Markov chain, where the ‘states’ on the chain were the RPOs, and the current state on a turbulent trajectory was assigned by measuring the distance to each RPO in the latent space of an autoencoder and identifying the closest solution. We adopt a similar approach here across ${\textit{Re}} \in [40, 300]$ , albeit without building Markov chain models. Instead, RPO weights are determined by fitting to a single representative statistic, and their robustness is then assessed by reconstructing other statistics using the same weights (see the approach of Redfern et al. Reference Redfern, Lazer and Lucas2024). We generate the various libraries of RPOs at each target ${\textit{Re}}$ by slicing the solution branches from the arclength continuation, as described in § 3.4. This allows us to track how the weights of each RPO change as ${\textit{Re}}$ is increased and assess the utility of RPOs converged at one ${\textit{Re}}$ value in making statistical predictions at higher values.

The number of unique solutions converged at each ${\textit{Re}}$ considered is shown in figure 9. There is a reduction in the number of solutions reported here compared with the number of monotonic-in- ${\textit{Re}}$ subsections of the solution curves reported in figure 2 due to multiple interpolated initial guesses in Newton sometimes converging to only one of the possible solutions on a branch repeatedly folding back and forth (e.g. both the lower and upper branch pair at a saddle node converging to the lower-branch solution). This could be improved by reducing the arclength step sizes in the original branch continuation. There are two distinct peaks in figure 9 at ${\textit{Re}} = 40$ and ${\textit{Re}} = 100$ , as the arclength continuation was seeded from solutions at these values. The substantial drop in numbers of solutions at the higher ${\textit{Re}}$ regimes commented on in § 3.1 is also apparent – there is a near-exponential drop off in the solutions. Many of the solutions counted in figure 9 are dynamically unimportant and will not contribute to the reconstruction of the turbulent statistics.

Figure 9. Number of unique RPOs converged by slicing the arclength continuation at discrete values of ${\textit{Re}}$ .

We compute a fixed set of weights at a given ${\textit{Re}}$ , $\{ w_{\!j}({\textit{Re}}) \}_{j=1}^{N_{\!p}({\textit{Re}})}$ , where $N_{\!p}({\textit{Re}})$ is the total number of RPOs found at that ${\textit{Re}}$ value. The ansatz, inspired by periodic orbit theory (Artuso et al. Reference Artuso, Aurell and Cvitanovic1990a ; Cvitanović Reference Cvitanović2013), is that the normalised distribution $p_{\boldsymbol{w}}$ of any observable $\gamma : \mathcal{M} \to \mathbb{R}^n$ (here, $\mathcal M$ is the state space) can be constructed via a linear superposition of the RPO statistics for that same observable:

(3.11) \begin{equation} p_{\boldsymbol{w}}(\gamma) = \frac {\sum _{j=1}^{N_{\!p}} w_{\!j} p_{\!j}(\gamma)}{ \sum _{j=1}^{N_{\!p}} w_{\!j} }, \end{equation}

where $p_{\!j}(\gamma)$ is the distribution of the $j\textrm{th}$ RPO and we are suppressing all dependencies on ${\textit{Re}}$ for clarity. We find the weights by minimising the Kullback–Leibler (KL) divergence between the reconstructed observable distribution and the distribution of the target turbulent statistic $q$ ,

(3.12) \begin{equation} D_{\textit{KL}}(p_{\boldsymbol{w}} \| q) = \int p_{\boldsymbol{w}}(\gamma) \log \left ( \frac {p_{\boldsymbol{w}}(\gamma)}{q(\gamma)}\right) \,\textrm{d}\gamma . \end{equation}

Optimisation of the weights proceeds via gradient-based optimisation of $D_{\textit{KL}}$ , using the Adam optimiser (Kingma & Ba 2015) with an initial learning rate of $10^{-2}$ . To ensure that the weights $w_{\!j} \in [0,\infty)$ , we write them in terms of a ‘softplus’ function,

(3.13) \begin{equation} w_{\!j} = \sigma _+(\hat {w}_j) := \log (1 + e^{\hat {w}_j}) \in \mathbb{R}_+, \end{equation}

and optimise for the latent representation $\hat {w}_j \in \mathbb{R}$ . The softplus function is a smooth approximation to the typical rectifier or ReLU function and ensures that the weights are strictly positive. The stopping criterion for convergence is set to

(3.14) \begin{align} \frac {D_{\textit{KL}}^i - D_{\textit{KL}}^{i-1000}}{D_{\textit{KL}}^i} \lt 10^{-8}, \end{align}

where $D_{\textit{KL}}^i$ denotes the KL divergence at the $i\textrm{th}$ optimiser iterate. This optimisation amounts to fitting to a one-dimensional curve and is almost instantaneous on modern hardware – negligible compared with the cost of converging the underlying solutions.

A key feature of the weights obtained via periodic orbit theory is that they are fixed by properties of the underlying exact solutions and do not vary between statistics. Accordingly, robustness of the weights determined in the data-driven approach can be checked by computing the weights using a single ‘training’ statistic, using these fixed weights to reconstruct other ‘test’ statistics at the same ${\textit{Re}}$ , and then comparing KL divergences. Three observables are considered here as training statistics – the normalised dissipation rate $D / D_{lam}$ , normalised production rate $I / D_{lam}$ and the normalised energy $E / E_{lam}$ .

The KL divergence of the three observables are shown as a function of ${\textit{Re}}$ in figure 10, setting the training observable (unfilled markers) to one of the three statistics in each panel. The KL divergences are normalised by the KL divergence of the training observable at ${\textit{Re}} = 100$ . Both $D$ and $E$ (figure 10 a,c) as training statistics lead to relatively robust reconstruction of statistics in the range $40\lesssim {\textit{Re}} \lesssim 150$ , with a noticeable increase in $D_{\textit{KL}}$ for all p.d.f.s beyond this point. This observation is consistent with the decreasing number of dynamically relevant solutions estimated in figure 2. The least robust observable appears to be $I / D_{lam}$ , exhibiting the largest discrepancy between the training and test reconstructions of typically two orders of magnitude across all values of ${\textit{Re}}$ . The reconstruction above ${\textit{Re}} = 300$ is not shown due to the lack of dynamically important solutions.

Figure 10. $D_{\textit{KL}}$ loss for the reconstruction of $D / D_{lam}$ (blue circles), $I / D_{lam}$ (orange squares) and $E / E_{lam}$ (green diamonds) as a function of ${\textit{Re}}$ , each normalised by $D_{\textit{KL}}$ at ${\textit{Re}} = 100$ of the training observable, which was set to $D$ , $I$ and $E$ respectively in panels (a)–(c). Unfilled markers denote this training observable in each panel. The weights were then fixed when reconstructing the distributions of the other two test observables. The normalisation constants are the laminar dissipation rate $D_{lam} = {\textit{Re}} / (2n^2)$ and the laminar energy $E_{lam} = {\textit{Re}}^2 / (4n^4)$ .

Unlike more standard error metrics, the KL divergence (which is not a metric) is not straightforward to interpret, and the performance of the reconstruction is assessed further in figure 11, where the p.d.f.s are compared with direct numerical simulation (DNS) ground truth, with $D$ as the training statistic. As expected, the reconstruction of $D$ itself is fairly well represented by the RPOs – particularly at ${\textit{Re}} \in \{40, 100\}$ where the solutions were originally converged – but with some gaps in the higher dissipation tails. At other ${\textit{Re}}$ values, the reproduction of the high-dissipation tails deteriorates. This can be attributed to the emerging dynamical irrelevance of solution class 1, which climb to increasingly high dissipation rates to connect to inviscid solutions and which may be outside of the turbulent attractor. Unsurprisingly, the reconstruction of the more probable production rates is also fairly robust, though the high- $I$ tail is only really apparent from the RPOs at ${\textit{Re}}=40$ . However, reconstruction of the kinetic energy p.d.f. is noticeably poor compared with $D$ and $I$ , particularly outside of the ${\textit{Re}}$ where solutions were originally converged.

Figure 11. Reconstruction of turbulence statistics for increasing ${\textit{Re}}$ from the top row to the bottom row, with weights obtained using $D / D_{lam}$ as the training observable. (a) From left to right, the reconstructed normalised dissipation rate $D / D_{lam}$ , normalised production rate $I / D_{lam}$ , and normalised energy $E / E_{lam}$ are shown in grey, while the true turbulent distributions are shown in blue. (b) In the left panels, the true turbulent mean velocity profile $\overline {U}(y)$ is shown in blue, while the reconstructed profile is shown in black. In the right panels, the true turbulent r.m.s. velocity fluctuations $u_{rms}$ (blue continuous line), $v_{rms}$ (orange dotted line), averaged over the streamwise direction, discrete symmetries and time, are shown, while the reconstructed profiles are shown in black.

Notably, the p.d.f. reconstruction fails to represent the high- $I$ tails or to reproduce the negative values of production observed at higher ${\textit{Re}}$ values. This is perhaps surprising given the relatively robust reconstruction of the dissipation and the fact that $\overline {D} = \overline {I}$ when averaging is over any relative periodic solution. However, note that the production p.d.f. is significantly wider than that of dissipation. For instance, peak values of $I$ for ${\textit{Re}}\gtrsim 100$ increase to more than twice the peak dissipation values (see p.d.f.s in figure 11). Therefore, a periodic orbit with time-averaged dissipation $\overline {D}$ that samples these extreme $I$ values will need to spend longer in the low- $I$ region. The absence of this behaviour in our solution library is therefore likely due to the relatively short periods of the converged cycles.

Prediction of a full p.d.f. is a particularly challenging task – though one common feature for all reconstructed distributions is that continuation of the RPOs from other ${\textit{Re}}$ values provides coverage for the more probable events. To understand the utility of the solution library at estimating point statistics, we also report in table 1 the relative error in reconstruction of the first two moments of the $D$ , $I$ and $E$ p.d.f.s using the RPOs. These estimates were computed from the distributions reported in figure 11, and show that estimates of the mean of all three quantities is relatively robust at all ${\textit{Re}}$ except ${\textit{Re}}=300$ : generally, errors of less than 5 % are observed. The variance is not so impressive, particularly for the energy, which is reflective of the fact that the solutions representing the rarer excursions to, for example, high-dissipation or high-energy regions of the state space appear to be relevant only in a small- ${\textit{Re}}$ range. This is true even in the ‘asymptotic’ regime beyond ${\textit{Re}} \gtrsim 100$ – for instance, see the relative lack of solutions in the tails at ${\textit{Re}}=160$ in figure 11 compared with those at ${\textit{Re}}=100$ .

Table 1. Relative error in predictions of the first two moments (mean $\hat {\mu }$ and variance $\hat {\sigma }^2$ ) in the distributions reported in figure 11, where the dissipation was used as the ‘training’ observable. Relative error of a statistical estimate is defined as $\varepsilon (\hat {\mu }) := |\hat {\mu } - \mu |/\mu$ , where the ground truth value (here $\mu$ ) is obtained from the ‘true’ turbulent distributions shown in blue in figure 11.

Similarly, predictions for the mean velocity and r.m.s. velocity fluctuations shown in figure 11 (averaged over the streamwise direction, the 16 discrete symmetries of $n=4$ Kolmogorov flow and time, see Chandler & Kerswell Reference Chandler and Kerswell2013; Farazmand Reference Farazmand2016) are also relatively robust over the range of ${\textit{Re}}$ considered. This is perhaps to be expected, as estimating the first moments of the velocities is a more straightforward task than a full p.d.f., as discussed previously. All three predictions at ${\textit{Re}} = 40$ are very close to the true DNS profiles and the mean velocity profile is predicted quite faithfully at all ${\textit{Re}}$ considered. Note the significant weakening of $\overline {U}$ , which was commented on in § 2.2, consistent with the increased decoupling from the background forcing (for reference, the laminar velocity profile is proportional to ${\textit{Re}}$ in our scaling). However, the fluctuating velocities grow with increasing ${\textit{Re}}$ . Predictions for the r.m.s. velocity fluctuations suffer slightly for ${\textit{Re}} \gt 40$ , but errors are consistently below $\mathcal{O}(10\,\%)$ . These velocity profile reconstructions are roughly in line with the observable-agnostic predictions at ${\textit{Re}} \in \{40, 100\}$ by Page et al. (Reference Page, Norgaard, Brenner and Kerswell2024b ).

Tracking how the weights for the RPOs change along their solution curves provides an indication of the changing dynamical importance of that solution as ${\textit{Re}}$ increases. This is shown for a sample of representative solution branches in figure 12(a), where $D$ was again used as the training statistic to fit the weights. The relative weight $w_{\!j} \!/\! \max _j w_{\!j}$ for each RPO is indicated by the size of its marker (circles), while points where the weight $w_{\!j} \!/\! \max _j w_{\!j} \lt 10^{-4}$ are indicated with a cross. The figure indicates the same broad trends identified previously, that weights are larger when the RPO sits closer to the centre of the turbulent distribution (see orange curve in figure 12) and drop as the RPO moves to higher dissipation values to become less dynamically relevant (see green and blue curves – these are class 1 solutions) – though note that the weights tend to vary non-monotonically along the branch, which is to be expected given the emergence/disappearance of solutions in folds and the occurrence of other bifurcations not documented/analysed here.

Figure 12. (a) Visualisation of the changing weights along the solution curves of some example RPOs, with $D$ as the training statistic. The RPOs are visualised via their time-averaged, scaled dissipation rate, and the relative size of each circle along the curves indicates the relative weight $w_{\!j} \!/\! \max _j w_{j}$ of the RPO at that particular ${\textit{Re}}$ . Each colour represents a different solution branch; larger circles indicate larger weights, while crosses denote RPOs with $w_{\!j} \!/\! \max _j w_{j} \lt 10^{-4}$ . (b–d) Dependence of the weights on the real part of the sum of unstable Floquet exponents $\sum _j \sigma _j, \, \sigma _j \gt 0$ , for each RPO at ${\textit{Re}} = 40, 100, 200$ , respectively. The dotted black line denotes the line of best fit, indicating an inverse dependence of $(\sum _i \sigma _i)^{-1.43}$ , $(\sum _i \sigma _i)^{-0.21}$ and $(\sum _i \sigma _i)^{-2.19}$ , at ${\textit{Re}}=40,100,200$ , respectively. Only the results at ${\textit{Re}}=40$ appear to show a clear trend. A total of 163, 22 and 28 RPOs are cut-off, respectively, due to the truncation of the weights at $10^{-4}$ .

The data-driven approach to assigning weights to solutions is driven by an acceptance that the solution library is incomplete. One major flaw in our incomplete library of solutions is that no search was made for pre-periodic solutions with a finite number of shift-reflects in the vertical. We have also seen that many solutions which were believed to be on the attractor (overlap with turbulence in low-dimensional projections) may actually be outside, but still somehow contain ‘turbulent’ dynamical processes. An examination of the dependence of the weights on the Floquet multipliers of the RPOs (figure 12 b–d) indicates the expected inverse correlation $w_{\!j} \sim (\sum _i \sigma _i)^{-1}$ (see also the discussion of Redfern et al. Reference Redfern, Lazer and Lucas2024) for moderate ${\textit{Re}} = 40$ , but shows little correlation at higher ${\textit{Re}}$ . The indication is a decoupling of the fitting from any underlying dynamical properties of the solutions, which is perhaps to be expected given both the fate of classes 1 and 3 as ${\textit{Re}} \to \infty$ and the diminished cover of the tails of the p.d.f.s.

It is pertinent at this point to question the utility of a periodic orbit statistical reconstruction in light of these results. The library of solutions at ${\textit{Re}}=100$ (for example) is, to our knowledge, the largest set of (apparently ‘dynamically relevant’) solutions assembled and is the result of thousands of GPU hours of computation (as documented by Page et al. Reference Page, Holey, Brenner and Kerswell2024a , Reference Page, Norgaard, Brenner and Kerswellb ) augmented by our continuation effort. The success of future efforts to accurately estimate statistics beyond mean values rests both on convergence of solutions with longer periods, bringing additional computational challenges, and an assessment of whether these states are embedded in the turbulent attractor, where we have seen that an inspection of low-dimensional projections is likely very misleading. One caveat here is that this analysis is based on two-dimensional turbulence and may not apply directly to the three-dimensional case. Given the apparent connection of many of the converged RPOs to inviscid solutions, we explore now whether the key dynamical features of the large set of ${\textit{Re}}=100$ solutions can be captured in a much simpler point vortex model.

4.1. Vortex initialisation

The first stage in the labelling procedure is the initialisation of a point vortex system from the reference turbulent RPO at ${\textit{Re}} = 100$ . First, $N_{S} = 50$ snapshots are sampled from the reference Kolmogorov RPO, equally spaced in time. The vortex cores in each snapshot are extracted using the vortex identification method outlined in § 2.2 and the modal number of vortex cores at each snapshot throughout the period of the RPO, $\bar {N}_v$ , is computed.

This procedure seeks to initialise a point vortex system that matches the strongest vortex cores over the period of the turbulent RPO. Trajectories of vortex cores in the reference RPO are tracked by connecting vortex cores $i$ and $j$ in sequential snapshots $k$ and $k+1$ which satisfy the similarity condition

(4.1) \begin{equation} \frac {\left\| \boldsymbol{X}_i^k - \boldsymbol{X}_j^{k+1} \right\|_2}{\| \boldsymbol{L} \|_2} \lt 0.1, \end{equation}

where $\boldsymbol{X}_i = (\tilde {x}_i, \tilde {y}_i, \varGamma _i, A_i)$ is a state vector collating the centre of vorticity, circulation and area of each vortex core. The normalisation $\boldsymbol{L}$ is set to $(L_x, L_y, \varGamma _i, A_i)$ , rather than $\boldsymbol{X}_i$ . This choice of normalisation is made so that the similarity condition is independent of the location of the vortex cores on the domain, and depends instead on the inter-core distance – note that periodic boundary conditions are accounted for when computing the difference between $\tilde {\boldsymbol{x}}_i^k$ and $\tilde {\boldsymbol{x}}_j^{k+1}$ . As the turbulent solutions are typically RPOs, there is a translational discontinuity between the first and final snapshot of the period of the solution. This discontinuity is taken into account by translating the vortex cores in the first snapshot by the reverse translational shift when attempting to connect a trajectory between the first and final snapshots.

Some vortex cores will not satisfy the vorticity thresholding condition in (2.6) for all snapshots along their trajectories. In general, very few vortex-core trajectories exist over the full period of the reference RPO due to the fluctuating strength of the vortices. The trajectories of the vortex cores are sorted by their average circulation over their existence,

(4.2) \begin{equation} \langle \varGamma _i \rangle = \frac {\int _{0}^T \varGamma _i(t) \mathbb{I}_{i}(t) \,\textrm{d}t}{\int _{0}^T \mathbb{I}_{i}(t)\, \textrm{d}t}, \end{equation}

where the indicator function $\mathbb{I}_{i}(t)$ equals 1 if vortex core trajectory $i$ exists at time $t$ , and 0 otherwise. The average strength $\langle \varGamma _i \rangle$ is used as a measure of the dynamical importance of the $i\textrm{th}$ vortex core trajectory. The $\bar {N}_v + 1$ most dynamically important trajectories are extracted and $\langle A_i \rangle$ is also computed for each of these trajectories. The areas and circulations of the point vortex state are initialised with $\langle A_i \rangle$ and $\langle \varGamma _i \rangle$ , respectively (the areas are required to compare to the reference RPO evolution – see later). If all of the extracted trajectories co-exist in some snapshot of the reference RPO, the centres of vorticity of the vortex cores in this snapshot are used to initialise the positions of the point vortices. Otherwise, the snapshot with the smallest number of absent vortices is used and the positions of the absent vortices are initialised with their average centres of vorticity $\langle \tilde {\boldsymbol{x}}_i \rangle$ .

4.2. Convergence

The dynamical initialisation algorithm outlined in the previous subsection outputs the position, circulation and area of $\bar {N}_v + 1$ point vortices. These point vortices correspond to the $\bar {N}_v + 1$ strongest vortex core trajectories extracted from the reference RPO. However, there is no guarantee that this initialisation will be near to an exact point vortex solution. The initialisation is also very sensitive to the vorticity extraction condition in (2.6). This is because vortex cores will not satisfy (2.6) for sections of their trajectories, directly impacting the modal number of vortex cores extracted $\bar {N}_v$ . For these reasons, we employ gradient-based optimisation to improve the point vortex fit and reduce the sensitivity to the extraction method, before converging onto an exact point vortex RPO with a Newton-GMRES-hookstep solver.

The gradient-based optimisation seeks to enforce two specific behaviours in the point vortex evolution: (i) that the point vortices must shadow the vortex cores from the reference RPO closely and (ii) that the point vortex evolution itself must also be (relative) time periodic. We construct a scalar loss function $\mathscr{L}$ to target these two effects. Unlike in the Newton-GMRES-hookstep solver subsequently used to converge exact solutions, the circulations of the point vortices are allowed to change during optimisation.

Matching to the reference vortex cores is done by a contribution $\mathscr{L}_{\textit{match}}$ to the total scalar loss function $\mathscr{L}$ . The contribution $\mathscr{L}_{\textit{match}}$ is included to attempt to minimise the difference between the reference RPO and a grid-based representation of the point vortex solution. The singular point vortex circulations of an $N_v$ -vortex configuration at position $\textbf{x}^* = [\boldsymbol{x}^*_1, \ldots, \boldsymbol{x}^*_{N_v}] \in \mathbb{R}^{2 N_v}$ (note $\boldsymbol x_{\alpha }^* := [x_{\alpha }^*, y_{\alpha }^*]$ ), with circulations $\boldsymbol{\varGamma } \in \mathbb R^{N_v}$ and areas $\boldsymbol{ A} \in \mathbb R^{N_v}$ , are used to define a doubly periodic observable by placing a Gaussian at each vortex location:

(4.3) \begin{equation} G\left(\boldsymbol{x}; \textbf{x}^*, \boldsymbol \varGamma\right) := \sum _{\alpha =1}^{N_v} g_\alpha \left(\boldsymbol{x}; \boldsymbol x^*_{\alpha }, \varGamma _{\alpha }\right)\!, \end{equation}

where

(4.4) \begin{align} g_\alpha \left(\boldsymbol{x}; \boldsymbol x^*_{\alpha }, \varGamma _{\alpha }\right) := \sum _{n,m = -\infty }^{\infty } \tilde {g}_{\alpha }\left(\boldsymbol{x}; \boldsymbol x^*_{\alpha } + (nL, mL), \varGamma _{\alpha }\right) \end{align}

is a doubly periodic function with Gaussian contributions from all image vortices, i.e.

(4.5) \begin{align} \tilde {g}_{\alpha }\left(\boldsymbol{x}; \boldsymbol x^*_{\alpha }, \varGamma _{\alpha }\right) = \frac {\varGamma _{\alpha }}{2\pi \sigma _{\alpha }^2} \exp \left (\frac {-\left(\boldsymbol x - \boldsymbol x^*_{\alpha }\right)^2 }{2\sigma _{\alpha }^2}\right)\!. \end{align}

In this expression, the properties of the vortex cores extracted from the reference turbulent RPO are used to define the standard deviation of the Gaussian, $\sigma _{\alpha }^2 = 0.1 A_{\alpha }$ . This smoothing of the vorticity delta functions centred on each point vortex allows for comparison between snapshots of the point vortex system and the vorticity field of the reference Kolmogorov RPO. It also allows for the possibility of vortex permutations if their strengths and areas are equal. To compute the doubly periodic Gaussian distribution in (4.3), we approximate its Fourier series with the Fourier transform of the Gaussian distribution on the infinite domain

(4.6) \begin{equation} \frac {1}{4\pi ^2}\int _{-\pi }^{\pi } \int _{-\pi }^{\pi } g_{\alpha }(\boldsymbol{x}) e^{-i\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{x}}\,\textrm{d}^2\boldsymbol{x} \approx \frac {1}{4\pi ^2}\int _{-\infty }^{\infty } \int _{-\infty }^{\infty } \tilde {g}_{\alpha }(\boldsymbol{x}) e^{-i\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{x}}\,\textrm{d}^2\boldsymbol{x}, \end{equation}

which we find to be numerically accurate to machine precision when $\sigma _{\alpha }^2 \lt 2$ . The Fourier transform approximation $\hat {g}_{\alpha }(\boldsymbol{k})$ is derived by applying the shift and stretching Fourier theorems to the standard Gaussian Fourier transform

(4.7) \begin{equation} \hat {g}_{\alpha }(\boldsymbol{k}) = \frac {1}{4\pi ^2}e^{-i\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{x}_{\alpha }} e^{-\boldsymbol{k}^2\sigma _{\alpha }^2 / 2} \, . \end{equation}

Taking the inverse Fourier transform of (4.7) yields the approximation for the doubly periodic Gaussian distribution in (4.3).

The function $\mathscr{L}_{\textit{match}}$ then minimises the average $L_2$ distance between snapshots extracted from the turbulent RPO $\{\omega (\boldsymbol{x}, t) : 0 \leqslant t \lt T_{\textit{RPO}}\}$ and the Gaussian observable of the point vortex solution over the period $T_{\textit{RPO}}$ of the turbulent solution,

(4.8) \begin{equation} \mathcal{L}_{\textit{match}}(\textbf{x}^*, \boldsymbol{\varGamma }) := \frac {1}{N_{S}} \sum _{i=1}^{N_{S}} \frac {\left \| G(\boldsymbol{x}; \boldsymbol{f}^{t_i}(\textbf{x}^*), \boldsymbol{\varGamma }) - \omega (\boldsymbol x, t_i)\right \|^2}{\left \| \omega (\boldsymbol x, t_i)\right \|^2}, \end{equation}

where $\boldsymbol{f}^t(\textbf{x}^*)$ is the time forward map of (2.15), $t_n = n T_{\textit{RPO}} / N_{S}$ are equally spaced times over the period and the norm $\| v \|^2 := (1/V)\int _V v^2 \, \textrm{d}^2 \boldsymbol x$ . The period and vortex areas $\boldsymbol{A}$ are kept fixed throughout the optimisation, so the only adjustable parameters for $\mathscr{L}_{\textit{match}}$ are the initial vortex positions $\textbf{x}^*$ and their circulations $\boldsymbol \varGamma$ .

The second contribution to the overall loss function $\mathscr{L}$ is designed to target point vortex RPOs, and is termed $\mathscr{L}_{\textit{RPO}}$ . To allow for permutations of the point vortices between the states $\textbf{x}^*$ and $\boldsymbol{f}^{T_{\textit{RPO}}}(\textbf{x}^*)$ , this loss function is also defined in terms of the periodic Gaussian observable in (4.3),

(4.9) \begin{equation} \mathcal{L}_{\textit{RPO}}(\textbf{x}^*, \boldsymbol{\varGamma }, \boldsymbol s) = \frac {\left \| G(\boldsymbol{x}; \mathcal{T}^{\boldsymbol s} \boldsymbol{f}^{T_{\textit{RPO}}}(\textbf{x}^*), \boldsymbol{\varGamma }) - G(\boldsymbol{x}; \textbf{x}^*, \boldsymbol{\varGamma })\right \|^2}{\left \| G(\boldsymbol{x}; \textbf{x}^*, \boldsymbol{\varGamma })\right \|^2}, \end{equation}

where $\mathcal{T}^{\boldsymbol s} = \mathcal{T}^{(s_x, s_y)}$ generates translational shifts in the $x$ - and $y$ -directions. While the turbulent RPOs cannot shift in $y$ (we do not have solutions in our library which shift-reflect), the point vortex states are allowed to drift in both directions. The adjustable parameters are $\textbf{x}^*$ and $\boldsymbol \varGamma$ (as in $\mathscr{L}_{\textit{match}}$ ), as well as the shifts $s_x$ and $s_y$ .

The full loss function is a combination of the two terms discussed previously:

(4.10) \begin{equation} \mathscr{L}(\textbf{x}^*, \boldsymbol{\varGamma }, \boldsymbol s) = \kappa \beta \mathscr{L}_{\textit{match}}(\textbf{x}^*, \boldsymbol{\varGamma }) + (1 - \kappa) \mathscr{L}_{\textit{RPO}}(\textbf{x}^*, \boldsymbol{\varGamma }, \boldsymbol s), \end{equation}

where the hyperparameter $\beta = 50$ is set to a large number to increase the relative importance of the matching loss term. This large weighting of $\mathscr{L}_{\textit{match}}$ places a heavy emphasis on finding point vortex solutions which closely mimic the dominant vortex dynamics, which are themselves time periodic. The second term $\mathscr{L}_{\textit{RPO}}$ should be thought of as a small correction to aid in this convergence towards exactly recurring solutions. The hyperparameter $\kappa$ is defined to conveniently set either of the loss terms to zero.

The gradients $\boldsymbol {\nabla} _{\textbf{x}^*, \boldsymbol{\varGamma }, \boldsymbol s}\mathscr{L}$ are passed to an Adam optimiser (Kingma & Ba 2015) with an initial learning rate of $\eta = 10^{-2}$ . The optimiser is allowed to run for a maximum $N_{\textit{opt}} = 1000$ steps, or until $\mathscr{L}_{\textit{RPO}} \leqslant 10^{-5}$ . The value of the hyperparameter $\kappa$ has a strong effect on the likelihood of satisfying this tolerance on $\mathscr{L}_{\textit{RPO}}$ . For example, for $\kappa = 0.5$ , just 4 % of initialisations reach this tolerance (the $\mathscr{L}_{\textit{match}}$ term is still heavily weighted due to our choice of $\beta$ ), compared with the 74 % of initialisations for $\kappa = 0$ . All the optimiser outputs are passed to a Newton solver to attempt convergence regardless of whether they reach this tolerance level.

Two non-differentiable operations are then performed on the output of the gradient-based optimisation. The first is to remove any vortices which satisfy $|\varGamma _{\alpha }| \lt \varepsilon$ , where $\varepsilon$ is a tolerance chosen as one of $\varepsilon \in \{0, 0.01, 0.05\}$ . For reference, vortex circulations were typically in the range $| \varGamma _{\alpha }| \in [1,5]$ after the gradient-based optimisation. As $\mathcal{L}_{\textit{RPO}}$ is computed in terms of the Gaussian observable, the optimiser tends to reduce the circulation of any point vortices with a large contribution to $\mathcal{L}_{\textit{RPO}}$ . Removing these weak vortices improves the final convergence rate of the Newton-GMRES-hookstep solver.

The second non-differentiable operation is the identification of any vortices which permute their positions over the period of the point vortex solution. As $\mathcal{L}_{\textit{RPO}}$ is computed in terms of the grid-based Gaussian observable, vortex cores with similar areas and circulations are allowed to permute over the period $T_{\textit{RPO}}$ while maintaining small $\mathscr{L}_{\textit{RPO}}$ loss. As such, we have not yet needed to be explicit about the permutation symmetry of point vortices with equal circulation. However, the Newton-GMRES-hookstep solver used in the final convergence step finds roots of the residual defined in terms of the vortex positions

(4.11) \begin{equation} \boldsymbol{F}(\textbf{x}^*, \boldsymbol{s}, T) = \mathcal{T}^{\boldsymbol{s}}\mathbb{P}\,\boldsymbol f^T(\textbf{x}^*) - \textbf{x}^*, \end{equation}

where the permutation $\mathbb{P}$ must be computed a priori. A modified Jonker–Volgenant variant (Crouse Reference Crouse2016) of the Hungarian algorithm (Kuhn Reference Kuhn1955) is used to compute the vortex permutation $\mathbb{P}$ between $\textbf{x}^*$ and $\boldsymbol f^{T_{\textit{RPO}}}(\textbf{x}^*)$ after the optimisation and prior to the initialisation of the Newton solver (Cleary & Page Reference Cleary and Page2023). The simplest permutation is a pairwise swap of two vortices $\boldsymbol{x}^*_{\alpha }$ and $\boldsymbol{x}^*_{\beta }$ , which can be written as a 2-cycle, $\mathbb{P} = (\alpha \, \beta)$ . However, higher-order $k$ -cycles which permute $k$ vortices are also possible in point vortex RPOs. For the state $\mathcal{T}^{\boldsymbol{s}}\mathbb{P}\,\boldsymbol f^T(\textbf{x}^*)$ to be a symmetric copy of $\textbf{x}^*$ , all vortices belonging to the $k$ -cycle must share the same circulation. Any permutation can be decomposed into disjoint cycles $\xi _i$ , such that $\mathbb{P} = \prod _i \xi _i$ and each vortex belongs to at most one $\xi _i$ . To ensure that the permutation symmetry is respected, the circulation of all vortices belonging to each $\xi _k$ is then manually set to their mean circulation. All vortex circulations are then fixed and not further updated by the Newton solver.

When converging solutions with the Newton solver, we apply a threshold on the Newton residual and a condition that the period is not significantly adjusted after optimisation:

(4.12) \begin{align} \frac { \left \| \mathcal{T}^{\boldsymbol{s}}\mathbb{P}\,\boldsymbol f^T(\textbf{x}^*) - \textbf{x}^* \right \|}{\| \textbf{x}^* \|} &\lt 10^{-8}, \end{align}

(4.13) \begin{align} \frac { \left | T - T_{\textit{RPO}} \right |}{ T_{\textit{RPO}} } &\lt 0.25. \end{align}

The optimal initial size of the trust region for the Hookstep solver was found to be $0.01 \times \| (\textbf{x}^*, \boldsymbol s, T) \|_2$ . This small trust region resulted in a greater rate of dynamically representative convergences – larger values resulted in point vortex RPOs that bore little resemblance to the reference solution; smaller initial trust regions required too many Newton steps for convergence.

4.3. Optimisation results

A number of different experiments were run, varying the circulation threshold $\varepsilon$ , controlling which vortices are retained, and the annealing parameter $\kappa$ in (4.10). These experiments are summarised in table 2, while examples of some converged solutions, overlaid over the reference Kolmogorov RPO, are shown in figure 14. The labelling procedure was initialised with the 252 RPOs converged at ${\textit{Re}} = 100$ (151 solutions from the original library, 101 new solutions from the continuation). Three different initialisations were run for each of these reference turbulent RPOs, by extracting the $N_v-1, \, N_v$ and $N_v + 1$ most dynamically important vortex cores over the period of the reference RPO (the overbar indicating the modal number of point vortices is dropped from now on). Vortex numbers are typically $2 \leqslant N_v \leqslant 8$ , with most periodic orbits being well described with $N_v = 4$ or 5. The ‘Total convergences’ column in table 2 reports the total number of convergences (larger than 252 for each experiment due to the three initialisations run for each RPO), while the ‘Total matches’ column reports how many of the reference RPOs yielded at least one converged point vortex RPO.

Figure 14. Out-of-plane vorticity (contours) of three RPOs from Kolmogorov flow at ${\textit{Re}} = 100$ . Snapshots in each row are sampled uniformly in time along the period of that solution. Overlaid on each snapshot are the point vortices at the corresponding (proportional) time along the period of the matched point vortex RPO. The size of the markers is proportional to the circulation of the vortex, and the lines denote the path swept out by each point vortex over its full evolution. Red (white) markers denote positive (negative) circulations. The reference RPO period $T_{\textit{RPO}}$ , reference RPO shift $s_{\textit{RPO}}$ , point vortex RPO period $T$ and point vortex shifts $\boldsymbol s$ for each solution are: (top row) $T_{\textit{RPO}} = 1.34$ , $s_{\textit{RPO}} = 0.0075$ , $T = 1.2$ , $\boldsymbol s = (-0.032, 0.012)$ ; (second row) $T_{\textit{RPO}} = 1.79$ , $s_{\textit{RPO}} = 0.085$ , $T = 1.755$ , $\boldsymbol s = (0.05, 0.006)$ ; (third row) $T_{\textit{RPO}} = 1.27$ , $s_{\textit{RPO}} = 0.69$ , $T = 1.27$ , $\boldsymbol s = (0.48, 0)$ . This point vortex solution is a vortex crystal (relative equilibrium); (bottom row) $T_{\textit{RPO}} = 4.16$ , $s_{\textit{RPO}} = -0.756$ , $T = 4.34$ , $\boldsymbol s = (0.04, -0.12)$ .

Table 2. Summary of point vortex RPO-fit experiments run using 252 turbulent RPOs at ${\textit{Re}} = 100$ to initialise the optimisation. Each reference RPO yields three different initialisations, by extracting the modal number of vortex cores, $N_v$ , as well as $N_v \pm 1$ , over the period of the reference RPO.

Overall, table 2 indicates that all parameter settings achieve very large numbers of converged solutions, though these convergences may in some cases differ substantially from the reference simulation – point vortex states which are ‘similar’ to the underlying turbulent periodic orbit are quantified by a dynamical relevance metric $\mathscr{L}_{\textit{rel}} \leqslant 1$ , which is introduced and discussed later. In some cases in table 2, the parameter $\kappa$ is labelled ‘Anneal’: we use this terminology to describe an optimisation that is initialised with $\kappa = 1$ , and in which $\kappa$ is then decreased by $1/N_{\textit{opt}}$ at each optimisation step. This annealing aims to initially closely match the reference turbulent RPO, while slowly increasing the relative importance of the $\mathscr{L}_{\textit{RPO}}$ throughout the optimisation procedure. This scheduling results in a slight increase in the total number of point vortex RPO convergences compared with a constant equal weighting ( $\kappa = 0.5$ ) of $\mathscr{L}_{\textit{RPO}}$ and $\mathscr{L}_{\textit{match}}$ . Ignoring the $\mathscr{L}_{\textit{match}}$ term (setting $\kappa = 1$ ) results in a significant increase in the total number of point vortex RPO convergences.

The solutions reported in the top three rows of figure 14 all represent cases where the converged point vortex dynamics closely track the dominant vortex cores in the original solution. They include cases featuring a multipolar structure (top row) and a bound state or triangular vortex (middle row) – in both of these examples, the point vortex RPO involves a permutation of vortices at each period. The final example is a crystal-like structure in which two rows of same-signed vortices translate to the right. In all three cases, the period of the converged point vortex solution matches the reference turbulent RPO to within a few percent or better.

However, by visual inspection, the dynamics of some converged point vortex RPOs differs qualitatively from the reference turbulent RPO. An example of which is shown in the bottom row of figure 14, in which a six vortex system has been converged which only slightly resembles the reference RPO – there are two trapped systems of $N_v=4$ and $N_v=2$ vortices which roughly align with large-scale vortices in the turbulent flow.

To quantify the dynamical relevance of the point vortex state to the original turbulent dynamics, we compute the translation-and-time-reduced metric:

(4.14) \begin{equation} \mathscr{L}_{\textit{rel}} := \min _{\boldsymbol{s}, \Delta t} \frac {1}{N_{S}} \sum _{j = 1}^{N_{S}} \frac {\left\| G(\boldsymbol{x}; \mathcal{T}^{\boldsymbol s} \boldsymbol f^{t_j^{v} + \Delta t} (\textbf{x}^*), \boldsymbol{\varGamma }) - \sum _{i} \mathbb{I}_i (\boldsymbol{x},t_j) \omega (\boldsymbol{x}, t_j) \right\|^2}{\left\| \sum _{i} \mathbb{I}_i (\boldsymbol{x},t_j) \omega (\boldsymbol{x},t_j) \right\|^2}, \end{equation}

where the indicator function $\mathbb{I}_i (\boldsymbol{x},t) = 1$ if $\boldsymbol{x} \in \mathcal{V}_{i}$ at time $t$ – i.e. we are comparing only the vortex cores in the original RPO to the point vortex reconstruction. The solutions are compared at equally spaced times $t_j^{v} := j T / N_{S}$ and $t_j = j T_{\textit{RPO}} / N_{S}$ along the periods of point vortex RPO $T$ and the reference turbulent RPO $T_{\textit{RPO}}$ , respectively. Due to potential small mismatches between the alignment of the vortices in the two solutions in space and time, we search in (4.14) over a restricted range of (constant) temporal and spatial shifts: $\Delta t \in \{0, \pm 10 \delta t, \pm 25\delta t, \pm 50\delta t, \pm 100\delta t\}$ , and $s_x, s_y \in \{0, \pm L/100, \pm L/40, \pm L/20, \pm L/10\}$ .

It is worth considering what a ‘good’ value of $\mathscr{L}_{\textit{rel}}$ would be. For unrelated point vortex observables and turbulent RPOs, we find $\mathscr{L}_{\textit{rel}} \approx 2$ (the expected result given we are comparing two near-orthogonal vectors in a high-dimensional space). For the dynamically similar states obtained in figure 14, we obtain values of $\mathscr{L}_{\textit{rel}} = 0.47, 0.49, 0.46$ for the first, second and third rows, respectively. The state in the bottom row of figure 14 only achieves a value of $\mathscr{L}_{\textit{rel}} = 1.58$ , indicating a poor representation of the turbulent dynamics. The relatively ‘large’ values of $\mathscr{L}_{\textit{rel}}$ in the well-matched configurations can be rationalised by the fact that the smoothed point vortex observable is constructed as a superposition of Gaussians, with the shape of the original vortex core modelled by a variance which is proportional to the vortex area, $\sigma _{\alpha }^2 = 0.1 A_{\alpha }$ : we are comparing fields which can differ quite significantly even in the best scenarios. As such, we find that $\mathscr{L}_{\textit{rel}}\leqslant 1$ is a sensible threshold to discern whether the converged point vortex RPO captures critical features of the original turbulent solution. We also checked all point vortex states with $\mathscr{L}_{\textit{rel}} \leqslant 1$ for uniqueness by rescaling variables such that $\max |\varGamma _j| = 1$ , before comparing their rescaled periods and values of the rescaled Hamiltonian. All solutions are unique.

Figure 15. Histograms of the dynamical relevance metric (4.14) for experiments 1c, 2b and 3b detailed in table 2, with a bin size of 0.05. The mean and mode of each experiment are indicated by the correspondingly coloured bold dashed line and filled bars, respectively.

A comparison of the distributions of $\mathscr{L}_{\textit{rel}}$ across three representative experiments from the list in table 2 (1c, 2b and 3b) is reported in figure 15. These point vortex RPO searches all enforced removal of vortices whose circulation shrunk in the optimisation below a threshold of $\varepsilon =0.05$ , and are distinguished from one another by the relative weight placed on the two contributions to the loss (4.10): In run 1c, the ‘RPO’ contribution (4.9) was included alongside the matching constraint (4.8), while it represented the entire loss in run 2b (no matching term). In run 3b, annealing was applied to incrementally increase the contribution of the RPO term over the course of the optimisation. All of the search configurations summarised in figure 15 show a wide spread in $\mathscr{L}_{\textit{rel}}$ . Notably the RPO-only results (2b) show much higher $\mathscr{L}_{\textit{rel}}$ when no requirement is made that the point vortices attempt to shadow the vortex cores. We focus for the remainder of the section on the results of experiment 3b, which has the lowest mean (and mode) in $\mathscr{L}_{\textit{rel}}$ , though the results are representative of the suite of experiments in general. Note that the value of the circulation cutoff, $\varepsilon$ , selected for analysis here (the largest, $\varepsilon = 0.05$ ) leads to an increased number of convergences relative to calculations which retain larger numbers of weak vortices, though the overall distribution of results are not substantially changed from those shown in figure 15.

4.4. Discussion

We now explore the properties of the library of converged point vortex states from experiment 3b and their relationship to the three classes of Kolmogorov RPO discussed in § 3.1. The number of convergences from this experiment satisfying $\mathscr{L}_{\textit{rel}} \lt 1$ is visualised as a function of the number of point vortices $N_v$ in figure 16(a). There is a preference for small numbers of vortices, $N_v \leqslant 5$ , with the most solutions found for $N_v=4$ . Intriguingly, there is no clear dependence in the reconstruction error, $\mathscr{L}_{\textit{rel}}$ , on the number of vortices, $N_v$ , other than an increase in the minimum observed value of $\mathscr{L}_{\textit{rel}}$ with increasing $N_v$ . This is perhaps unexpected as intuitively, one might think that placing more point vortices in the system would result in a more faithful representation of the turbulent dynamics. However, the solutions in the turbulent RPO library at ${\textit{Re}} = 100$ tend to be dominated by a small number of vortex cores (see for example figures 5 and 6), and additional point vortices are attempting to represent weak, small-scale filamentary structures which are poorly represented with Gaussians (discussed further later).

Figure 16. (a) Total number of point vortex RPO convergences at each value of $N_v$ for experiment 3b in table 2 which satisfy $\mathscr{L}_{\textit{rel}} \lt 1$ . (b) Histogram showing the distribution of $\mathscr{L}_{\textit{rel}}$ for experiment 3b as a function of $N_v$ . The bin size for $\mathscr{L}_{\textit{rel}}$ is 0.1.

Some examples of solutions in the most common $N_v=4$ class of point vortex RPOs are reported in figure 17. Most $N_v=4$ solutions appear as an approximately equispaced lattice, with two dominant vortices of opposite circulation. An example of one of these ‘quasi-crystal’ configurations is shown in figure 17(a). Other representative solutions for the $N_v = 4$ system are shown in the remaining two panels, and show a ‘tripole’ structure and an oscillatory lattice, featuring rows of equal-signed vortices, respectively.

Figure 17. Evolution of three samples point vortex RPOs with $N_v = 4$ . Trajectories are plotted as red (blue) points denoting positive (negative) circulation, such that increasing opacity denoting increasing time along the period and the separation of points along a trajectory gives an indication of the speed of that vortex. Streamlines of the final (most opaque) vortex configuration are also shown, such that darker streamlines denote greater local induced speed. Each system is simulated for one complete period. (a) Vortex crystal configuration with $T = 2.16$ , $\boldsymbol s = (-0.008, 0.002)$ . (b) Tripolar configuration with $T = 1.75$ , $\boldsymbol s = (-0.056, 0.073)$ . (c) Uniform crystal lattice configuration with $T = 1.27$ , $\boldsymbol s = (0.48, 0)$ .

A key question is whether the successful fitting of a point vortex solution to a Kolmogorov RPO says anything about the possible connection of the original RPO to a solution of the (unforced) Euler equation. To explore this point, we first examine how many of the subset of RPOs which were successfully labelled with a ‘representative’ point vortex solution at ${\textit{Re}}=100$ (quantified by $\mathscr{L}_{\textit{rel}} \lt 1$ ) can be continued to higher ${\textit{Re}}$ . This is done via the histogram in figure 18(a), which compares the number of labelled branches to the total number of solutions as ${\textit{Re}}$ varies; both scale similarly with ${\textit{Re}}$ .

Figure 18. (a) Histogram of all the distinct RPO branches which cross ${\textit{Re}} = 100$ as a function of ${\textit{Re}}$ in grey. The histogram of the branches which were successfully matched at ${\textit{Re}} = 100$ is shown in red. (b) Energy dissipation rate $D$ against energy production rate $I$ at ${\textit{Re}} = 100$ , both normalised by the laminar dissipation $D_{lam} = {\textit{Re}} / (2n^2)$ , for the turbulent RPOs which were successfully matched with a point vortex RPO. The RPOs are coloured according to their branch class: 1, red; 2, blue and 3, green. The grey background is the p.d.f. computed from a trajectory of $10^5$ samples, separated by 1 advective time unit. The contour levels of the p.d.f. are spaced logarithmically.

More revealing is figure 18(b), which identifies point-vortex-labelled RPOs on a production/dissipation plot at ${\textit{Re}}=100$ , where the solutions are coloured according to the class assigned in § 3.1. The figure appears to be dominated by solutions from classes 1 (red, dissipation stronger than turbulent p.d.f.) and 2 (blue, dissipation within turbulent p.d.f.). However, this figure is distorted by the large number of solutions in class 2, while the labelling success rate as a percentage of the cardinality of the class was significantly larger for class 1 and class 3. These statistics are summarised in table 3, which reports the number and percentage of successfully labelled solution branches (one branch may contain multiple solutions at a given ${\textit{Re}}$ ) for each of the three distinct branch classes. The branches which could be continued the highest in ${\textit{Re}}$ by the arclength continuation (class 3) were the most likely to be successfully labelled with a point vortex RPO (albeit with only a small sample size). The RPOs in class 1 were the second most likely to be successfully labelled.

Table 3. Number and percentage of successfully labelled ( $\mathscr{L}_{\textit{rel}} \lt 1$ ) branches for each of the three distinct classes of solution branches identified in § 3.1. The ${\textit{Re}} = 100$ columns report these figures for experiment 3b, run on the library of 252 RPOs at ${\textit{Re}} = 100$ . The ${\textit{Re}} \gt 200$ columns report these figures for the same labelling procedure, but the reference solutions considered are instead all the terminal branch solutions which were continued to above ${\textit{Re}} = 200$ .

Even if a solution is successfully labelled, another important consideration is whether the extracted ‘vortices’ which are used to define the point vortex initialisation (and subsequently measure the success of the ‘reconstruction’) remain as ${\textit{Re}}$ increases, or whether they thin and shrink. The latter behaviour is expected to occur for filamentary structures where point vortices are a poor modelling choice (see figure 1).

To probe this, we report in figure 19 p.d.f.s of the extracted vortex areas, $\{A\}$ , for the Kolmogorov RPOs for which a point vortex solution was also found. P.d.f.s are shown at both ${\textit{Re}}=100$ where the extraction + fit was performed (colour) and at the highest ${\textit{Re}}$ values on the solution curve (grey). The results show that the area of the vortices for class 1 solutions is roughly maintained (and on average increases slightly) on continuation. However, the class 2/3 solutions exhibit on average decreasing vortex area – though the continuation does not generally proceed particularly high in ${\textit{Re}}$ for most class 2 states. Class 1 also indicates a bimodal distribution with a peak at a large vortex area.

Figure 19. P.d.f.s of vortex areas (thresholding used to define this variable is described in § 2.2) in the converged RPOs for which a point vortex RPO was found. (a) Class 1 solutions: red is the p.d.f. at the starting value ${\textit{Re}}=100$ , grey is the p.d.f. obtained using solutions at the highest ${\textit{Re}}$ on the branch. (b) Classes 2 and 3: cyan is the p.d.f. at the starting value of ${\textit{Re}}=100$ , grey the p.d.f. at the highest ${\textit{Re}}$ on the branch (note only states which were continued successfully beyond ${\textit{Re}} \gt 150$ were included in the class 2/3 results here).

Figure 20. (a–d) Contours of vorticity for four RPOs in Kolmogorov flow. Top row shows a snapshot at ${\textit{Re}}=100$ along with the point vortex RPO overlaid as in previous figures. The bottom row is the final state obtained after arclength continuation up in ${\textit{Re}}$ . The two leftmost RPOs belong to ‘class 1’, while the two rightmost RPOs belong to ‘class 3’ as described in § 3.1. The reference RPO period and shift, $T_{\textit{RPO}}$ , shift $s_{\textit{RPO}}$ , point vortex RPO period $T$ and shifts $\boldsymbol{s}$ for each labelled solution on the top row are: (a) $T_{\textit{RPO}} = 2.99$ , $s_{\textit{RPO}} = -0.44$ , $T = 2.97$ , $\boldsymbol s = (-0.041, -0.023)$ ; (b) $T_{\textit{RPO}} = 1.97$ , $s_{\textit{RPO}} = -0.43$ , $T = 1.89$ , $\boldsymbol s = (0.01, -0.01)$ ; (c) $T_{\textit{RPO}} = 4.61$ , $s_{\textit{RPO}} = 0.071$ , $T = 4.695$ , $\boldsymbol s = (-0.008, -0.001)$ ; (d) $T_{\textit{RPO}} = 4.67$ , $s_{\textit{RPO}} = -0.17$ , $T = 4.497$ , $\boldsymbol s = (0.016, -0.018)$ . The point vortex solution in panel (d) is a travelling wave. The final solutions on the bottom row were converged at ${\textit{Re}} = 185.69, 175.003, 345.297$ and $277.23$ , from panel (a–d), respectively.

The results in figure 19 should be interpreted alongside the labelled vortex areas shown earlier in figure 1: filamentary structures at ${\textit{Re}}=100$ tend to have $A \lesssim 0.2$ with some exceptions, while these structures are not strong enough to be flagged at the higher value of ${\textit{Re}}=400$ . The behaviour of the class 1 area p.d.f.s under increasing ${\textit{Re}}$ is consistent with an absence of filaments and the dominance of larger-area vortex cores, a result supported by two examples from both classes 1 and 3 shown in figure 20, where the point vortex solution is overlaid. Also included is a higher- ${\textit{Re}}$ state from the same solution branch. In the cases from class 1, the point vortex periodic orbits correspond to features of the turbulent solution which remain as ${\textit{Re}}$ grows while the thin background filaments weaken in comparison. For instance, note the multipole structure in the second column approximating ring of vorticity, which amplifies under increasing ${\textit{Re}}$ , while the weaker triangular vortex in the same snapshot is not represented in the point vortex solution, and weakens as ${\textit{Re}}$ increases. Similar effects are also apparent in the earlier examples in figure 14 from class 1 (top 3 rows of that figure).

In contrast, the solutions in class 3 have resulted in the initialisation of weaker point vortices on small filamentary regions (e.g. see the weaker pair of vortices in the second class 3 example in figure 20(d)), a feature hinted at by the relative importance of low area vortices in figure 19. In these cases, it appears we have found a point vortex RPO which can mimic the weaker viscous dynamics, without implying anything stronger about the behaviour as ${\textit{Re}}\to \infty$ .

The identification by the point vortex fit of structures in class 1 which persist as ${\textit{Re}}$ increases is consistent with the connection of these RPOs to exact solutions of the Euler equation as ${\textit{Re}}\to \infty$ , as suggested by their asymptotic scalings. Under continuation, these RPOs tend to be increasingly dominated by large-scale, coherent vortices while weaker filamentary structures lessen. The vortex extraction algorithm largely avoids the filamentary components for class 1 states and the point vortex fit captures the large-scale flow features which are retained in the high- ${\textit{Re}}$ limit. Interestingly, these class 1 solutions contain examples with many interacting vortices (typical of high dissipation states seen in past work (Chandler & Kerswell Reference Chandler and Kerswell2013; Page et al. Reference Page, Holey, Brenner and Kerswell2024a ; Redfern et al. Reference Redfern, Lazer and Lucas2024)) but also examples of condensate-like structures dominated by a pair of opposite signed vortices, which are reminiscent of the Euler solutions discussed in previous work exploring Euler/Navier–Stokes connections (Kim & Okamoto Reference Kim and Okamoto2015; Kim et al. Reference Kim, Miyaji and Okamoto2017; Zhigunov & Grigoriev Reference Zhigunov and Grigoriev2023). The weaker performance of the point vortex fit in classes 2 and 3, and failure to extract the structure retained as ${\textit{Re}}$ increases, reflects the fact that – particularly in class 2 – the flow field is dominated by a large scale vortex pair and thin filamentary structures.