2 Numerical set-up and notation

For numerical simulations, we use Openpipeflow (Willis Reference Willis2017), which integrates the Navier–Stokes equations for fluctuations $\boldsymbol{u}$ around the base (Hagen–Poiseuille) solution. The axial pressure gradient is adjusted throughout the simulation in order to ensure a constant flux equal to that of the base flow at a given $\mathit{Re}=UD/\unicode[STIX]{x1D708}$ , where $U$ is the mean axial velocity, $D$ is the pipe diameter and $\unicode[STIX]{x1D708}$ is the kinematic viscosity of the fluid. The flow field satisfies the incompressibility condition $\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}=0$ throughout the volume; there are periodic boundary conditions $\boldsymbol{u}(z,\unicode[STIX]{x1D703},r)=\boldsymbol{u}(z+L,\unicode[STIX]{x1D703},r)$ and $\boldsymbol{u}(z,\unicode[STIX]{x1D703},r)=\boldsymbol{u}(z,\unicode[STIX]{x1D703}+2\unicode[STIX]{x03C0},r)$ in axial and azimuthal directions, and the no-slip boundary condition $\boldsymbol{u}(z,\unicode[STIX]{x1D703},r=D/2)=0$ on the pipe wall. For spatial discretization, flow fields are expanded in Fourier series in axial and azimuthal directions and finite differences are used to evaluate derivatives in the radial direction. In the numerical work presented here, the computational domain is $L=25D$ long and Fourier series in axial and azimuthal directions are respectively truncated at $K=192$ and $M=16$ , and $N=64$ finite difference points are used in the radial direction. Nonlinear terms are evaluated following the $3/2$ -rule for dealiasing resulting in $64\times 576\times 48$ spatial grid points. We adopted our resolution from Chantry, Willis & Kerswell (Reference Chantry, Willis and Kerswell2014), who studied bifurcations of the localized solutions considered here from travelling waves as the pipe length is increased. We checked our results for robustness at higher spatial resolution $(N,K,M)=(128,320,32)$ . Re is set to 1700, at which puffs have long lifetimes (up to $1000\,D/U$ ) in the symmetry subspace considered here.

Pipe flow is equivariant under streamwise translations $g_{z}(l)\boldsymbol{u}(z,r,\unicode[STIX]{x1D703})=\boldsymbol{u}(z-l,r,\unicode[STIX]{x1D703}),$ azimuthal rotations $g_{\unicode[STIX]{x1D703}}(\unicode[STIX]{x1D719})\boldsymbol{u}(z,r,\unicode[STIX]{x1D703})=\boldsymbol{u}(z,r,\unicode[STIX]{x1D703}-\unicode[STIX]{x1D719})$ and the azimuthal reflection $\unicode[STIX]{x1D70E}[u,v,w](z,r,\unicode[STIX]{x1D703})=[u,v,-w](z,r,-\unicode[STIX]{x1D703}),$ where $u,v,w$ are velocity field components in the axial, radial and azimuthal directions respectively. Following Avila et al. (Reference Avila, Mellibovsky, Roland and Hof2013), we restrict our study to the velocity fields that are symmetric under rotation by $\unicode[STIX]{x03C0}$ , $g_{\unicode[STIX]{x1D703}}(\unicode[STIX]{x03C0})\boldsymbol{u}=\boldsymbol{u}$ , and the reflection $\unicode[STIX]{x1D70E}\boldsymbol{u}=\boldsymbol{u}$ . Imposing reflection invariance breaks the continuous rotation symmetry of the system, allowing only for half-domain rotations. Since we are considering the subspace invariant under rotation by $\unicode[STIX]{x03C0}$ , half-domain rotations are given by $g_{\unicode[STIX]{x1D703}}(\unicode[STIX]{x03C0}/2)$ , hence the symmetry group of the system becomes $G=\{g_{z}(l),g_{\unicode[STIX]{x1D703}}(\unicode[STIX]{x03C0}/2)\}$ .

