3. Search for URPOs in 3-D plane Poiseuille flow

In this section, we present invariant solutions obtained by Newton searches initialised using SRDMD reconstructions for incompressible channel flow. The computational domain consists of two stationary parallel plates at wall-normal positions $y=\pm h$ , which are periodically extended in the streamwise and spanwise directions with periods $L_x$ and $L_z$ , respectively. The flow is driven by a constant mass flux and is characterised by the Reynolds number $\textit{Re}=U_0h/\nu$ based on the laminar centreline velocity $U_0$ and the kinematic viscosity $\nu$ . The flow is governed by the Navier–Stokes equations, which in their non-dimensional form read

(3.1) \begin{equation} \partial _t\boldsymbol{u}+(\boldsymbol{u}\boldsymbol{\cdot} \boldsymbol{\nabla} )\boldsymbol{u}=-\boldsymbol{\nabla} p+\dfrac {1}{\textit{Re}}\,\Delta \boldsymbol{u}, \end{equation}

together with the continuity equation

(3.2) \begin{equation} \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u}=0, \end{equation}

where $\boldsymbol{u}$ and $p$ are the velocity and pressure, respectively. The Navier–Stokes equations are subject to no-slip and no-penetration boundary conditions at the walls,

(3.3) \begin{equation} \boldsymbol{u}(x,y=\pm h,z,t) = \boldsymbol{0}, \end{equation}

and periodic boundary conditions in $x$ and $z$ . This system admits the one-dimensional laminar velocity $\boldsymbol{u}_l=U(y)\,\hat{\boldsymbol {x}}$ , with the parabolic profile

(3.4) \begin{equation} U(y) = 1-y^2. \end{equation}

This profile corresponds to the choices $h=1$ and $U_0=1$ , which are used throughout this paper. We study the flow in a channel with a streamwise length $L_x=2\unicode{x03C0}$ and a spanwise width $L_z=0.77\unicode{x03C0}$ . For numerical integration of the Navier–Stokes equations and Newton searches for invariant solutions, we used the channelflow 2.0 library (Gibson Reference Gibson2014; Gibson et al. Reference Gibson2019). Numerical details for this study are summarised in table 1.

Table 1. Numerical details. The Reynolds number is denoted by $\textit{Re}$ , and the friction Reynolds number by $\textit{Re}_{\tau }$ . The streamwise, spanwise and wall-normal dimensions of the computational domain are denoted by $L_x$ , $L_z$ and $2h$ , respectively. The number of grid points in the streamwise, spanwise and wall-normal directions are denoted by $N_x$ , $N_z$ and $N_y$ , respectively. The test case is adopted from the work of Rawat et al. (Reference Rawat, Cossu and Rincon2014).

The DNS trajectories were produced using an initial flow field based on the work of Rawat et al. (Reference Rawat, Cossu and Rincon2014). The initial flow field consists of the laminar profile (3.4), perturbed by streamwise counter-rotating vortices of amplitude $A_1$ and a sinusoidal perturbation of amplitude $A_2$ in the spanwise direction:

(3.5) \begin{equation} \boldsymbol{u}_0=\{U(y),0,0\}+A_1\left\{0,\frac {\partial \psi }{\partial z},-\frac {\partial \psi }{\partial y}\right\}+A_2\{0,0,W_{\textit{sin}}\}, \end{equation}

where

(3.6) \begin{equation} \psi (y,z)=(1-y^2)\sin {(\unicode{x03C0} y)}\sin {(2\unicode{x03C0} z/L_z)} \end{equation}

and

(3.7) \begin{equation} W_{\textit{sin}}(x,y)=(1-y^2)\sin {(2\unicode{x03C0} x/L_x)}. \end{equation}

The initial flow field (3.5) has the discrete symmetries

(3.8) \begin{align}&\quad\,\,\, [u,v,w](x,-y,z)=[u,-v,w](x,y,z), \end{align}

(3.9) \begin{align}& [u,v,w](x+L_x/2,y,-z)=[u,v,-w](x,y,z), \end{align}

(3.10) \begin{align}& [u,v,w](x+L_x/2,-y,z)=[u,-v,w](x,y,z), \end{align}

which were explicitly imposed during the simulations. Although the governing equations and boundary conditions are equivariant under continuous translations in either of the periodic wall-parallel directions, the imposed discrete symmetries prevent the solution from drifting in the spanwise direction $z$ . However, the flow is not constrained in the streamwise direction $x$ . The flow’s continuous symmetry in this direction allows for time-dependent drift, leading to poor performance of the DMD method, and necessitating the reduction of the continuous symmetry. In order to reduce the continuous symmetry, we applied the first Fourier mode slicing technique described in § 2.3. To ensure that all involved scalar products in (2.17) and (2.18) are well defined at all times when applying the slicing technique, we utilise the symmetries (3.8)–(3.10) for all velocity components and their relation to the symmetries of the template functions (Engel et al. Reference Engel, Ashtari and Linkmann2023).

