2. Objective functionals for unstable periodic orbits

The well-known incompressible Navier–Stokes equations, in non-dimensional form, describe the evolution of the 2-D flow velocity $(u,v)$ with time $t$,

(2.1)$$\begin{gather} \partial_t u + u \partial_x u + v \partial_y u ={-} \partial_x p + \frac{1}{Re} \Delta u + f_x, \end{gather}$$

(2.2)$$\begin{gather}\partial_t v + u \partial_x v + v \partial_y v ={-} \partial_y p + \frac{1}{Re} \Delta v + f_y, \end{gather}$$

subject to the incompressibility constraint

(2.3)\begin{equation} \partial_x u + \partial_y v = 0,\end{equation}

where the pressure $p$ acts as a Lagrange multiplier to enforce the incompressibility. Here $Re$ is the usual Reynolds number, obtained by non-dimensionalising the equation with respect to the forcing amplitude (note that different definitions of $Re$ are used for Kolmogorov flow in the literature), and $f_x$ and $f_y$ are the $x$- and $y$-components of a body force. Henceforth, following previous authors (Chandler & Kerswell Reference Chandler and Kerswell2013; Lucas & Yasuda Reference Lucas and Yasuda2021), we will take $f_x=\sin {4y}$, $f_y=0$, and consider a doubly periodic box of size $2{\rm \pi} \times 2{\rm \pi}$. Throughout, we take $Re=40$, which is sufficiently large to allow extensive chaotic behaviour, but permits the use of relatively low resolutions to expedite computations. The system (2.1)–(2.3) with this choice of forcing exhibits several symmetries which permit special exact coherent structures. Continuous invariance under the family of translations $x\mapsto x+\phi$ for $\phi \in \mathbb {R}$ allows travelling waves and relative periodic orbits (RPOs) with a drift velocity $c$ in this direction. The system is also invariant under the discrete transformation $(x,y)\mapsto (-x,y+{\rm \pi} /4)$, which permits RPOs, though these are equivalent to stationary periodic orbits with eight times the period.

By differentiating (2.1) and (2.2), we can derive a governing equation for the vorticity $\omega = \partial _x v - \partial _y u$,

(2.4)\begin{equation} \partial_t \omega +\partial_y\psi\partial_x\omega -\partial_x \psi\partial_y\omega = \frac{1}{Re}\Delta\omega +4\cos{4y}, \end{equation}

where the streamfunction $\psi$ is defined such that $u=\partial _y\psi$ and $v=-\partial _x\psi$ and so $\omega = - \Delta \psi.$

These two formulations, called respectively primitive variables (PV) and streamfunction vorticity (SV), are equivalent in two dimensions. The latter uses only two variables rather than three, since incompressibility is directly built into the formulation, but this advantage is lost when the equations are extended to three spatial dimensions. It is further possible to eliminate $\omega$ to give one equation solely in terms of $\psi$. As variational methods based on the resulting fourth-order equation were found to be very inefficient, we do not discuss this formulation further.

Care must be taken when comparing the boundary conditions between the two formulations. Assuming a doubly periodic streamfunction in the SV formulation implies that $u$ and $v$ have no mean component over the domain, i.e. no net flow, which is not a priori imposed by the PV formulation with periodic boundary conditions on $u$ and $v$. However, it is straightforward to show that (2.1)–(2.3) conserve the net flow, which is implicitly specified by the choice of initial conditions. Therefore, the system modelled by the SV formulation represents an invariant subspace of the PV system. For initial conditions with no net flow, both formulations are equivalent, and we henceforth make this choice.

Let us consider closed ‘loops’ in the phase space of the system, where each component ($u$, $v$ and $p$ or $\psi$ and $\omega$) is defined analogously to

(2.5)\begin{equation} u(x,y,s):[0,2{\rm \pi})\times[0,2{\rm \pi})\times[0,2{\rm \pi})\to\mathbb{R}, \end{equation}

