3.1. Spectral submanifolds of the LB travelling wave

In the symmetry-reduced phase space, the LB travelling wave becomes a hyperbolic fixed point with a single unstable eigenvalue. By the classic unstable manifold theorem, there exists a one-dimensional unstable manifold tangent to the unstable eigenvector and a codimension-one stable manifold. The unstable manifold of the LB travelling wave forms a heteroclinic connection with the laminar state. Kaszás et al. (Reference Kaszás, Cenedese and Haller2022) demonstrated that such heteroclinic orbits can serve as one-dimensional reduced-order models, as also supported by the experiments of Suri et al. (Reference Suri, Tithof, Grigoriev and Schatz2017).

To capture more of the transient dynamics near the heteroclinic orbit and obtain an approximation for the basin boundary of the laminar state, we extract here a two-dimensional invariant manifold that can intersect the edge. In order to construct an attracting manifold, we need to take the slowest two-dimensional SSM of the LB travelling wave. This SSM is tangent to both its single unstable eigenvector and its slowest stable eigenvector and hence contains both the heteroclinic orbit and the slowest submanifold of the edge of chaos.

Such a manifold is often called a pseudounstable manifold (de la Llave & Wayne Reference de la Llave and Wayne1995), whose existence can also be concluded from the recent results of Haller et al. (Reference Haller, Kaszás, Liu and Axås2023) or Buza (Reference Buza2023). Haller et al. (Reference Haller, Kaszás, Liu and Axås2023) also find a large variety of generalized SSMs of (2.5), including ones of limited smoothness (fractional SSMs) and others of mixed stability type (mixed-mode SSMs). These results apply in finite dimensions and give $C^{\infty }$-smooth, mixed-mode SSMs under the same non-resonance conditions as those of Sternberg (Reference Sternberg1958).

As discussed by Haller et al. (Reference Haller, Kaszás, Liu and Axås2023), exact resonances among complex eigenvalues are unlikely in a typical finite-dimensional system, such as the numerical discretization (2.5). We thus expect that a unique, $C^\infty$, mixed-mode SSM of the LB travelling wave exists. We list the corresponding non-resonance conditions in the Appendix A and verify that they hold for a subset of the spectrum. Alternatively, one may invoke the results of Buza (Reference Buza2023), who showed that a $C^1$-smooth pseudounstable manifold exists for the Navier–Stokes equations (2.1).

Based on these recent technical developments, we can apply the data-driven methodology of Cenedese et al. (Reference Cenedese, Axås, Bäuerlein, Avila and Haller2022) to discover this mixed-mode SSM from simulation data. We follow the approach of Kaszás et al. (Reference Kaszás, Cenedese and Haller2022) and introduce the square roots of the variables $I$ and $D$ as $J := \sqrt {|I|}$ and $K := \sqrt {|D|}$. This is necessary, because the velocity field $\boldsymbol {u}$ is a non-differentiable function of $I$ and $D$ at the origin, due to their quadratic dependence on the components of $\boldsymbol {u}$. Thus, we parametrize the SSM with the variables $J$ and $K$. Specifically, we seek a two-dimensional invariant manifold in the phase space of (2.14) of the form

(3.1)\begin{equation} \boldsymbol{x}_R = \sum_{1\leq n + m \leq M_p} \boldsymbol{c}_{nm}K^mJ^n, \quad \boldsymbol{c}_{nm}\in \mathbb{R}^d. \end{equation}

The polynomial-type dependence up to order $M_p$ on the reduced coordinates in (3.1) is justified because polynomials are universal approximators (Rudin Reference Rudin1976). It is also motivated by the success of local Taylor approximations used in the original methods (Haller & Ponsioen Reference Haller and Ponsioen2016) for SSM-based model-order reduction.

To determine the coefficient vectors $\boldsymbol {c}_{nm}$ in (3.1) we initialize training trajectories that lie on the mixed-mode SSM of the LB travelling wave. The parametrization of the manifold in physical coordinates can then be recovered from (3.1) using the discretization (2.4) performed by Openpipeflow (Willis Reference Willis2017). This results in