Figure 2. The $L_2$ -norm of the velocity as a function of time for the DNS of the chaotic trajectory (blue) and the trajectory that shadows the URPO for approximately four cycles between $t\approx 1000$ and $4000$ (orange). Both DNS correspond to Reynolds number $\textit{Re}=3000$ .

3.1. Method validation

To validate SRDMD, we applied both slicing and DMD to the URPO in the channel flow identified by Rawat et al. (Reference Rawat, Cossu and Rincon2014). This solution is obtained with the friction Reynolds number $\textit{Re}_{\tau }\approx 151$ ( $\textit{Re}=3000$ ), and has period $T=739$ . We used the same computational domain size and resolution as Rawat et al. (Reference Rawat, Cossu and Rincon2014) (see table 1). With one unstable direction in the invariant subspace formed by the discrete symmetries (3.8)–(3.10), this solution serves as an ideal candidate for validating the methods, since it can be identified using a bisection method, as described in Rawat et al. (Reference Rawat, Cossu and Rincon2014). The bisection method iteratively adjusts the perturbation amplitudes $A_1$ and $A_2$ in the initial condition (3.5) to generate progressively longer trajectories that neither become turbulent nor decay to the laminar state. This is accomplished by performing a bisection on the amplitudes from two preceding bisection steps: one that leads to decay to the laminar state, and the other that leads to turbulent behaviour after a sufficiently long time. In this way, we tracked a trajectory shadowing the URPO of interest for approximately four cycles. Figure 2 illustrates the evolution of the velocity $L_2$ -norm for both a chaotic trajectory and the trajectory shadowing the URPO. Inspection of the data reveals that the corresponding flow field shows a streamwise drift, since the plane Poiseuille flow has a continuous translation symmetry in the streamwise direction. However, the system is prevented from drifting in the spanwise direction through the imposed shift–reflect symmetry (3.9).

To validate the methodology introduced in § 2, we first applied the slicing technique to reduce the drift from the flow. We used a slicing template function $\hat{\boldsymbol{u}}^{\prime}_x$ aligned in the spanwise direction $z$ , with a profile function $f(y)=2.5\,T_2(y)-1.25\,T_3(y)-1.25\,T_5(y)$ (see § 2.3). Our choice of alignment direction and Chebyshev polynomials was based on the flow field’s discrete symmetries, ensuring that the slice velocity (2.19) is well defined at all times. Figure 3 shows the slicing phase velocity obtained via reconstruction equation (2.19) for the last 3000 advective time units. We observe that the slicing phase velocity is well defined over the whole time window. A failure of the slicing method at a specific moment in time would be indicated by discontinuities in $\dot {\phi }_x(t)$ , resulting in non-physical jumps of the flow field through the computational domain in the streamwise direction.

Figure 3. Phase velocity for the global drift in the streamwise direction as a function of time, obtained by applying the first Fourier mode slicing technique to the DNS dataset. Snapshots were sampled with time resolution $\Delta t=1$ . The flow does not have a global drift in the spanwise direction, as it is constrained by a discrete symmetry.

Figure 4. State-space projections of the trajectory following a URPO over one period $T_p\approx742$ . The state space is projected onto the first three dominant POD modes, with $a_1, a_2, a_3$ representing the corresponding amplitudes. (a) Trajectory of the original dataset. (b) The symmetry-reduced trajectory after applying the first Fourier mode slicing technique. The black dots mark the start and end points of the trajectories.

The result of the symmetry reduction is best illustrated by inspecting a 3-D state-space projection of the trajectory. Figure 4(a) presents the trajectory before symmetry reduction, and figure 4(b) after symmetry reduction, following the URPO for one of its cycles. The flow fields are projected onto the first three dominant modes of proper orthogonal decomposition (POD). The POD modes in the non-symmetry-reduced case were calculated for a data matrix corresponding to a segment of the original trajectory following the URPO for one full cycle. Columns of the data matrix contain the grid values of instantaneous velocity fields in the physical domain, separated by a $\Delta t = 1$ advection time unit. For the symmetry-reduced case, we applied the same procedure, but the POD modes were calculated from the symmetry-reduced segment of the flow fields. Clearly visible is the more complicated state-space trajectory for the DNS with no post-processing in comparison to the symmetry-reduced case, due to the additional degree of freedom. In contrast to the DNS without post-processing, the start and end points of the symmetry-reduced trajectory (indicated by black dots) coincide after a time period that matches the URPO period reported in Rawat et al. (Reference Rawat, Cossu and Rincon2014).

After ensuring the success of the symmetry reduction, we proceeded with the DMD analysis of the DNS and symmetry-reduced data. All parameters used for producing the period guess and the low-dimensional flow field representation in this case were determined through experimentation. This includes the observation time window $T_w=680$ advective time units as well as the snapshot sampling separation $\Delta t$ and spacing $\delta t$ of 17 advective time units. The sparsity-promoting SRDMD algorithm utilises a bisection method to determine the penalty parameter $\gamma$ – calculated as 12.5 in this case – ensuring a low-dimensional reconstruction by selecting five optimal DMD modes. This choice is based on the study of Page & Kerswell (Reference Page and Kerswell2020), where five DMD modes were found to be sufficient to construct initial guesses for UPOs that converge successfully in a Newton search. Figure 5 compares the DMD spectra obtained by applying sparsity-promoting DMD to the DNS and symmetry-reduced data. The red line indicates the unit circle around which the eigenvalues are distributed. Eigenvalues selected by the sparsity-promoting DMD are highlighted in blue.