periodic in all three dimensions, and assumed to be sufficiently smooth. Similarly, we can define loops of $v$, $p$, $\psi$ and $\omega$. Introducing a period $T>0$, such loops are exact periodic orbit solutions of the system if they satisfy the governing equations with $\partial _t = ({2{\rm \pi} }/{T})\partial _s$. In this perspective, $s$, which parameterises the loop in phase space, is a scaled time variable. From the governing equations of the two different formulations, we can define objective functionals based on the PV formulation

(2.6)\begin{align} J_{PV}[u,v,p,T] &= \frac{1}{2}\int \left\{ \left(\frac{2{\rm \pi}}{T}\partial_s u + u\partial_x u +v\partial_y u +\partial_x p -\frac{1}{Re}\Delta u - \sin{4y}\right)^{2}\right.\nonumber\\ &\quad \left.+\left(\frac{2{\rm \pi}}{T}\partial_s v + u\partial_x v +v\partial_y v +\partial_y p -\frac{1}{Re}\Delta v\right)^{2} +\left(\partial_x u + \partial_y v\right)^{2} \right\}\mathrm{d}V , \end{align}

and for the SV formulation,

(2.7)\begin{align} J_{SV}[\psi,\omega,T] &= \frac{1}{2}\int \left\{ \left(\frac{2{\rm \pi}}{T}\partial_s \omega + \partial_y\psi\partial_x\omega -\partial_x \psi\partial_y\omega - \frac{1}{Re}\Delta\omega -4\cos{4y}\right)^{2}\right.\nonumber\\ &\quad \left.\vphantom{\left(\frac{2{\rm \pi}}{T}\partial_s \omega + \partial_y\psi\partial_x\omega -\partial_x \psi\partial_y\omega - \frac{1}{Re}\Delta\omega -4\cos{4y}\right)^{2}}+\left(\omega+\Delta\psi\right)^{2} \right\}\mathrm{d}V , \end{align}

where the integrals are taken over $(x,y,s)\in [0,2{\rm \pi} )^{3}$. Both functionals are necessarily greater than or equal to zero, with equality if and only if the governing equations are satisfied everywhere on the loop. Therefore, finding zeros of these functionals gives us, in general, periodic orbits of the system. Equilibria are a special case, for which $\partial _s u = 0$ etc., and the value of $T$ is arbitrary. While there is an exact correspondence between zeros, i.e. global minima, of both objective functionals, local minima do not generally translate from one formulation to the other.

To minimise these objective functionals, we derive gradients. For the scalar variable $T$, these are simple partial derivatives, for example,

(2.8)\begin{equation} \frac{\partial J_{SV}}{\partial T} = \int -\frac{2{\rm \pi}}{T^{2}}\partial_s \omega\left(\frac{2{\rm \pi}}{T}\partial_s \omega + \partial_y\psi\partial_x\omega -\partial_x \psi\partial_y\omega - \frac{1}{Re}\Delta\omega -4\cos{4y}\right) \mathrm{d}V. \end{equation}

For the loop variables $u$, $v$, $p$, $\psi$ and $\omega$, we must instead derive variational derivatives which are defined at all $(x,y,s)\in [0,2{\rm \pi} )^{3}$, such as

(2.9)\begin{equation} \frac{\delta J_{SV}}{\delta \omega} ={-}\frac{2{\rm \pi}}{T}\partial_s I - \partial_x(I\partial_y \psi)+ \partial_y(I\partial_x \psi) - \frac{1}{Re}\Delta I + (\omega+\Delta\psi), \end{equation}

where for clarity, we have defined

(2.10)\begin{equation} I = \frac{2{\rm \pi}}{T}\partial_s \omega + \partial_y\psi\partial_x\omega -\partial_x \psi\partial_y\omega - \frac{1}{Re}\Delta\omega -4\cos 4y. \end{equation}

Since we have no physical boundaries in this problem, such variational derivatives are relatively straightforward to derive.

2.1. Relative periodic orbits

As mentioned above, there are symmetries in both the $x$ and $y$ directions of the system, which permit the existence of RPOs. We choose not to study RPOs which exploit the discrete symmetry in the $y$-direction, though this would be a simple extension of our method, where necessary considering two periods of the orbit to eliminate the discontinuity caused by the change of sign in the $x$-direction, as in Chandler & Kerswell (Reference Chandler and Kerswell2013).