(3.2)\begin{equation} \boldsymbol{u}(r,\varphi,z) = \boldsymbol{W}(J,K) = \boldsymbol{u}_{HP}(r) + \sum_{1\leq n+m\leq M_p}\boldsymbol{w}_{nm}(r, \varphi,z)K^mJ^{n}, \end{equation}

where the coefficient-functions $\boldsymbol {w}_{nm}(r,\varphi,z)$ are determined by $\boldsymbol {c}_{nm}$.

To generate initial conditions, we start with the LB travelling wave discussed by Kerswell & Tutty (Reference Kerswell and Tutty2007) at $Re = 2,400$. The travelling wave is found using a Newton–Krylov scheme (Viswanath Reference Viswanath2007). This method also returns the leading eigenvalues and eigenvectors of the travelling wave (computed in a comoving frame) (see Willis Reference Willis2019).

For training, we take a total of six trajectories, four of which were randomly initialized at a distance of $10^{-4}$ from the LB travelling wave. The distance is measured as the Euclidean distance in the phase space of (2.14). In this metric, the distance between the laminar state and the LB travelling wave is 0.72. We reserve an additional randomly initialized trajectory for validation. After the initial transients, the trajectories approach the unstable manifold of the LB travelling wave. To also capture the dynamics in the slowest stable direction, we construct two more training trajectories that are approximately constrained to the edge. We use the edge-tracking algorithm of Itano & Toh (Reference Itano and Toh2001) started from initial guesses lying along the subspace spanned by the slowest stable eigenvector. The energy input and dissipation variables are computed using Openpipeflow from the spectral representation of the flow fields.

We show the set of training trajectories in figure 2. Figure 2(a) shows the original trajectories as a projection onto the three most dominant spatial modes, which are displayed in figure 2(d). In all further computations, we work with the symmetry-reduced phase space and preprocess the trajectories using the method of slices. The symmetry-reduced trajectories can be seen in figure 2(b). We also show the reduced coordinates, $J=\sqrt {|I|}$ and $K=\sqrt {|D|}$, which coincide for the full-space trajectories and the symmetry-reduced trajectories.

Note that the method of slices generally can only be applied over a bounded domain. When trajectories cross the border of the chart associated with the template (Froehlich & Cvitanović Reference Froehlich and Cvitanović2012), finite jumps are observed in the symmetry-reduced phase space. Figure 2 indicates no such singularities in our training trajectories, which prompts us to represent them by a single chart. Willis et al. (Reference Willis, Cvitanović and Avila2013) found that a global atlas defined using multiple template states was necessary to construct the symmetry-reduced representation of the turbulent state. Since our training trajectories avoid the turbulent state and remain in the neighbourhood of the heteroclinic orbit, using the LB travelling wave as the only template state was sufficient.

Given the training trajectories, we identify the coefficient vectors $\boldsymbol {c}_{nm}$ via a sparsity promoting ridge regression (Hastie, Tibshirani & Friedman Reference Hastie, Tibshirani and Friedman2009; Brunton, Proctor & Kutz Reference Brunton, Proctor and Kutz2016), minimizing the squared norm of the deviation from these trajectories. Once the SSM geometry is identified, we seek the SSM-reduced dynamics in the form

(3.3)\begin{equation} \frac{{\rm d}}{{\rm d}t}\begin{pmatrix}J \\ K \end{pmatrix} = \boldsymbol{f}(J,K)= \sum_{1\leq n+m\leq M_r} \begin{pmatrix} R^{(J)}_{nm} K^m J^n \\ R^{(K)}_{nm} K^m J^n\end{pmatrix}. \end{equation}