For clarity, we are going to use state-space notation. Let $a(t)$ be a vector containing all $(3\times 64\times 384\times 32)$ numerical degrees of freedom of the flow fields $[u,v,w]$ . Evolution under Navier–Stokes equations implies a finite-time flow mapping $f^{\unicode[STIX]{x1D70F}}(a(0))=a(\unicode[STIX]{x1D70F})$ that takes a solution at time $0$ to a new one at time $\unicode[STIX]{x1D70F}$ . In this formalism, equivariance under $g\in G$ means that the flow and the symmetry operation commute, i.e. $gf^{\unicode[STIX]{x1D70F}}(a(0))=f^{\unicode[STIX]{x1D70F}}(ga(0))$ . Throughout this article, our notation will not distinguish between an abstract group element and its particular representation. Thus, whenever a group action is present, its appropriate representation on the corresponding velocity fields is implied. Finally, we need to introduce an inner product to use in our calculations. Our choice is the standard ‘energy norm’: let $a$ and $a^{\prime }$ be state-space vectors corresponding to velocity fields $\boldsymbol{u}$ and $\boldsymbol{u}^{\prime }$ . We define the $L_{2}$ inner product as $\langle a,\,a^{\prime }\rangle =(1/2)\int \boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{u}^{\prime }\,\text{d}V$ ; hence $\Vert a\Vert ^{2}=\langle a,\,a\rangle$ is the kinetic energy of $a$ .

A relative periodic orbit is a recurrence (after the orbit’s period $T_{p}^{\prime }$ ) up to a symmetry operation

(2.1) $$\begin{eqnarray}f^{T_{p}^{\prime }}(a_{p})=g_{p}^{\prime }a_{p}.\end{eqnarray}$$

For the pair of relative periodic orbits found by Avila et al. (Reference Avila, Mellibovsky, Roland and Hof2013), $g_{p}^{\prime }=g_{\unicode[STIX]{x1D703}}(\unicode[STIX]{x03C0}/2)g_{z}(l_{p}^{\prime })$ . For simplicity, we are going to treat this orbit as if its period were twice its fundamental period $T_{p}=2T_{p}^{\prime }$ , and its only symmetry is a streamwise shift since $g_{p}=(g_{p}^{\prime })^{2}=g_{z}(l_{p}=2l_{p}^{\prime })$ .

A relative periodic orbit with a one-parameter compact continuous symmetry defines a $2$ -torus

(2.2) $$\begin{eqnarray}\{g_{z}(l)f^{\unicode[STIX]{x1D70F}}(a_{p})\,|\,\unicode[STIX]{x1D70F}\in [0,T_{p}),l\in [0,L)\}\end{eqnarray}$$

in the state space, which is parametrized by shifts $l$ and time $\unicode[STIX]{x1D70F}$ . The stability of a relative periodic orbit is determined by the eigenvalues of the Jacobian matrix

(2.3) $$\begin{eqnarray}\unicode[STIX]{x1D645}_{p}=g_{p}^{-1}\,\text{d}f^{T_{p}}(a_{p})/\text{d}a_{p},\end{eqnarray}$$

which are known as Floquet multipliers (Chossat & Lauterbach Reference Chossat and Lauterbach2000; Cvitanović et al. Reference Cvitanović, Artuso, Mainieri, Tanner and Vattay2015). Since this Jacobian matrix (2.3) is very large, in practice, we compute the leading Floquet multipliers $\unicode[STIX]{x1D6EC}_{i}$ and the corresponding Floquet vectors $V_{i}$ via Arnoldi iteration (Trefethen & Bau Reference Trefethen and Bau1997). The lower-branch relative periodic orbit of Avila et al. (Reference Avila, Mellibovsky, Roland and Hof2013) has only one positive real Floquet multiplier $\unicode[STIX]{x1D6EC}_{1}$ greater than one, two marginal $\unicode[STIX]{x1D6EC}_{2,3}=1$ multipliers corresponding to the disturbances as axial and temporal shifts; the rest $|\unicode[STIX]{x1D6EC}_{i>3}|$ of the multipliers have absolute values smaller than one. Counting axial and temporal shift directions, the unstable manifold of this relative periodic orbit is three-dimensional. However, these marginal directions are of no dynamical importance and our next step is to cancel them.