For an RPO with a drift velocity $c$ in the $x$-direction relative to the frame in which there is no net flow, we can transform into a frame moving with the RPO, so that it becomes a simple periodic orbit. This is achieved by replacing $({2{\rm \pi} }/{T})\partial _s u$ by $({2{\rm \pi} }/{T})\partial _s u-c\partial _x u$ etc., so that extra terms appear in the gradients. We also now must optimise with respect to $c$, using the partial derivative

(2.11)\begin{align} \frac{\partial J_{SV}}{\partial c} = \int - \partial_x \omega\left(\frac{2{\rm \pi}}{T}\partial_s \omega-c\partial_x\omega + \partial_y\psi\partial_x\omega -\partial_x \psi\partial_y\omega - \frac{1}{Re}\Delta\omega -4\cos 4y\right)\mathrm{d}V \end{align}

and similarly for $J_{PV}$. Since a periodic orbit is a special case of an RPO with $c=0$, and travelling waves and equilibria are also captured, we use these extended objective functionals for all computations. The full objective functionals and all the gradients are given in the supplementary materials. It is important to note that the gradients for the PV formulation have no mean component, assuming that $u$ and $v$ do not, and so the (zero) net flow is conserved along gradients.

With the inclusion of the phase speed drift velocity $c$, the space of loops in each of the formulations is represented as a vector space of tuples

(2.12)\begin{equation} X = (u,v,p,T,c) \quad\mathrm{or}\quad X = (\omega,\psi, T,c), \end{equation}

with inner products defined on these in the obvious way,

(2.13)\begin{gather} \left\langle X_1,X_2\right\rangle_{PV} = \int \left(u_1u_2+v_1v_2+p_1p_2\right)\mathrm{d}V + T_1T_2 + c_1c_2,\end{gather}

(2.14)\begin{gather} \left\langle X_1,X_2\right\rangle_{SV} = \int \left(\omega_1\omega_2 + \psi_1\psi_2\right)\mathrm{d}V + T_1T_2 + c_1c_2.\end{gather}

3. Optimisation methods for minimizing the objective functionals

Minimising the value of the objective functionals we have defined, which are single valued functions of (at least in a discretised sense) a very large number of unknown variables with known derivatives, is best performed through gradient-based methods. This contrasts with methods by which we would directly find zeros of the integrands through Newton iteration and related methods, which require the same number of equations as unknowns. Descent methods find only local minima of the function, but to find exact solutions of the governing equations, we need to find zeros, i.e. global minima, of these non-negative functionals. The methods we present give monotonically decreasing and/or globally convergent results, but this does not mean we are guaranteed to converge on exact solutions to the governing equations.

The simplest method is gradient descent with a fixed step size $\epsilon$ at each iteration. Previous authors (Farazmand Reference Farazmand2016; Azimi et al. Reference Azimi, Ashtari and Schneider2022) have considered a dynamical system by the introduction of a fictitious time $\tau$ such that

(3.1)\begin{equation} \frac{{\rm d} u}{{\rm d} \tau} ={-}\frac{\delta J}{\delta u}. \end{equation}

Simple gradient descent is equivalent to discretising this system with a forward-Euler scheme, and some improvements may be found by considering terms implicitly or using higher-order schemes, or adaptive time steppers. This differs from the method of Lan & Cvitanović (Reference Lan and Cvitanović2004), whose fictitious system evolves in a direction equivalent to an infinitesimal step of a Newton–Raphson method, as opposed to in the steepest descent direction for the cost function. This offers much faster convergence, but requires the computation of an explicit Jacobian matrix.