Here the coefficients $R^{(J)}_{nm}$ and $R^{(K)}_{nm}$ are also determined by ridge regression onto $\dot {J}$ and $\dot {K}$ obtained by finite differencing along the training trajectories. We also introduce a constraint to this optimization, forcing the LB travelling wave and the laminar state to be fixed points of (3.3), as in Kaszás et al. (Reference Kaszás, Cenedese and Haller2022).

The parameters of the regression are the polynomial orders $M_p$ and $M_r$, as well as the weight of the ridge-type penalty term. For simplicity, we take $M_p=M_r$. This choice is motivated by the usual Taylor expansion representation of SSMs, where the parametrization and the reduced dynamics are computed up to the same order. The value of the polynomial orders and the ridge-type penalty term are determined by cross-validation on a trajectory initialized in the same way as the training set.

We minimize the overall mean-prediction error on the validation trajectory, defined as

(3.4)\begin{equation} \text{Error} = \sum_{i=0}^{n_t} \left \Vert \boldsymbol{u}_{true}(t_i) - \boldsymbol{W} \circ F^{t_i}(J_{true}(0), K_{true}(0))\right \Vert, \end{equation}

where $F^t$ is the time $t$ flow-map of the reduced dynamics (3.3) and $n_t$ is the number of time snapshots available along the validation trajectory. To avoid overfitting and promote simpler models, as long as the value of the error (3.4) is similar, we prioritize lower polynomial orders and larger penalty terms. The results are reported with $M_r=M_p=5$, and the penalty term is chosen as $10^{-5}$. We refer to the supplementary JFM Notebooks for further details.

We also note that the approximation of the SSM benefits from a richer training set. The accuracy of the reduced-order model increases if more training trajectories are used for the regression or if the temporal resolution of the trajectories is increased, as explained by Cenedese et al. (Reference Cenedese, Axås, Bäuerlein, Avila and Haller2022).

The LB travelling wave is a fixed point of the reduced-order model by construction, that is $\boldsymbol {f}(J_{LB}, K_{LB}) = \boldsymbol {0}$. The Jacobian, $D\boldsymbol {f}(J_{LB}, K_{LB})$, has eigenvalues ${\lambda ^{(+)}_{red.} = 0.0199}$ and ${\lambda ^{(-)}_{red.} = -0.0107}$, matching the eigenvalues of the LB in the full-order model obtained by Krylov iteration in Openpipeflow, whose values are ${\lambda ^{(+)} = 0.0198}$ and ${\lambda ^{(-)} = -0.0010}$.

3.2. Predictions of the reduced-order model

To assess the predictive power of the two-dimensional, SSM-reduced model obtained above, we also generate an ensemble of test trajectories distributed over a sphere of radius $10^{-4}$ around the LB travelling wave in the phase space. Some of these trajectories transition to turbulence while others laminarize, but we can make predictions based on the reduced-order model for all of them. In figure 3 we compare the time evolution of the $J$ coordinate of the test trajectories with their predicted counterparts. While the reduced-order model accurately differentiates between turbulent and laminar trajectories, the turbulent trajectories are only reliably modelled for 200 time units, since the model was not trained on trajectories in the turbulent state. Figure 3(b) shows the overall relative error of the predictions (3.4) for the ensemble of test trajectories.

Figure 3.

Predictions of the reduced-order model on an ensemble of 20 test trajectories initialized near the LB travelling wave. In panel (a), the time evolution of their $J$ coordinate is compared with the predictions made by the reduced-order model. For turbulent trajectories, only the first 200 time units of the trajectories are shown. In panel (b), the corresponding relative errors are plotted, defined as the time-dependent quantity $\Vert \boldsymbol {u}_{true}(t_i) - \boldsymbol {W} \circ F^{t_i}(J_{true}(0), K_{true}(0)) \Vert / \max \Vert \boldsymbol {u}_{true} \Vert$. Ensemble averages of the relative error are computed. The errors of turbulent and laminar trajectories are averaged separately.

In figure 4 we show the SSM-reduced model vector field $\boldsymbol {f}$ on the $J$–$K$ plane. We also show the approximation of the stable manifold of the LB travelling wave, obtained by integrating initial conditions of the form