3 Reductions

For a $(1+1)$ -dimensional partial differential equation under a periodic boundary condition, Budanur et al. (Reference Budanur, Cvitanović, Davidchack and Siminos2015) show that the translational symmetry can be reduced by fixing the phase of the first Fourier mode to $0$ , and this transformation can be interpreted as a ‘slice’, that is a codimension-1 hyperplane where all symmetry-equivalent state-space points are represented by a single point. Extension of this method to the canonical shear flows requires a particular choice for the slice template’s dependence on the non-homogeneous directions. Willis, Short & Cvitanović (Reference Willis, Short and Cvitanović2016) adapted this idea to the study of a short $(L\approx 1.85)$ pipe with imposed shift-and-reflect symmetry: they experimentally selected a typical turbulent state for this purpose, retained its dependence on the radial and azimuthal coordinates, while imposing a $\cos (2\unicode[STIX]{x03C0}z/L)$ dependence in the axial direction. Here, we take a much simpler approach and define a first Fourier mode slice template $\hat{a}^{\prime }$ as the state vector corresponding to the three-dimensional field

(3.1) $$\begin{eqnarray}[\hat{u} ^{\prime },\hat{v}^{\prime },{\hat{w}}^{\prime }](z,\unicode[STIX]{x1D703},r)=\text{J}_{0}(\unicode[STIX]{x1D6FC}r)\cos (2\unicode[STIX]{x03C0}z/L),\end{eqnarray}$$

where $\text{J}_{0}$ is the zeroth Bessel function of the first kind and $\unicode[STIX]{x1D6FC}$ is chosen such that $\text{J}_{0}(\unicode[STIX]{x1D6FC}D/2)=0$ . For a given trajectory $a(\unicode[STIX]{x1D70F})$ , we can now define a symmetry-reduced trajectory as

(3.2) $$\begin{eqnarray}\hat{a}(\unicode[STIX]{x1D70F})=g_{z}(L\unicode[STIX]{x1D719}/2\unicode[STIX]{x03C0})a(\unicode[STIX]{x1D70F}),\end{eqnarray}$$

where

(3.3) $$\begin{eqnarray}\unicode[STIX]{x1D719}(\unicode[STIX]{x1D70F})=\arg (\langle a(\unicode[STIX]{x1D70F}),\,\hat{a}^{\prime }\rangle +i\langle a(\unicode[STIX]{x1D70F}),\,g_{z}(-L/4)\hat{a}^{\prime }\rangle ).\end{eqnarray}$$

The pair of state-space vectors $\hat{a}^{\prime },g_{z}(-L/4)\hat{a}^{\prime }$ are orthogonal to each other and span a two-dimensional subspace. As depicted in figure 1(a), the slice-fixing phase $\unicode[STIX]{x1D719}$ in (3.3) is the polar angle when a state is projected onto this hyperplane. Transformation (3.2) fixes this angle to $0$ , defining a unique symmetry-reduced $\hat{a}(\unicode[STIX]{x1D70F})$ for all $a(\unicode[STIX]{x1D70F})$ . We have chosen the $r$ dependence of the slice template in (3.1) to be a Bessel function because of the cylindrical geometry. In practice, many other choices can be equally good for the purpose of symmetry reduction.

Figure 1. Sketches of two different formulations of continuous symmetry reduction by the first Fourier mode slice. (a) Projection of a full state-space trajectory $a(\unicode[STIX]{x1D70F})$ on a two-dimensional hyperplane spanned by $(\hat{a}^{\prime },g_{z}(-L/4)\hat{a}^{\prime })$ . Symmetry-fixing phases $\unicode[STIX]{x1D719}(\unicode[STIX]{x1D70F})$ (3.3) are shown for instances $\unicode[STIX]{x1D70F}_{1}$ and $\unicode[STIX]{x1D70F}_{2}$ . (b) Slice hyperplane defined by (3.4) and the transformation of full state-space trajectory $a(\unicode[STIX]{x1D70F})$ onto the slice by the action of $g_{z}(L\unicode[STIX]{x1D719}(\unicode[STIX]{x1D70F})/2\unicode[STIX]{x03C0})$ .