Figure 5. Comparison of (a) the classical DMD spectra and (b) the SRDMD spectra for the trajectory shadowing the oscillating structure. The observation window when applying the SRDMD was set to be 680 advective time units, which represents approximately $92\,\%$ of the converged period. The snapshot separation was set to 17 advective time units. The red line represents the unit circle. The eigenvalues selected by the sparsity constraint for the reconstruction are highlighted with filled markers.

For both the original and the symmetry-reduced data, the neutral DMD mode corresponding to the temporal mean is selected, as well as the DMD mode corresponding to the fundamental frequency. In sparsity-promoting DMD without symmetry reduction, the third pair of eigenvalues lies off the unit circle and relates to fast and small-scale dynamics as its imaginary part is large in absolute value. In contrast, the sparsity-promoting SRDMD selects the third DMD mode corresponding to the first higher harmonic of the fundamental frequency. This is plausible because the processed trajectory, which closely shadows the URPO, already represents a low-dimensional flow. The symmetry reduction enabled the DMD-based analysis to capture this characteristic as expected. By selecting the fundamental frequency and its higher harmonics, both the fundamental frequency and its residual can be calculated, as given in (2.11) and (2.12), yielding $T_g=741.979$ and $\epsilon _{\omega }\approx 10^{-9}$ .

The low-dimensional reconstruction of the flow (2.7) provides a closed curve in state space, which is not a solution trajectory of the governing equations but approximates the DNS data over the processed time window. Any time section along this curve can serve as an initial condition for the Newton solver (see figure 1). However, we followed Page & Kerswell (Reference Page and Kerswell2020) and initialised the Newton searches using the reconstructed flow field closest to the original trajectory, considering the $L_2$ -norm as the metric. The drift in the streamwise direction, which needs to be taken into account when applying a Newton search, can be derived from the slice phase tracked during the symmetry reduction by

(3.11) \begin{equation} \delta x = \frac {\phi _{t_0} - \phi _{t_p}}{2\unicode{x03C0} }, \end{equation}

where $t_0$ refers to the starting time of the SRDMD observation time window, and $t_p = t_0 + T_g$ corresponds to the time instant after one estimated period. The Newton search converges after seven iterations, reaching a residual $\mathcal{O}(10^{-14})$ , and yielding the predicted URPO with period $T_{\textit{URPO}}=742.017$ . The period of the converged URPO matches the period of the URPO reported in Rawat et al. (Reference Rawat, Cossu and Rincon2014) by a relative error $0.4 \,\%$ .

Figure 6. The $L_2$ -norm of the velocity along a trajectory at the friction Reynolds number $\textit{Re}_{\tau }\approx 51$ as a function of time. The time series shows chaotic behaviour for approximately 2600 advective time units.

3.2. Search for previously unknown URPOs

In this subsection, we present previously unknown URPOs buried in the chaotic repeller of channel flow at friction Reynolds number $\textit{Re}_{\tau }\approx 51$ (Zammert & Eckhardt Reference Zammert and Eckhardt2015). The choice of the Reynolds number is motivated by producing a weakly turbulent but long enough trajectory with a state space feasible for a SRDMD study. The discrete symmetry (3.9) breaks the continuous translational symmetry in the spanwise direction, preventing the solution from drifting in that direction. As a result, the present study focuses on handling the drifts due to the mean advection in the streamwise direction only. The spatial resolution of the simulation was increased from the test case in § 3.1 to $(N_x,N_y,N_z)=(48,65,48)$ . We tested the existence of all obtained solutions with a higher resolution $(N_x,N_y,N_z)=(60,81,60)$ , and found changes in their period to be below relative error $0.01\,\%$ .

We obtained sufficiently long, weakly chaotic trajectories by applying the same bisection method used in § 3.1 to an initial condition according to (3.5). This procedure was necessary, since initial conditions of random perturbations can quickly laminarise. Figure 6 shows the $L_2$ -norm of the simulated velocity fields as a function of time, exhibiting a chaotic character for approximately 2600 advective time units. Figure 7 shows velocity contour plots of the snapshot taken from the DNS at time $t=1000$ . The colour code indicates the deviation of the streamwise velocity from the laminar profile (3.4). Figure 7(a) shows the 3-D flow field, and figure 7(b) shows the flow field averaged in the streamwise direction. The black lines indicate the streamwise-averaged $y$ – $z$ streamfunction. The flow is dominated by a low-momentum region around its centre, which extends in the streamwise direction and is accompanied by four near-wall high-momentum streamwise streaks. Similar to the test case, the flow drifts fast in the streamwise direction, while it is prohibited from drifting in the spanwise direction due to the imposed shift–reflect symmetry (3.9). To remove the streamwise drift from the flow, we followed our procedure from the test case, using the same template profile.