However, we are not interested in solving the fictitious system accurately, but merely in finding its stable fixed points. Faster convergence is found using the nonlinear conjugate-gradient method (Hager & Zhang Reference Hager and Zhang2006), which requires only the first-order gradients. If the Hessian matrix were also calculated, Newton's method could be employed. We implemented this with the Fletcher & Reeves (Reference Fletcher and Reeves1964) choice of the conjugate-gradient parameter $\beta$. The optimal step size $\alpha$ was chosen by a line search, from an initial guess of $\alpha =10^{-5}$ for $J_{PV}$ and $\alpha =10^{-7}$ for $J_{SV}$. The strong Wolfe (Reference Wolfe1969) conditions refine this initial guess, using the tunable parameters $0< c_1< c_2<1$, where a smaller $c_2$ enforces a sufficiently large step that the objective functional decreases rapidly, while a larger $c_1$ enforces that we do not overshoot and take too large a step. Through trial and error, we found acceptable values for these parameters to be $c_1=10^{-5}$ and $c_2=0.999$, though we make no claims that these are optimal. With this choice of parameters, typically only one iteration on the line search algorithm was necessary before the Wolfe conditions were satisfied. The significantly smaller step sizes required for the conditions to be satisfied for SV mean that this method takes much longer to converge, as discussed in § 4.

Since $x$, $y$ and $s$ are all periodic, the fields can be easily expressed as Fourier series, using a pseudo-spectral approach – with 2/3 dealiasing – for the nonlinear terms. A new code was developed for this work in C$++$, with OpenMP parallelisation. Following Chandler & Kerswell (Reference Chandler and Kerswell2013), we perform time integrations using Heun's method, with Crank–Nicolson on the viscous terms at every substep. We concentrate on the specific case $Re=40$, which is sufficiently high to allow fully chaotic behaviour, but low enough so that a spatial resolution of $64\times 64$ grid points gives accurate results when compared with previous work. A resolution of $64$ grid points was also used in the temporal dimension, which was found to be more than sufficient for even the longest periodic orbits studied. A $2/3$ dealiasing rule was applied in the spatial but not temporal dimensions.

3.1. Leray projection

In § 2 we gave objective functionals for two different formulations of the system. The PV formulation does not assume the velocity field is divergence free, and instead includes the incompressibility condition within the objective functional, whereas the alternative SV formulation automatically enforces incompressibility. Since we know that any converged solution must be divergence free, it may also be desirable to enforce this during the optimisation within the PV formulation. However, care must be taken to ensure that the resulting method still ensures the monotonic decrease of the objective functional. For this, we recall the 2-D Leray projection operator (Temam Reference Temam2001)

(3.2)\begin{equation} \boldsymbol{P}:(u,v)\mapsto(u,v)-(\partial_x,\partial_y)\varDelta^{{-}1}(\partial_x u + \partial_y v), \end{equation}

where $\varDelta ^{-1}$ is the inverse Laplacian, here taken with periodic boundary conditions. For any $(u,v)$, observe that $\boldsymbol {P}(u,v)$ is divergence free. Then for the usual inner product,

(3.3) \begin{align} \left\langle\boldsymbol{P}(u,v),(u,v)\right\rangle&=\left\langle\boldsymbol{P}(u,v),\boldsymbol{P}(u,v)\right\rangle+\langle\boldsymbol{P}(u,v),\nabla\varDelta^{{-}1}(\partial_x y + \partial_y v)\rangle \nonumber\\ &=\left\|\boldsymbol{P}(u,v)\right\|^{2}-\langle\boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{P}(u,v),\varDelta^{{-}1}(\partial_x u + \partial_y v)\rangle\nonumber\\ &=\left\|\boldsymbol{P}(u,v)\right\|^{2} \geqslant 0, \end{align}

where we have used the fact that $\boldsymbol {\nabla }$ is adjoint to $-\boldsymbol {\nabla }\boldsymbol {\cdot }$ with respect to the inner product. Therefore, using an update direction of $-\boldsymbol {P}({\delta J_{PV}}/{\delta u},{\delta J_{PV}}/{\delta v})$ instead of $-({\delta J_{PV}}/{\delta u},{\delta J_{PV}}/{\delta v})$ will still converge to a minimum of the objective functional, since the projected gradient is never directed opposite to the unconstrained gradient.