Our next step is to redefine the transformation (3.2) as a slice, a codimension-1 half-hyperplane

(3.4a,b ) $$\begin{eqnarray}\langle \hat{a}(\unicode[STIX]{x1D70F})-\hat{a}^{\prime },\,t_{}^{\prime }\rangle =0,\quad \langle t(\hat{a}),\,t_{}^{\prime }\rangle >0,\end{eqnarray}$$

where $t(\hat{a})=T_{z}\hat{a}$ is the group tangent, $t_{}^{\prime }=t(\hat{a}^{\prime })$ is the slice tangent, $T_{z}$ is the generator of infinitesimal translations, satisfying $g_{z}(l)=\exp (T_{z}l)$ . In our particular case, $T_{z}\boldsymbol{u}=-\text{d}\boldsymbol{u}/\text{d}z$ , hence $t_{}^{\prime }=(2\unicode[STIX]{x03C0}/L)g_{z}(-L/4)\hat{a}^{\prime }$ . After defining the slice by (3.4), one looks for the phases $\unicode[STIX]{x1D719}$ such that $g_{z}(L\unicode[STIX]{x1D719}/2\unicode[STIX]{x03C0})a(\unicode[STIX]{x1D70F})$ satisfies (3.4) as sketched in figure 1(b).

At first sight, redefining the polar coordinate transformation (3.2) as a slice (3.4) might seem an overcomplication; however, this reformulation provides some important tools. In particular, one can derive (see appendix A) the projection operator

(3.5) $$\begin{eqnarray}\unicode[STIX]{x1D643}(\hat{a})=1-\frac{t(\hat{a})\otimes t_{}^{\prime }}{\langle t(\hat{a}),\,t_{}^{\prime }\rangle }\end{eqnarray}$$

that projects infinitesimal perturbations to $\hat{a}$ from full state space to the slice. With (3.4) and (3.5), we can now reduce the torus (2.2) to a closed curve and define its stability. Let $a_{p}$ be a point on a relative periodic orbit and $V_{i}$ be the Floquet vectors computed at this point. If $\hat{a}_{p}=g(L\unicode[STIX]{x1D719}_{p}/2\unicode[STIX]{x03C0})a_{p}$ is the symmetry-reduced state-space point corresponding to $a_{p}$ , then the symmetry-reduced Floquet vectors are $\hat{V}_{i}=\unicode[STIX]{x1D643}(\hat{a}_{p})g(L\unicode[STIX]{x1D719}_{p}/2\unicode[STIX]{x03C0})V_{i}$ . Note that $\unicode[STIX]{x1D643}(\hat{a})t(\hat{a})=0$ , thus the continuous symmetry direction is eliminated in the slice.

It can be shown (Guckenheimer & Holmes Reference Guckenheimer and Holmes1983; Cvitanović et al. Reference Cvitanović, Artuso, Mainieri, Tanner and Vattay2015) that in the vicinity of a periodic orbit of a dynamical system, one can define a Poincaré map, which contains a fixed point, whose stability multipliers are equal to the Floquet multipliers of the periodic orbit, except for the marginal multiplier corresponding to the time direction that is eliminated by the Poincaré section. Moreover, the stability eigenvectors in the Poincaré section can be obtained from Floquet vectors by a projection. We define such a Poincaré section as

(3.6a,b ) $$\begin{eqnarray}\langle \hat{a}_{{\mathcal{P}}}-\hat{a}_{p},\,\hat{v}(\hat{a}_{p})\rangle =0,\quad \langle \hat{v}(\hat{a}_{{\mathcal{P}}}),\,\hat{v}(\hat{a}_{p})\rangle >0,\end{eqnarray}$$