(3.5)\begin{equation} \begin{pmatrix} J_{LB} \\ K_{LB} \end{pmatrix} \pm 10^{-6} \boldsymbol{v} \end{equation}

backward in time, where $\boldsymbol {v}\in \mathbb {R}^2$ is the eigenvector of $D\boldsymbol {f}(J_{LB}, K_{LB})$ associated with the negative eigenvalue. In figure 4, the predicted stable manifold matches the edge trajectory well, even outside the domain of the training data for $J>J_{LB}$.

Figure 4.

Two-dimensional reduced-order model obtained by restricting the dynamics to the two-dimensional mixed-mode SSM of the LB travelling wave. This manifold is parametrized by the variables $J$ and $K$ that define the reduced dynamics. Training trajectories are shown in green. The directory, including the data and the Jupyter notebook that generated this figure can be accessed at https://multimediafluids.uam.elogim.com/S0022112023009564/JFM-Notebooks/files/Figure4/fig4.ipynb.

Having computed the parametrization $\boldsymbol {W}$, we can also display the global shape of the mixed-mode SSM. In figure 5, we highlight the domain in the $J$–$K$ plane, over which we visualize the SSM by plotting its relative kinetic energy,

(3.6)\begin{equation} \Delta E = \tfrac{1}{2}\left \Vert \boldsymbol{u}_{HP} - \boldsymbol{W}(J,K) \right \Vert^2, \end{equation}

in figure 5(b).

Figure 5.

Visualization of the reduced-order model and the geometry of the mixed-mode SSM carrying the model flow. Panel (b) shows the mixed-mode SSM (green surface) in the $(J,K,\Delta E)$ plane obtained as the image of the domain shown in (a) under the parametrization $\boldsymbol {W}$. The relative energy $\Delta E$ is defined as (3.6). The fixed points, the base state and the LB travelling wave are indicated with coloured dots. The dashed line denotes the predicted stable manifold of the LB travelling wave in the SSM-reduced model. The directory, including the data and the Jupyter notebook that generated this figure can be accessed at https://multimediafluids.uam.elogim.com/S0022112023009564/JFM-Notebooks/files/Figure5/fig5.ipynb.

3.3. Capturing the edge of chaos in the reduced-order model

The mixed-mode SSM shown in figure 5(b) is two-dimensional, and the edge manifold (i.e. the stable manifold of the LB travelling wave) has codimension one. Adding up their dimensions results in $d+1$, where $d$ is the dimension of the phase space of (2.14). Therefore, these two manifolds generically intersect transversely along a one-dimensional curve. This implies that the intersection is robust under small perturbations to system (2.5).

The generically expected transverse intersection of the SSM with the edge manifold is displayed in figure 6, where the mixed-mode SSM is the same as in figure 5 and a schematic representation of the edge manifold is added. In the three-dimensional $(J, K, \Delta E)$ space, both the edge manifold and the mixed-mode SSM appear as two-dimensional surfaces, even though in the full phase space the edge manifold is a much higher-dimensional object.

Figure 6.

Mixed-mode SSM and the edge manifold shown in the ($J,K,\Delta E$) space. The edge manifold is sketched as a two-dimensional surface in the three-dimensional space that intersects the edge manifold along the stable manifold of the LB travelling wave in the reduced-order model (dashed black line).

From the spectrum of the LB travelling wave in figure 1, we conclude that its slowest eigenvalue with negative real part is real. Therefore, the stable manifold of this travelling wave has a one-dimensional slow SSM tangent to the eigenspace of that real eigenvalue. This SSM coincides with the stable manifold of the LB travelling wave in the reduced-order model, and hence must be the line of intersection shown in dashed lines in figure 6.