Figure 7. Velocity contour plots of the flow field taken from the DNS at $t=1000$ advective time units (see figure 6). The colour code indicates the deviation of the streamwise velocity from the laminar profile. (a) The velocity contour in a 3-D plot. (b) The velocity field averaged in the streamwise direction. The black lines indicate the streamwise-averaged $y$ – $z$ streamfunction.

Figure 8. Phase velocity for the global drift in the streamwise direction as a function of time for the chaotic trajectory over a duration of 2600 advective time units. The zoomed inset highlights a peak in the phase velocity, illustrating that the projection onto the symmetry-reduced slice remains well defined at this point. In contrast, the five dominant peaks correspond to singularities in the phase velocity, occurring when the denominator in the reconstruction equation (2.19) vanishes as the trajectory approaches the slice border.

Figure 8 shows the slice phase velocity in the streamwise direction as a function of time. The inset magnifies a peak in the slice phase velocity, demonstrating that the projection onto the symmetry-reduced slice remains well defined at this point. In contrast, the five dominant peaks correspond to singularities in the phase velocity, which occur when the denominator in the reconstruction equation (2.19) vanishes as the trajectory approaches the slice border. However, these peaks are rare cases. Therefore, the template choice is sufficient, as we can reject initial guesses from SRDMD whose observation windows include a change of sign in the denominator of (2.19), resulting in a singularity in the phase velocity.

The DMD analysis was performed by scanning the symmetry-reduced flow with observation time window sizes $T_w\in \{60,80,100,120\}$ advective time units. The scan started at $t=0$ , with each time window shifted in steps of one unit to $t=2600$ , and snapshot spacing $\Delta t=3$ was used. To ensure the convergence of the DMD algorithm, we focused on the distribution of the resulting eigenvalues, expecting them to lie on the unit circle in the complex plane when the observation windows shadow a URPO, as these are stationary on average. In practice, a range of parameter values typically satisfies this criterion. By setting the residual acceptance threshold for the estimated fundamental frequency to $\epsilon _\omega \lt 10^{-3}$ , we constructed 354 initial guesses, each comprising a flow field, an estimated period, and the shift information from the slicing analysis. In the subsequent Newton search, 35 out of the 354 initial conditions converged to a URPO, resulting in a success rate of approximately 10 %. Among these 35 cases, we identified six distinct URPOs, which were found at several positions in the time series, suggesting their importance for the flow dynamics. Table 2 summarises the properties of the distinct URPOs found within the chaotic flow.

Table 2. Distinct URPOs converged using SRDMD to estimate periods and construct initial flow fields for Newton searches. We denote by $t_{i}$ the time at which the $L_2$ -distance between the DNS and the SRDMD reconstruction is minimal. The reconstructed flow field at time $t_{i}$ was used as the initial flow field for the Newton search. The length of the time window in the respective SRDMD analysis is denoted by $T_w$ . The estimated period is $T_g$ , and $T_{\textit{URPO}}$ is the period of the converged URPO. We give the deviation of $T_g$ from $T_{\textit{URPO}}$ in per cent. The number of distinct initial data points resulting in the same converged URPO is denoted by $N_{\textit{Hits}}$ . Also shown are $N$ , the number of unstable directions, $\sum _{j=1}^N\text{Re} {(\lambda _j)}$ , the sum of the real parts of all the unstable eigenvalues (growth rates), and $\max \,\text{Re} {(\lambda _j)}$ , the real part of the eigenvalue corresponding to the dominant unstable direction.

Figure 9 shows a 2-D state-space projection of the URPOs listed in table 2. The projection shows the total energy input $I$ against the total dissipation $D$ , normalised by their respective laminar values $I_0$ and $D_0$ . The grey shaded area represents the probability density function (p.d.f.) of the dynamics obtained from the DNS, where density decrease is represented by a grey scale from black to white. The majority of URPOs cover the dense part of the p.d.f. and therefore represent the typical dynamics. This is important since we aim to provide appropriate predictions of flow properties based on those of the identified URPOs. The URPO with the shortest period $T\approx 67$ is embedded in a low-probability region of low energy input and dissipation, and is only barely visited by the chaotic trajectory.

Figure 9. The 2-D state-space projection of the identified URPOs with respect to the total energy input $I$ and dissipation $D$ , normalised by their laminar values $I_0$ and $D_0$ , respectively. The dashed line represents the diagonal in the $D$ – $I$ plane, indicating states where energy input $I$ and dissipation $D$ are exactly balanced. The grey shaded area corresponds to the p.d.f. of the dynamics, where the density decrease is represented by a grey scale from black to white.