where $\hat{v}(\hat{a})=\unicode[STIX]{x1D643}(\hat{a})v(\hat{a})$ and $v(a)=\lim _{\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D70F}\rightarrow 0}(\,f^{\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D70F}}(a)-a)/\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D70F}$ is the state-space velocity. Continuous-time symmetry-reduced flow $\hat{a}(\unicode[STIX]{x1D70F})=\hat{f}^{\unicode[STIX]{x1D70F}}(\hat{a}(0))$ induces a discrete-time dynamical system $\hat{a}_{{\mathcal{P}}}[n]=\hat{f}_{{\mathcal{P}}}^{n}(\hat{a}_{{\mathcal{P}}}[0])$ , where $n$ counts the number of intersections of trajectories with the Poincaré section as the discrete-time variable. Finally, we define the operator that reduces the infinitesimal perturbations to $\hat{a}_{{\mathcal{P}}}$ from the symmetry-reduced state space to the Poincaré section as

(3.7) $$\begin{eqnarray}\unicode[STIX]{x1D64B}(\hat{a}_{{\mathcal{P}}})=1-\frac{\hat{v}(\hat{a}_{{\mathcal{P}}})\otimes \hat{v}(\hat{a}_{p})}{\langle \hat{v}(\hat{a}_{{\mathcal{P}}}),\,\hat{v}(\hat{a}_{p})\rangle },\end{eqnarray}$$

which allows us to project the symmetry-reduced Floquet vectors $\hat{V}_{i}$ onto the Poincaré section as $\hat{V}_{i,{\mathcal{P}}}=\unicode[STIX]{x1D64B}(\hat{a}_{p})\hat{V}_{i}$ . Similar to (3.5), the projection operator (3.7) follows from (A 6).

Figure 2. (a) Unstable manifold of LB (grey dots) and a small perturbation (cyan circles) to the laminar flow that develops into a puff on the Poincaré section (3.6). The direction of discrete time from $n=0$ to $n=9$ is indicated by arrows for the perturbation. Intersections of relative periodic orbits LB and UB with the Poincaré section are marked yellow and black respectively. Magenta markers correspond to the orbit with the initial condition $\hat{a}_{{\mathcal{P}}}(\unicode[STIX]{x1D6FF}=0)$ . Few points on the unstable manifold to the left of LB are shown because this direction consists of laminarizing orbits. (b) Three-dimensional projection of the trajectories (grey) that populate the unstable manifold of LB along with the LB (yellow) and UB (red) within the slice. The laminar equilibrium is annotated with $L$ .

4 The unstable manifold

Budanur & Cvitanović (Reference Budanur and Cvitanović2017) demonstrate the approximation of one- and two-dimensional unstable manifolds of relative periodic orbits on Poincaré sections for the Kuramoto–Sivashinsky system. The main idea is to select trajectories that approximately cover the linear unstable manifold, hence their forward integration approximately covers the nonlinear unstable manifold; see Budanur & Cvitanović (Reference Budanur and Cvitanović2017) for details. Here, we are interested in computing the unstable manifold of the lower-branch relative periodic orbit (LB) of Avila et al. (Reference Avila, Mellibovsky, Roland and Hof2013) at $\mathit{Re}=1700$ . At this Re, the only unstable Floquet multiplier of LB is $\unicode[STIX]{x1D6EC}_{1}=2.8291523$ . Initial conditions on the Poincaré section that approximately cover the associated unstable manifold are

(4.1) $$\begin{eqnarray}\hat{a}_{{\mathcal{P}}}(\unicode[STIX]{x1D6FF})=\hat{a}_{p,{\mathcal{P}}}\pm \unicode[STIX]{x1D716}\unicode[STIX]{x1D6EC}_{1}^{\unicode[STIX]{x1D6FF}}\hat{V}_{1,{\mathcal{P}}},\end{eqnarray}$$

where $\unicode[STIX]{x1D6FF}\in [0,1)$ and $\unicode[STIX]{x1D716}$ is a small constant. Note that under discrete-time linear dynamics, the first return of the initial condition $\hat{a}_{{\mathcal{P}}}(0)$ corresponds to $\hat{a}_{{\mathcal{P}}}(1)$ ; therefore, forward integrating these initial conditions populates the one-dimensional unstable manifold on the Poincaré section. Openpipeflow normalizes Floquet vectors such that their magnitudes are equal to that of the initial point on the relative periodic orbit, i.e. $\Vert V_{i}\Vert =\Vert a_{p}\Vert$ . Projection onto the Poincaré section and setting $\unicode[STIX]{x1D716}=10^{-4}$ yields the relative amplitude of the initial perturbation as $\unicode[STIX]{x1D716}\Vert \hat{V}_{1,{\mathcal{P}}}\Vert /\Vert \hat{a}_{p,{\mathcal{P}}}\Vert \approx 0.68\times 10^{-4}$ . In order to approximate the unstable manifold, we chose nine equidistant points in $[0,1)$ for $\unicode[STIX]{x1D6FF}$ and forward integrated these initial conditions while recording their intersections with the Poincaré section (3.6).