This leaves us with three potential variational methods: optimising $J_{SV}$ (SV), optimising $J_{PV}$ (PV) or optimising $J_{PV}$ with Leray projection of the gradient at each step to enforce the divergence-free condition (PV-LP).

4. Comparison of methods

To compare the performance of the SV, PV and PV-LP methods, we discuss the convergence starting from two initial candidates, which were found through a recurrent flow analysis (see § 4.2). Figure 1 depicts a phase space projection of these two candidates as well as the converged solutions, and figure 2 shows the flow structure at single snapshots in time. The first candidate converged to a simple periodic orbit of period $T=5.38$, given as ‘P1’ in Chandler & Kerswell (Reference Chandler and Kerswell2013), and the second to the RPO ‘R19’ with $T=12.2$, although the algorithm has converged to two periods of this orbit, to give $T=24.4$. Note how vastly more complicated the initial guess is than the converged solution in this latter case, at least in this simple projection, hinting at the ability of the variational methods to converge to solutions from very distant initial guesses. The initial candidates were taken from a time series, with the 64 temporal collocation points for each loop evenly distributed along the series between the starting point and the nearest recurrence. The loop is then closed, which gives the appearance of a sharp discontinuity visible in figure 1. Since temporal derivatives are calculated after a Fourier transform, this discontinuity is not problematic.

Figure 1. Projections of the guess (thin) and converged (thick) loops for the two examples discussed. The energy input is defined as $I=8\int u\sin 4y\,\mathrm {d} x\,\mathrm {d} y/{\rm \pi} ^{2} Re$, and the dissipation as $D=8\int (|\boldsymbol {\nabla } u|^{2}+|\boldsymbol {\nabla } v|^{2})\,\mathrm {d} x\,\mathrm {d} y/{\rm \pi} ^{2} Re^{2}$, both of which have been normalised by the laminar solution (see Chandler & Kerswell Reference Chandler and Kerswell2013). The labels correspond to the snapshots shown in figure 2.

Figure 2. Slices at $s=0$ of the vorticity field for the guess (a,c) and converged (b,d) loops discussed in the text, at the snapshots labelled in figure 1. Movies are available in the supplementary materials available at https://doi.org/10.1017/jfm.2022.299.

Figure 3 shows the convergence of the relevant residual against the number of conjugate-gradient iterations, for the first of these two candidates. The algorithms were each run for three days of wall-clock compute time on a 28 core CPU. Due to the line search requiring more evaluations, each SV iteration is on average around 40 % slower than the two PV methods, which leads to fewer iterations in 72 h. More importantly, SV shows a significantly slower convergence rate. All three methods show an initial fast improvement of the objective function, followed by a much slower period, consistent with Azimi et al. (Reference Azimi, Ashtari and Schneider2022); the convergence rate in this latter region is much faster with the PV methods, as much smaller steps were required to be taken for SV. The precise reason for this is unclear, and we make no claim that our conjugate-gradient algorithm uses optimal parameters, but this result was robust after significant trial-and-error tweaking of them. This strongly suggests that the equation (3.1) of the ‘fictitious’ system is much stiffer for SV than PV. Whether Leray projection is used has negligible impact on the rate of convergence, but the path taken towards the solution is certainly different as PV allows for compressible intermediates. Consequently, the methods may converge to different solutions.

Figure 3. Convergence of the respective residual $J_{PV}$ or $J_{SV}$ for the PV formulation, with (PV-LP) and without (PV) projection and the SV methods applied to the first guess in figure 1. In the SV case, $J_{SV}$ starts at a significantly larger value, but also converges much more slowly.