Figure 10. Recurrent flow analysis of the DNS. The greyscale colour code represents the $L_2$ -distance of DNS snapshot pairs separated by the time interval $T$ . Dark regions reflect close distances, suggesting possible recurrence events. Black circles show URPOs converged using the sparsity-promoting SRDMD method as a function of time $t$ , when initial flow fields were passed to the Newton solver, and their converged periods $T$ . Red-filled black circles indicate a selection of distinct URPOs from table 2. Multiple guesses have successfully converged to a URPO despite the absence of a recurrence signal, highlighting the complementary role of the DMD approach, which eliminates the need to track the entire URPO over a complete cycle.

Figure 10 shows a classical recurrent flow analysis on the symmetry-reduced dataset (see § 2.2). The grey scale represents the $L_2$ -distance between future flow states and previous ones. Possible recurrences are indicated in dark regions as a function of simulation time on the horizontal axis and possible periods on the vertical axis. The scattered circles represent URPOs converged using the sparsity-promoting SRDMD method. Markers filled in red correspond to distinct URPOs listed in table 2. Each marker shows the period of the respective converged solution and the time location of the initial flow field used for the Newton search, taken from the SRDMD reconstruction of the DNS. At several points in time, URPOs with periods $T_{\textit{URPO}}\in \{86.08,93.86\}$ are identified in the absence of a recurrence event. This highlights the complementary role of symmetry-reduced DMD in the search for URPOs compared to recurrent flow analysis, as a URPO would likely have gone undetected in these cases. The complementary role of the DMD-based method is further supported by the results of Newton searches initialised using guesses from recurrent flow analysis. Our guesses from recurrent flow analysis consist of the velocity field at time $t$ , and the period $T$ corresponding to the $(t,T)$ coordinate where $R$ takes a local minimum value. The initial guess for the spatial shift is taken as the difference between the shifts applied to $\boldsymbol{u}(t)$ and $\boldsymbol{u}(t+T)$ to bring them to the slice (see (2.16)). We selected the local minima where $R_{\textit{min}}\lt 0.2$ , yielding twelve guesses, of which two successfully converged (not shown in figure 10). Both successful searches converged to the URPO with period $T_{\textit{URPO}}=86.08$ , to which the DMD-based method converged 27 times. Therefore, guesses from recurrent flow analysis do not result in any additional URPOs compared to the DMD-based approach for the analysed time series.

In what follows, we connect the set of successfully predicted and converged invariant solutions to the physical properties of the flow.

4. The URPOs and statistical properties of the flow

In order to use the set of URPOs to predict the statistical properties of the flow, we investigated two different periodic orbit averaging approaches. First, we utilised POT and followed Cvitanović et al. (Reference Cvitanović, Artuso, Mainieri, Tanner and Vattay2016) by applying a simplified averaging equation for bounded flows:

(4.1) \begin{equation} \Gamma _{\textit{Prediction}}:=\frac {\langle \Gamma \rangle }{\langle T\rangle }, \end{equation}

where

(4.2) \begin{equation} \langle T\rangle :=\sum _{\pi }(-1)^{k+1}\frac {\sum _{i=1}^kT_{p_i}}{\Lambda _{p_1}\Lambda _{p_2}\cdots \Lambda _{p_k}} \end{equation}

is the mean cycle period, and

(4.3) \begin{equation} \langle \Gamma \rangle :=\sum _{\unicode{x03C0} }(-1)^{k+1}\frac {\sum _{i=1}^{k}T_{p_i}\Gamma _{p_i}}{\Lambda _{p_1}\Lambda _{p_2}\cdots \Lambda _{p_k}}. \end{equation}

Here, we sum over so-called pseudo-cycles $\unicode{x03C0} =(p_1,p_2,\ldots ,p_k)$ , each consisting of combinations of prime cycles $p_{i}$ . For each $i{\text{th}}$ URPO, the quantity $\Gamma _{i}$ represents the temporal average of the desired physical quantity, such as the mean flow, over one cycle. Additionally, (4.2) and (4.3) incorporate the URPOs’ periods $T_{i}$ and their corresponding Floquet multipliers

(4.4) \begin{equation} \Lambda _k^{(i)}=\exp {\left (\text{Re} {\left(\lambda _k^{(i)}\right)}\,T_{\textrm{i}}\right )}, \end{equation}

where we take into account unstable directions only. As a truncation criterion for the sums in (4.2) and (4.3), Cvitanović et al. (Reference Cvitanović, Artuso, Mainieri, Tanner and Vattay2016) suggest using the pseudo-cycle period and the stability ordering of the pseudo-cycle expansion. The latter is easily feasible using

(4.5) \begin{equation} \lvert \Lambda _{p_1}\cdots \Lambda _{p_k}\rvert \leqslant \Lambda _{\textit{max}} \end{equation}

as the truncation criterion, ignoring pseudo-cycles of greater instability than $\Lambda _{\textit{max}}$ . In the second approach, we followed Chandler & Kerswell (Reference Chandler and Kerswell2013) and applied the averaging formula

(4.6) \begin{equation} \Gamma _{\textit{Prediction}}:=\frac {\sum _{i=1}^N\omega _{i}\Gamma _{i}}{\sum _{i=1}^N\omega _{i}}, \end{equation}

where the weights