Figure 2 shows the unstable manifold approximated this way. In panel (a), we show the first $15$ intersections of each trajectory with the Poincaré section, projected onto bases formed by the unstable Floquet vector $\hat{V}_{1}$ and the real part of the least stable Floquet vector $\hat{V}_{4}$ of LB, orthonormalized by the Gram–Schmidt procedure, namely $e_{1}=\langle \hat{a}_{{\mathcal{P}}}-\hat{a}_{p},\hat{V}_{1,\bot }\rangle ,e_{2}=\langle \hat{a}_{{\mathcal{P}}}-\hat{a}_{p},\text{Re}\hat{V}_{4,\bot }\rangle$ , where subscript $\bot$ implies that vectors are orthonormalized. We observed that as the trajectories leave the neighbourhood of LB (the origin of figure 2) they approach the upper branch (UB). In figure 2(b) we visualize same trajectories as three-dimensional projections from the first Fourier mode slice along with LB and UB, where for the third projection direction we used the symmetry-reduced state-space velocity vector evaluated at $\hat{a}_{p}$ , i.e. $e_{3}=\langle \hat{a}(\unicode[STIX]{x1D70F})-\hat{a}_{p},\hat{v}(\hat{a}_{p})/\Vert \hat{v}(\hat{a}_{p})\Vert \rangle$ . Here, relative periodic orbits are closed periodic orbits as expected; and the unstable manifold appears two-dimensional. The manifold’s approach to the UB is also clearly visible on this projection. Additionally, we are also able to visualize laminarizing part of the unstable manifold.

Another trajectory shown in figure 2(a) is an initially small perturbation developing into a puff. We generated this perturbation by taking a typical puff state and scaling its amplitude by $0.5$ . Nine discrete-time steps displayed in Figure 2(a) show its approach to UB before redeveloping into a puff. We observed this behaviour in numerical experiments prepared with similar initial conditions. Of course, this is not extremely surprising: if we started with an initial condition that is only large enough for triggering turbulence, then the trajectory would have closely followed the unstable manifold of LB since the numerical evidence strongly suggests that the laminar–turbulent boundary is set by its stable manifold. Nevertheless, it is worth noting that the signatures of UB could be seen in puff formation without the need to prepare very precise initial conditions through a bisection procedure.

Figure 3. (a) $L_{2}$ distance from UB on the Poincaré section for $20$ returns. (b) Time evolution of kinetic energy in units of the kinetic energy $k_{HP}$ of the laminar flow. Magenta (solid): a trajectory on the unstable manifold of LB, cyan (dashed): puff transition initiated by a small perturbation to the laminar state.

In order to confirm apparent approaches to the UB, we measured the trajectories’ distance from it on the Poincaré section. Figure 3(a) shows the distance of the orbit with initial condition $\hat{a}_{{\mathcal{P}}}(\unicode[STIX]{x1D6FF}=0)$ on the unstable manifold of LB (shown magenta in figure 2) and the perturbation we show in figure 2 from the UB. For both trajectories, we see a clear initial drop before they move away following the unstable manifold of the UB. For further comparison, we show the time evolution of turbulent kinetic energy for both trajectories in figure 3(b). Note for the trajectory on the unstable manifold that before the kinetic energy goes up to puff levels ( ${\sim}1.4k_{HP}$ ) it oscillates around $0.5k_{HP}$ during $\unicode[STIX]{x1D70F}\in (125,175)(D/U)$ . Similarly, for the small perturbation, after an initial increase, the kinetic energy stays close to $0.6k_{HP}$ during $\unicode[STIX]{x1D70F}\in (25,75)(D/U)$ before further increasing to puff levels. Both episodes correspond to the approach of trajectories to UB. Note that the time interval shown in figure 3(a) is not necessarily the same for each orbit, nor are the discrete-time intervals equal to each other. For the perturbation, the interval $n\in [0,20]$ corresponds to $\unicode[STIX]{x1D70F}\in (0,146.25)(D/U)$ , whereas for the trajectory on the unstable manifold, it corresponds to $\unicode[STIX]{x1D70F}\in (0,193.67)(D/U)$ .