Since the intersection of the edge manifold and the mixed-mode SSM is a curve, it will not capture the dynamics within the whole edge manifold. Nevertheless, this intersection identifies a clear footprint of the edge manifold in the SSM-reduced model. To demonstrate this, we construct 10 trajectories constrained to the edge manifold using the bisection algorithm of Itano & Toh (Reference Itano and Toh2001). These trajectories converge to the LB travelling wave, which is a relative attractor within the edge. In the reduced phase space $(J,K)$, these trajectories clearly approach the LB travelling wave along the predicted stable manifold, as shown in figure 7(a).

Figure 7.

(a) Edge trajectories constructed using the bisection algorithm in the full-order model approaching the LB travelling wave displayed in the reduced-phase space. A single representative member of this ensemble of trajectories is shown with a red dashed line to guide the eye. The stable manifold of the LB travelling wave predicted from the reduced-order model is shown in black. In panel (b), trajectories of the full-order model are initialized by placing initial conditions (ICs) on both sides of the predicted edge in the reduced-order model.

The distinguishing property of the edge manifold is that it divides the phase space locally around the LB travelling wave. We show that its footprint in the reduced-order model also has this dividing role by initializing a set of trajectories in the reduced phase space at a fixed distance $\delta =0.01$ from the stable manifold. Initial conditions are then prepared using the parametrization $\boldsymbol {W}$ and are supplied to the full-order model Openpipeflow. Figure 7(b) shows that the initial conditions prepared in this way are clearly separated in both the SSM-reduced model and the full-order model.

In figure 8(a), we show a grid of initial conditions in a rectangle around the LB travelling wave in the reduced-order model. The initial conditions for the full-order model are prepared using the identified parametrization $\boldsymbol {W}$ and are integrated forward in time. Different colours indicate the different long-time behaviours of the grid of initial conditions. The stable manifold of the reduced-order model clearly separates the trajectories of the full-order model, as expected.

Figure 8.

(a) Phase portrait of the SSM-reduced model with a grid of test trajectories. Initial conditions are placed in a grid around the LB travelling wave and are integrated forward in time using the full-order model. The initial conditions are coloured according to their long-time behaviour, as seen in the insets. The predicted edge manifold footprint is shown in black. (b) Time spent by these trajectories near the edge as a function of their initial distance from the edge measured in the SSM-reduced phase space. This time is defined as the time needed to reach a distance of $D_{max} = 0.1$ from the LB travelling wave. A least-squares fit of the form $T=C - ({1}/{\lambda })\log \xi _0$ is also shown with $C=-234 \pm 14$ and $\lambda = 0.0175 \pm 0.0005$.

Another indication that this separation is indeed caused by the saddle-type edge state is the slowdown of trajectories, which is already present in a linear system. Consider a linear system with a saddle-type fixed point that has a one-dimensional unstable manifold. In diagonal form, the dynamics along the unstable subspace are

(3.7)\begin{equation} \dot{\xi} = \lambda \xi,\quad \text{for } \lambda > 0, \end{equation}

while for the remaining coordinates we have

(3.8)\begin{equation} \dot{\eta}_i=-\kappa_i \eta_i,\quad \text{for } \kappa_i>0,\ i=1,\ldots, d-1. \end{equation}

For sufficiently long times, the distance of a trajectory from the saddle is well approximated by its $\xi$ coordinate component, i.e. the distance $\delta (t)$ can be written as $\delta (t)=\xi _0{\rm e}^{\lambda t}$. Therefore, the time $T$ it takes for a trajectory to leave the $\delta _{max}$-neighbourhood of the saddle satisfies $\delta _{max}=\xi _0 {\rm e}^{\lambda T}$, or, equivalently

(3.9)\begin{equation} T = \frac{1}{\lambda}\log \delta_{max} - \frac{1}{\lambda}\log \xi_0, \end{equation}

where $\xi _0$ measures the initial distance from the stable subspace.