(4.7) \begin{equation} \omega _i=\frac {T_i}{\sum _{k\in \mathcal{K}_{i}}\text{Re} {\left(\lambda _k^{(i)}\right)}} \end{equation}

incorporate the period $T_{i}$ and positive real parts of the $i{\text{th}}$ URPO’s spectra. Though purely heuristic in nature, this approach, referred to as escape-time weighting, is intuitively understandable. First, it states that chaotic trajectories should spend less time in the vicinity of more unstable URPOs. Second, it assigns greater weight to URPOs with longer periods, as these occupy a larger proportion of the state space compared to those with shorter periods.

Figure 11. (a) Mean streamwise velocity profile $\langle u\rangle$ , and (b) Reynolds stress $\langle uu\rangle$ , both for the DNS (grey line) in comparison with the predictions using POT (blue) and the escape-time weighting (orange).

In the following, we present predictions of statistical flow properties provided by both approaches. We start with the mean streamwise velocity profile $\langle u\rangle$ and Reynolds stresses $\langle u_{i} u_j\rangle$ , where the angle brackets indicate averaging in the $x$ and $z$ directions as well as in time. These quantities, as evaluated from the DNS along with their predictions, are shown in figures 11 and 12. Both the POT and the escape-time weighting show promising predictions for each quantity. In a direct comparison, the escape-time weighting outperforms the POT approach. To explain why POT provides less accurate predictions than the escape-time weighting, we follow the argument of Chandler & Kerswell (Reference Chandler and Kerswell2013) and refer to a very strong variation of the amplitudes $\Lambda _{i}$ used in the pseudo-cycle expansion, which can differ by a factor of $10^4$ ; in contrast, the terms $\sum _{k\in \mathcal{K}_i}\text{Re} {(\lambda _k^{(i)})}$ used in the escape-time weighting lie much closer together. Furthermore, the accuracy of all predictions is restricted by the number of URPOs present in this study. More challenging is the p.d.f. prediction for energy input $I$ and dissipation $D$ . Figure 13 shows the p.d.f.s of $I$ and $D$ , normalised by their respective laminar values, for the DNS as well as their predictions using POT and escape-time weighting. We observe that the p.d.f. predictions for $I/I_0$ and $D/D_0$ exhibit shapes similar to those obtained from the DNS. However, they over-predict high energy input ( $I/I_0 \gt 2.4$ ) and dissipation ( $D/D_0 \gt 2.4$ ) regions, while under-predicting the mid-range energy input ( $1.8 \lt I/I_0 \lt 2.1$ ) and dissipation ( $1.8 \lt D/D_0 \lt 2.1$ ) regions. This discrepancy arises from the absence of URPOs that correspond to these specific energy input and dissipation ranges, as illustrated in figure 9. Similarly, the predictions based on the single URPO with $T_p = 67.02$ in the low energy input and dissipation regions fail to accurately capture the p.d.f. tails in those low-value ranges, as expected. Notably, the p.d.f. predictions based on POT and escape-time weighting show qualitatively similar results, reinforcing the consistency of these heuristic approaches.

Figure 12. (a) Reynolds stress $\langle uv\rangle$ and (b) Reynolds stress $\langle vv\rangle$ , both for the DNS (grey line) in comparison with the predictions using POT (blue) and the escape-time weighting (orange).

Figure 13. P.d.f.s for normalised (a) energy input $I/I_0$ and (b) dissipation $D/D_0$ . The solid grey line represents the p.d.f. for the DNS. The dash-dotted lines represent the prediction using POT in blue and the escape-time weighting in orange. The p.d.f.s were calculated via kernel density estimation.

5. Enabling multi-shooting to find URPOs

The state-space complexity of turbulent flow makes finding URPOs with the Newton method challenging, as initial conditions $\{\boldsymbol{u}_0, T\}$ that are not close to the expected solution diverge further when integrated over time. A multi-shooting method can stabilise initial deviations by dividing the time interval $[0,T]$ into smaller intervals, and therefore leads to improvements over the single-shot method (Sánchez & Net Reference Sánchez and Net2010). In this context, SRDMD provides a complementary approach by producing an entire object shadowing the suspected URPO, even in the absence of recurrences. Rather than choosing a single-shot method and taking as the initial guess the reconstructed flow field closest to the DNS trajectory with respect to a specific norm, we enable multi-shooting by incorporating reconstructed flow fields that are potentially not covered by the processed trajectory segment. Figure 14 shows a state-space sketch illustrating multi-shooting with three shots. The grey solid line represents a segment of the symmetry-reduced trajectory shadowing a symmetry-reduced URPO (not shown), whose period guess is $T$ . The blue dashed line represents the SRDMD reconstruction of the global object incorporating initial guesses that lie beyond the trajectory. The dashed grey line represents the trajectory resulting from time-forward integration $f^T$ of the initial flow field $\boldsymbol{u}_0$ . Its deviation from the symmetry-reduced trajectory represents the effect of the system’s underlying continuous symmetry. The pull-back operator $\sigma (t)$ maps between the original and the symmetry-reduced trajectories. For an $M$ -shot search, the pull-back operator $\sigma _m$ , associated with the $m$ th flow field, represents the spatial shift corresponding to the phase difference given by