Figure 4. Isosurfaces of streamwise vorticity at $\unicode[STIX]{x1D714}_{z}=\pm 2.0\,U/D$ (purple and green) and streamwise velocity at $u=\pm 0.4\,U$ (red and blue) for the trajectory on LB’s unstable manifold at $n=0,9,20$ (a,c,f), for the small perturbation at $n=0,4$ (b,d) and for UB (e). Shown here is the pipe region $z\in [5,20]D,\unicode[STIX]{x1D703}\in [0,\unicode[STIX]{x03C0}/2]$ since the rest of the structures can be obtained from reflection and rotation by $\unicode[STIX]{x03C0}$ symmetries that are present.

For further comparison, we visualized streamwise velocity and vorticity isosurfaces for three snapshots $(n=0,9,20)$ on the orbit with initial condition $\hat{a}_{{\mathcal{P}}}(\unicode[STIX]{x1D6FF}=0)$ (magenta on figure 2) respectively in figure 4(a,c,f). At the closest approach ( $n=9$ , figure 4 c) the resemblance of flow structures of the unstable manifold to those of UB (figure 4 e) is very close. In addition, we show flow structures for two snapshots of the orbit shown with cyan markers in figure 2(a) at times $n=0,4$ respectively in figure 4(b,d). While not as dramatic, at the closest approach ( $n=4$ , figure 4 d) to the UB, we also see structural resemblance. The same isosurfaces for the LB are visually indistinguishable from the initial point ( $n=0$ , figure 4 a) on the unstable manifold, and hence not shown separately in figure 4. Note that puffs are structurally (figure 4 f) much more complicated than LB and UB.

In order to compare sizes of turbulent structures, we plotted the kinetic energy of fluctuations as a function of axial position, i.e. $k(z)=(1/2)\int _{0}^{D/2}\int _{0}^{2\unicode[STIX]{x03C0}}\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{u}r\,\text{d}r\,\text{d}\unicode[STIX]{x1D703}$ , for the trajectory on the unstable manifold of LB. Figure 5(a) shows the time interval $n\in [0,9]$ , during which the orbit leaves the neighbourhood of LB and approaches UB. During the time interval $n\in [10,20]$ shown in figure 5(b), the orbit moves away from UB and becomes a puff.

5 Conclusion and outlook

Figure 5. Fluctuating kinetic energy $k(z)$ contained at an axial location $(z)$ for $n=0,1,\ldots ,9$ (a) and $n=10,11,\ldots ,20$ (b). THe darkest shade of blue corresponds to initial time $n=0$ , and the plotting colour becomes lighter as time increases. Dashed yellow and red curves correspond to the LB and the UB. Values of $k(z)$ are normalized by its maximum at $n=19$ .

In this paper, we presented strong numerical evidence that the unstable manifold of a localized relative periodic orbit of pipe flow approaches the neighbourhood of the corresponding upper-branch solution. It is remarkable that this dynamical structure in the state space persists at $\mathit{Re}=1700$ , which is much higher than $\mathit{Re}\approx 1430$ , at which these orbits are born of a saddle-node bifurcation. Even though the UB is stable (within the symmetry subspace) initially after the bifurcation, it quickly undergoes a series of bifurcations. At $\mathit{Re}=1700$ , the Arnoldi method yields seven unstable Floquet multipliers ( $\unicode[STIX]{x1D6EC}>1$ ) (three pairs of complex conjugates and one real), rendering its unstable manifold nine-dimensional including marginal directions. We repeated our calculation of the unstable manifold at $\mathit{Re}=1900$ and obtained a qualitatively similar picture.