To investigate the range of convergence, we systematically vary the initial guess by linearly interpolating between the converged solution and the candidate extracted from flow recurrences, $u_{initial} = (1-\gamma ) u_{solution} + \gamma u_{candidate}$, with $\gamma =0$ being the exact (relative) periodic orbit (P1 and R19) and $\gamma =1$ being the candidate. As shown on the top left in figure 4, for the first candidate, the PV formulation, both with and without Leray projection, shows good convergence over the full range of $\gamma$, with $J_{PV}$ consistently reaching less than $10^{-13}$. However, extrapolating to $\gamma >1$, where the initial guess is very different from the periodic orbit P1, the converged period $T$ is shown to rapidly vary. This is because, in this case, the algorithm is no longer converging to P1 but onto a travelling wave solution, for which the period is indeterminate. Note that around $\gamma =1.2$, the Leray-projected method does converge onto P1 whereas without it, it does not, which confirms that Leray projection does indeed have a significant role, and that without it, the PV formulation loses the divergence-free property of the guess in the course of converging before again becoming divergence free as the solution converged, a fact which was confirmed when the divergence was examined. In contrast with the PV formulations, the SV formulation performed poorly, with the residual not reaching an acceptably low level within the allotted time for all but the smallest values of $\gamma$, and the final value $T$ being noticeably wrong for $\gamma >0.3$, without the rapid fluctuations indicative of convergence onto a different solution.

Figure 4. Comparison of methods applied to an initial guess linearly interpolated/extrapolated between the converged periodic orbit ($\gamma =0$) and the recurrent flow analysis candidate ($\gamma =1$). The top plots show the final value of the residual $J_{PV}$, where a value of less than $10^{-10}$ is well converged. The lower plots show the final value of the period $T$. Left and right as per the candidates in figure 1.

For the second, longer period candidate, the results are similar, with convergence good for the PV formulations but much worse for the SV formulation. However, in this case, for $\gamma \gtrsim 1.3$, the PV formulations appear to converge to a consistent value of $T$ and a consistent value of $J_{PV}$, which is considerably greater than zero. This, therefore, is a local minimum of the objective functional. This demonstrates the importance of ensuring that the results of the algorithms we find are genuine solutions, rather than just local minima.

4.1. Comparison with Newton-based shooting

As the SV formulation is significantly slower and has worse convergence properties than the PV formulation, in this section we concentrate on the PV formulations, with (PV-LP) and without (PV) Leray projection and compare its performance to a traditional shooting method.

We compared against a Newton shooting (SN) method. This differs significantly from the loop-based methods discussed in this paper. A single point in phase space is integrated in time, using the SV formulation of the equations in a comoving frame, up to a guess of the period $T$, with a guess of the phase speed $c$. The SN residual is then given by the squared distance between the terminal and initial points, defined in our case, due to the implementation, as a slightly unusual metric on the Fourier coefficients of the vorticity. A Newton–Krylov solver attempts to minimise this so that the trajectory closes back on itself. Since the handling of time between the methods is very different, a numerically converged solution in SN does not automatically translate to a converged solution in PV, and the converged periods were observed to differ by up to ${\pm }0.001$. We expect these discrepancies to disappear as the resolutions are increased, though we did not test this. We used the Newton solver from Channelflow (channelflow.ch) (Gibson et al. Reference Gibson2021). With the Hookstep procedure of Channelflow, the algorithm is theoretically globally convergent, but in practice with a poor initial condition the steps can become so small that numerically the algorithm fails to converge, after a few hundred Newton iterations. The method was run for a wall-clock time of three days in each case.

It is instructive to apply SN to the two candidates discussed above. In the first case, the SN procedure described above does not converge starting from the initial guess. In the second case, it does converge, but converges to a simple travelling wave solution (‘T1’) rather than the RPO which the variational methods found.

We then consider the effects of a partial convergence of the variational methods before applying SN, a hybrid between the two. The procedure we follow is, after a fixed number of PV iterations, simply to take the first point in the loop as the initial condition for the shooting method, discarding information about all other points on the loop. Figure 5 shows the results of this on the two candidates discussed before. In both cases, after a number of PV iterations, SN converges rapidly to the same solution as the complete PV method. Once sufficiently many PV iterations have been performed for the result to be within the region of convergence for SN, performing more iterations did not appear to significantly improve the convergence time of SN. Note the differences between the two solutions: in the first case, many more SN iterations were required, but an order of magnitude fewer PV iterations than in the second case.