(5.1) \begin{equation} \Delta \phi _m:=\int _0^{mT/M}\dot {\phi }(t)\,\mathrm{d}t,\quad m=0,1,2,\ldots ,M. \end{equation}

The relevant phase velocity $\dot {\phi }$ is given by (2.19) or (2.20), and results from symmetry reduction using the first Fourier slicing technique. Suppose that $\boldsymbol{u}_0,\boldsymbol{u}_1\ldots ,\boldsymbol{u}_{M-1}$ are $M$ snapshots from the reconstructed trajectory (2.7), equally spaced over the time interval $t\in [0,T]$ . The multiple-shot Newton method solves the root-finding problem

(5.2) \begin{equation} \begin{cases} f^{T/M}(\boldsymbol{v}_0) - \boldsymbol{v}_1= \boldsymbol{0}, \\ f^{T/M}(\boldsymbol{v}_1) - \boldsymbol{v}_2= \boldsymbol{0}, \\ \quad \vdots \\ f^{T/M}(\boldsymbol{v}_{M-1}) - \sigma _M\boldsymbol{v}_0= \boldsymbol{0}, \end{cases} \end{equation}

Figure 14. Initialisation of a multiple-shot Newton search using the SRDMD technique. The schematic sketch illustrates the state space, where a symmetry-reduced DNS trajectory segment (solid grey line) serves as the observation window for the SRDMD method. This method reconstructs the trajectory as a periodic function of time, appearing as a closed loop on the slice (dashed blue line). Integration of the reconstruction equation along the loop (5.1) determines the relative phase of each snapshot with respect to a reference snapshot, denoted as $\boldsymbol{u}_0$ . This provides an approximation of the URPO in the full state space (dashed grey line), which does not close on itself. Filled markers represent three snapshots equally spaced in time along the reconstructed trajectory. These snapshots, when transformed to the unconstrained state space (empty markers), are used to initialise a multiple-shot Newton search. For an exact URPO, these snapshots lie on a single integral curve $f^t(\boldsymbol{u}_0)$ of the vector field induced by the Navier–Stokes equations.

where $\boldsymbol{v}_m:=\sigma _m\boldsymbol{u}_m,\ m=0,1,\ldots ,M-1$ . The unknowns in this search include $\{\boldsymbol{v}_0,\ldots ,\boldsymbol{v}_{M-1},\sigma _M,T\}$ . For this vector of unknowns, a closed system of nonlinear equations is obtained by augmenting (5.2) with appropriate phase-locking constraints (Viswanath Reference Viswanath2009).

We demonstrate our approach with an example using a two-shot Newton search. We consider a guess identified by the SRDMD method described in § 2, with guessed period $T_g=62.1673$ . The standard single-shot Newton search failed to converge when initialised with snapshots corresponding to different phases along the reconstructed periodic function $\boldsymbol{u}(t)$ , defined by (2.7). We did a two-shot search using snapshots at $t=0$ and $t=T_g/2$ from the reconstructed flow. The choice of the snapshots is not unique, and any pair of snapshots separated in time by $T_g/2$ can be used for a two-shot search, as (2.7) provides an explicit continuous function of time. The two-shot Newton search converges within $28$ iterations to an exact solution. Figure 15 shows the convergence result. In this plot, $G$ represents the to-be-zeroed vector in the respective root-finding problem, and $x$ represents the vector of unknowns in the discretised problem. We consider an exact solution to be achieved within machine accuracy if the $L_2$ -norm of $G(x)$ falls below $10^{-13}$ . The single-shot search starting from the snapshot at $t=0$ could not converge within $200$ iterations, after which the search was terminated. The single-shot search starting from the snapshot at $t=T_g/2$ stopped after $117$ iterations because the hook-step trust-region method could no longer improve the residual. These two failed searches are shown too in figure 15, where the residual stagnates after approximately $50$ Newton iterations. To save computational time, the multiple-shot Newton search was done at a lower resolution $(N_x,N_y,N_z)=(16,41,16)$ (convergence shown in figure 15). The converged URPO was then continued to the original resolution $(N_x,N_y,N_z)=(48,65,48)$ to confirm that it is structurally stable and to determine its accurate period $T=67.0226$ .

Figure 15. Results of two variants of Newton iterations for computing a URPO. The searches are initialised with a guess generated using the SRDMD method, that has guessed period $T_g=62.1673$ . The standard single-shot Newton iterations fail to converge when initialised with snapshots corresponding to $t=0$ and $t=T_g/2$ along the reconstructed periodic trajectory (2.7) (red and grey circle markers). The two-shot Newton search, initialised with these two snapshots, converges successfully in $28$ iterations (blue square markers). The vertical axis shows the $L_2$ -norm of the to-be-zeroed vector $G(x)$ in the respective root-finding problem, where $x$ denotes the vector of unknowns.

6. Search for unstable TWs in 3-D plane Poiseuille flow