While the qualitative features of the unstable manifold’s approach towards the UB appear to be robust, proving whether or not is there an exact connection between the unstable manifold of LB and the stable manifold of UB is a highly challenging computational task, which we did not address here. Since the sum of the dimensions of LB’s unstable manifold and UB’s stable manifold does not exceed the dimension of the state space, one cannot expect a generic intersection of these manifolds without fine-tuning the parameters Re and $L$ . Our computational method is ultimately prone to the dynamical instabilities of UB, which become effective in its vicinity. Experimenting with $\unicode[STIX]{x1D6FF}$ values in (4.1) yielded only a marginal change in the closest distance between the unstable manifold and the UB. The numerical representation of the manifold can be improved by implementing van Veen, Kawahara & Atsushi (Reference van Veen, Kawahara and Atsushi2011)’s formulation of the boundary value problem after symmetry reduction; however, this would be very heavy on computational resources because the problem studied here has more than two orders of magnitude more computational degrees of freedom than the ones they consider.

Previous studies (Halcrow et al. Reference Halcrow, Gibson, Cvitanović and Viswanath2009; van Veen & Kawahara Reference van Veen and Kawahara2011) reported homoclinic and heteroclinic connections in Couette flow, which can accommodate equilibria and periodic orbits without spatial drifts. This is not the case for pressure-driven flows since all invariant solutions except the laminar equilibrium have non-zero streamwise drifts. In this regard, the current work fills an important technical gap and provides a new set of tools for study of pressure-driven flows. While the methods we used here appeared in different publications (Budanur et al. Reference Budanur, Cvitanović, Davidchack and Siminos2015; Ding et al. Reference Ding, Chaté, Cvitanović, Siminos and Takeuchi2016; Budanur & Cvitanović Reference Budanur and Cvitanović2017), where they are applied to much lower-dimensional systems, this is their first application to the full Navier–Stokes equations.

While Avila et al. (Reference Avila, Mellibovsky, Roland and Hof2013) show that the chaotic motion emerges from the bifurcations of the UB, as the Reynolds number is increased from $1430$ to $1545$ , the role of the UB in transition is not obvious at $\mathit{Re}=1700$ . Our result shows that far from the bifurcation point, the UB takes the role of mediating the transition. Note that the UB, whose kinetic energy swings around $k=0.5k_{HP}$ , is energetically separated from turbulent puffs, which have typical kinetic energies $k\approx 1.2k_{HP}$ . This is also clearly visible from the axial distribution of kinetic energies in figure 5 and flow structures in figure 4, where the UB clearly has a much simpler structure than puffs.

The time evolution of the spatial distribution of the kinetic energy in figure 5 suggests a two-stage transition scenario, where the spatial complexity of a puff forms as trajectories follow the unstable manifold of the UB (figure 5 b). Mellibovsky et al. (Reference Mellibovsky, Meseguer, Schneider and Eckhardt2009) studied the transition to turbulence in a long computational domain and observed a two-stage process similar to ours at $\mathit{Re}=2800$ . Via ‘edge tracking’ (Skufca, Yorke & Eckhardt Reference Skufca, Yorke and Eckhardt2006) they found that a localized disturbance with chaotic dynamics resides in between the laminar and turbulent regions of the state space. When they studied the transition to turbulence from this state, they observed that the initial linear energy amplification is followed by the spread of turbulent structures accompanied by an exponential increase in the kinetic energy. Our results show that, in our symmetry-restricted setting at lower Re, this change from local amplification to spatial expansion takes place when the unstable manifold of the LB reaches the vicinity of the UB.

Ritter, Mellibovsky & Avila (Reference Ritter, Mellibovsky and Avila2016) recently reported a detailed numerical study, where they used turbulent kinetic energy and pressure gradient as indicators to support the hypothesis that spatial complexity in this system arises as different chaotic regions in the state space merge with the neighbourhood of the UB. These observations along with ours motivate a detailed study of the UB’s unstable manifold in order to understand the spatial expansion of chaotic spots. The tools we introduce here can be useful for such a study.