To demonstrate the slowdown, starting from the initial conditions shown in figure 8(a), we record the time needed for the full-order model trajectory to develop a distance of $\delta _{max} =0.1$ from the LB travelling wave. Figure 8(b) shows the dependence of this time on the initial distance from the edge manifold in the reduced model. This relationship can be well described by a logarithmic function of the form (3.9), which would hold for a linear system. Therefore, we expect that $T=C -({1}/{\lambda })\log \xi _0$ is satisfied approximately for some constants $C$ and $\lambda$. We obtain $\lambda = 0.0175$ from a least squares fit to the data, which reasonably matches the true unstable eigenvalue of the LB travelling wave, $\lambda ^{(+)}$.

Although the one-dimensional footprint of the edge manifold cannot act as a global barrier in the phase space, it is still an influential curve for trajectories near the LB travelling wave. Since the mixed-mode SSM is constructed to be tangent to the slowest dynamics, it is locally attracting. Therefore, one may project nearby trajectories onto the mixed-mode SSM to decide their long-term dynamics: if the projection lies above the footprint of the edge manifold within the SSM, then the trajectory is expected to become turbulent. Otherwise, the trajectory is expected to quickly laminarize. As one moves away from the vicinity of the LB, this correspondence will gradually break down due to the overall complicated shape of the edge (Schneider et al. Reference Schneider, Eckhardt and Yorke2007; Mellibovsky et al. Reference Mellibovsky, Meseguer, Schneider and Eckhardt2009).

3.4. Parameter dependent SSM-reduced models

The spectrum of the LB travelling wave changes smoothly as the Reynolds number is varied, therefore the slowest mixed-mode SSM remains two-dimensional for a wide range of Reynolds numbers. As in Kaszás et al. (Reference Kaszás, Cenedese and Haller2022), we may then seek a parametrization of the Reynolds number-dependent family of mixed-mode SSMs. Similarly, the reduced dynamics on the family of manifolds then also depend on the Reynolds number. Formally, this means that we seek $\boldsymbol {W}$ and $\boldsymbol {f}$ in a Reynolds number-dependent form

(3.10)\begin{equation} \boldsymbol{u}(r,\varphi,z) = \boldsymbol{W}(J,K, {Re}) = \boldsymbol{u}_{HP}(r) + \sum_{l=0}^{N_p}\sum_{1\leq n+m\leq M_p}\boldsymbol{w}_{nml}(r, \varphi,z)K^mJ^{n} {Re}^l \end{equation}

and

(3.11) \begin{equation} \frac{{\rm d}}{{\rm d}t}\begin{pmatrix}J \\ K \end{pmatrix} = \boldsymbol{f}(J,K, {Re})= \sum_{l=0}^{N_p}\sum_{1\leq n+m\leq M_r} \begin{pmatrix} R^{(J)}_{nml} K^m J^n \, {Re}^l \\ R_{nml}^{(K)} K^m J^n \, {Re}^l\end{pmatrix}. \end{equation}

The general dependence of these expressions on the Reynolds number is restricted by the requirement that the laminar state must be a fixed point at $J=K=0$ for all values of the Reynolds number, therefore $\boldsymbol {w}_{00l} =0$ and $R_{00l}=0$ must hold for all $l$.

To construct a parameter-dependent SSM-reduced model, we select the range ${{Re}\in [2400, 2520]}$, where the slowest eigenvalues of the LB travelling wave do not change considerably and hence the mixed-mode SSMs are only expected to show minor variation. In this case, we can take a linear approximation for the Reynolds-number dependence in (3.10) and (3.11), i.e. $N_p=1$. We use three sets of training trajectories initialized at ${Re} = 2400; 2550; 2520$. These are prepared in the same way as described in the case of ${Re}=2400$. Figure 9 shows this parametrized family of reduced dynamics in the (${Re}, J, K$) space. We find that the phase portraits of the reduced models are qualitatively the same for all Reynolds numbers in this range.

Figure 9.

Family of SSM-reduced models parametrized by the Reynolds number. Phase portraits at the sections $Re = 2400$ and $Re = 2520$ are shown, along with the LB travelling wave and its stable and unstable manifolds computed from the SSM-reduced models.