The simplest category of invariant solutions to the Navier–Stokes equations (3.1) and (3.2) that can provide new insights into the flow are TWs. In this section, we demonstrate using two examples how the sparsity-promoting SRDMD can be used to construct initial conditions for the Newton search of TWs.

When a chaotic trajectory visits the neighbourhood of a steady state, the evolution is nearly linear; hence SRDMD is able to reconstruct trajectory segments accurately in the vicinity of a TW. In addition, when the trajectory is repelled from (or attracted to) the TW’s neighbourhood along its unstable (stable) manifold, the DMD spectrum is expected to be dominated by a few leading eigenvalues of the linearised dynamics. Therefore, we consider the trajectory segment processed by SRDMD to potentially be in the vicinity of a TW if the low-dimensional linear reconstruction is sufficiently accurate while the selected set of SRDMD eigenvalues includes one neutral eigenvalue with the remaining ones being repelling (attracting). In such cases, we can approximate the nearby TW solution by sending $t\to -\infty$ (or $t\to +\infty$ ) in the reconstruction.

Here, we specifically consider the case of purely real eigenvalues. If complex-valued eigenvalues were considered too, then SRDMD could detect a spiralling trajectory, suggesting oscillatory behaviour while the trajectory is being repelled from or attracted to a steady solution. However, since only the SRDMD spectra with real eigenvalues are considered, the trajectory instead resembles a star-type or a saddle-type solution, without oscillations. In this case, (2.7) for the SRDMD reconstruction simplifies to

(6.1) \begin{equation} \boldsymbol{u}(t_j)\approx \sum _{k=0}^{l-1}a_k\,{\rm e}^{\sigma _kt_j}\,\boldsymbol{\psi }_k, \end{equation}

where only the growth or decay rate $\sigma _k$ determines the temporal behaviour of the SRDMD reconstruction. The amplitudes $a_k$ of the purely real dynamic modes $\boldsymbol{\psi }_k$ are again determined via sparsity-promoting optimisation as introduced in § 2. With $l$ , we denote the number of dynamic modes being selected for the reconstruction. This approach is motivated by the idea that the sparsity-promoting optimisation selects dynamic modes hierarchically in accordance with their contribution to an optimal flow reconstruction, and therefore selects purely real dynamic modes prior to dynamic modes with fundamental frequencies, when the flow can be well captured by a non-oscillating SRDMD approximation.

To provide a proof of concept, we chose two distant SRDMD observation time windows near state-space locations where energy input and dissipation balance out, and thus TW solutions can be expected. The SRDMD observation time windows were set to be short: each window spanned 7.5 advective time units, with snapshot sampling $\Delta t$ and spacing $\delta t$ of 0.5 advective time units, to ensure that the flow could be well captured by the linear approximation. Furthermore, we found two real dynamic modes to be sufficient for reconstructing successful initial conditions. Figure 16(a) shows the resulting SRDMD spectra for both cases, where purely real eigenvalues selected through sparsity-promoting optimisation for the flow reconstruction are coloured in blue and orange. In both cases, the neutral eigenvalue $\lambda _0=1$ ( $\sigma _0 = 0$ ) corresponds to the time-independent part of the reconstruction (6.1), while the other eigenvalue $\lambda _1\lt 1$ ( $\sigma _1 \lt 0$ ) corresponds to a decaying mode as $t\to +\infty$ . Taking into account the phase shift $\Delta \phi$ from the symmetry reduction, the reconstructed initial conditions converged quickly to TW solutions by Newton iteration. Figure 16(b) shows the location of both TWs in a 2-D $I$ – $D$ state-space projection, with both quantities normalised by their respective laminar values $I_0$ and $D_0$ . The colour coding of the markers relates the TWs to the dynamic modes in figure 16(a) that were used for constructing the initial conditions. In table 3, we summarise the numerical details and stability properties of both TWs.

Figure 16. (a) The SRDMD spectra obtained from the analysis of two segments of the chaotic trajectory in figure 6, each located at different times within the DNS. Marked in blue and orange are purely real eigenvalues corresponding to the two cases, automatically selected by the sparsity-promotion optimisation when reconstructing the flow. (b) The location of the converged unstable TWs in a state-space projection of energy input $I$ and dissipation $D$ , normalised by their respective laminar values $I_0$ and $D_0$ . The grey-shaded area represents the p.d.f. of the dynamics, as described in figure 9. The marker colours relate the dynamic modes from the SRDMD in (a), selected to construct an initial condition for the Newton search, to the respective converged solutions indicated in (b).

Table 3. The unstable TWs converged from initial data constructed through SRDMD. The spatial resolutions of the computational domain in the streamwise, wall-normal and spanwise direction are $N_x$ , $N_y$ and $N_z$ , respectively, and $c$ is the phase velocity of the TW. Stability properties of the TWs are summarised with the number of unstable directions $N$ , the sum $\sum _{j=1}^N\text{Re} {(\lambda _j)}$ of real parts of the corresponding eigenvalues, and their maximal values $\max \,\text{Re} {(\lambda _j)}$ .