Figure 5. The PV method was paused at intervals and the SN method applied to the result. Top: the SN residual after varying numbers of PV iterations, the (differently defined) residual for which is also shown. Bottom: the number of SN iterations required for convergence to a residual of $10^{-5}$ after the PV method, for those which converged. Where no point is shown, the method did not converge. Left and right as per the candidates in figure 1.

The partially converged solutions from PV do not satisfy incompressibility, and so when converting to the SV representation used by the SN time stepper (by setting $\omega =\partial _xv-\partial _yu$ and $\psi =-\varDelta ^{-1}\omega$), some information is lost. This is not the case with the PV-LP formulation, which was found to modestly improve this hybrid approach, with initial SN residuals around one half of the equivalent PV case.

4.2. Recurrent flow analysis

To test the methods, we implemented a recurrent flow analysis to generate candidates for exact coherent solutions. Our implementation follows that of Chandler & Kerswell (Reference Chandler and Kerswell2013): at regular intervals in a time series, the current state at time $t$ is compared against previous snapshots at time $t-T$. The phase space distance, minimised over the phase velocity $c$ in the $x$-direction, is plotted against $t$ and $T$, and local minima below a threshold value are stored as candidate solutions. We made no attempt to search for RPOs which exhibit a shift in the $y$-direction, exploiting the discrete symmetry of the equations. There is no reason to suspect that including this would affect the methods we study. Subsequent authors (Lucas & Kerswell Reference Lucas and Kerswell2015; Page & Kerswell Reference Page and Kerswell2020) have improved upon this method to find more candidates, but the focus of the present work is merely to investigate the properties of the variational methods. A time series of length $200$ was searched every $0.1$ time unit for near recurrences up to period $T=50$, giving 106 candidate solutions for which the distance $\int (\omega (x,y,t)-\omega (x+cT,y,t-T))^{2} \,\mathrm {d} x\,\mathrm {d} y/\int \omega (x,y,t)^{2} \,\mathrm {d} x\,\mathrm {d} y$ was less than the threshold of $0.5$. This is a significantly greater threshold value than that used by Chandler & Kerswell (Reference Chandler and Kerswell2013).

Table 1 summarises the results from trying to converge each of these candidates via PV, PV-LP and SN. We counted as converged any solution for which the residual $J_{PV}$ was less than $10^{-8}$, as the known local minima all had $J_{PV}>10^{-6}$. For all cases except those with particularly long periods, in practice the final value of $J_{PV}$ was significantly smaller than $10^{-8}$. Each attempted convergence for each of the algorithms was given 72 h of wall-clock time on a 28 core CPU, and terminated as soon as the convergence criterion was met, and otherwise counted as a failed convergence. All three methods were successful at converging solutions, including several that had been found previously and a number of new solutions. Even with this relatively small sample, it is clear that our new variational methods converge to a greater variety of solutions within the chaotic attractor than the traditional shooting method. Instances were found where the Newton method converged for a given candidate and the variational methods did not, and vice versa, hinting at the complex, fractal regions of convergence that are believed to exist for the different methods.

Table 1. Results of different methods applied to candidates from the recurrent flow analysis. The names are from Chandler & Kerswell (Reference Chandler and Kerswell2013), for those solutions which are not new. All values are given to three significant figures. Note that the sign of $c$ is irrelevant, due to the symmetries of the system.

Both our variational methods and the shooting method tend to struggle more in converging orbits with longer periods compared with short ones, though the reasons for this differ. As the shooting method time marches the orbit, exponential error amplification leads to sensitive dependence of the recurrence condition on initial conditions and, thus, an ill-conditioned root search. It is for this reason that multi-shooting methods have been developed (Sánchez & Net Reference Sánchez and Net2010). In our variational methods, the loop structure means that recurrence is ‘built in’. Nevertheless, we observe that the algorithm converges more slowly for longer orbits. Indeed, a large number of candidates appeared to be converging to long-period orbits when the algorithms were terminated and may well have converged had they been allowed to continue. In other cases, definite local minima were found by the variational methods, at which point convergence is no longer possible. However, switching to a different variational method (e.g. from PV to PV-LP) at this point was found to subsequently lead to convergence to a solution.