Modified snaking in plane Couette flow with wall-normal suction

A specific family of spanwise-localised invariant solutions of plane Couette flow exhibits homoclinic snaking, a process by which spatially localised invariant solutions of a nonlinear partial differential equation smoothly grow additional structure at their fronts while undergoing a sequence of saddle-node bifurcations. Homoclinic snaking is well understood in the context of simpler pattern forming systems such as the one-dimensional Swift-Hohenberg equation with cubic-quintic nonlinearity, whose solutions remarkably well resemble the snaking solutions of plane Couette flow. We study the structural stability of the characteristic snakes-and-ladders structure associated with homoclinic snaking for flow modifications that break symmetries of plane Couette flow. We demonstrate that wall-normal suction modifies the bifurcation structure of three-dimensional plane Couette solutions in the same way, a symmetry-breaking quadratic term modifies solutions of the one-dimensional Swift-Hohenberg equation. These modifications are related to the breaking of the discrete rotational symmetry. At large amplitudes of the symmetry-breaking wall-normal suction the connected snakes-and-ladders structure is destroyed. Previously unknown solution branches are created and can be parametrically continued to vanishing suction. This yields new localised solutions of plane Couette flow that exist in a wide range of Reynolds number.


Introduction
Invariant solutions of the Navier-Stokes equations play a key role for the dynamics of transitional shear flows (Kawahara et al. 2012). These solutions, in the form of equilibria, travelling waves and periodic orbits, have been computed for many canonical shear flows including pipe flow (Faisst & Eckhardt 2003), plane Couette flow (Gibson et al. 2009), plane Poisuille flow (Waleffe 2003), and asymptotic suction boundary layer flow (Kreilos et al. 2013). Invariant solutions are mostly found in small periodic domains, or 'minimal flow units' (Jiménez & Moin 1991). Later investigations have considered extended domains and identified localised invariant solutions that capture large scale flow patterns like turbulent spots, stripes and puffs (Avila et al. 2013;Brand & Gibson 2014;. The first family of spatially localised invariant solutions in shear flows were calculated by Schneider et al. (2010b) in plane Couette flow. This includes equilibria and travelling waves, which are equilibria in a frame of reference moving relative to the lab frame. These localised invariant solutions in plane Couette flow are specifically noteworthy because they exhibit the characteristic behaviour of homoclinic snaking (see review by Knobloch 2015) under parametric continuation (Schneider et al. 2010a). Homoclinic snaking is a process previously observed in many dissipative pattern-forming systems such as binaryfluid convection systems (Batiste et al. 2006) and optical systems (Firth et al. 2007), by which localised solutions grow in one direction while undergoing a sequence of successive saddle-node bifurcations (Woods & Champneys 1999). Homoclinic snaking manifests itself by a characteristic snakes-and-ladders structure in the bifurcation diagram (Burke & Knobloch 2007).
A well-studied one-dimensional model system which supports localized solutions that exhibit homoclinic snaking is the Swift-Hohenberg equation ∂ t u = ru+(1+∂ 2 x ) 2 u+N (u), for a real-valued function u(x) on the real axis with r the bifurcation parameter and a nonlinearity N (u) (Burke & Knobloch 2006Beck et al. 2009;Knobloch et al. 2019). Several variants of the Swift-Hohenberg equation with differing forms of the nonlinear term have been considered. Most studied are both a quadratic-cubic N = b 2 u 2 −u 3 and a cubic-quintic N = b 3 u 3 − u 5 nonlinearity, where b 2 and b 3 are adjustable parameters. For both nonlinear terms, the Swift-Hohenberg equation supports localised solutions arranged in a snaking bifurcation structure. The localised invariant solutions of the Navier-Stokes equations in plane Couette geometry share remarkably similar properties with solutions of the Swift-Hohenberg equation with the cubic-quintic nonlinearity (Schneider et al. 2010a): (i) First, the bifurcation diagram with its characteristic snakes-and-ladders structure is almost indistinguishable from that of the Swift-Hohenberg equation.
(ii) Second, the three-dimensional snaking solutions of the Navier-Stokes equations very closely resemble one-dimensional snaking solutions of Swift-Hohenberg, when threedimensional velocity fields are averaged in the downstream direction and zero-sets of downstream velocity are visualised as a function of the spanwise coordinate. The almost perfect resemblance of three-dimensional Navier-Stokes and one-dimensional Swift-Hohenberg solutions is observed for the downstream wavelength of 4π studied in Schneider et al. (2010a). Subtle modifications related to internal deformations of the three-dimensional flow field appear when the downstream wavelength changes, but the characteristic snakes-and-ladders structure remains intact (Gibson & Schneider 2016). Why localised solutions of the three-dimensional Navier-Stokes equations closely resemble solutions of the one-dimensional Swift-Hohenberg equation and are organized in a snakesand-ladders bifurcation structure is not fully understood. The detailed analysis of model systems including Swift-Hohenberg using concepts such as spatial dynamics indicates the importance of discrete symmetries for an equation to support homoclinic snaking (Champneys 1998;Burke et al. 2009;Knobloch 2015). In both the Navier-Stokes solutions and in solutions of the Swift-Hohenberg equation with cubic-quintic nonlinearity, the snakes-and-ladders bifurcation structure is composed of two pairs of intertwined snaking branches along which the solutions are invariant under discrete symmetries that are part of the equivariance group of the respective system. The snaking branches of symmetric solutions are connected by so-called rungs which emerge in symmetry-breaking pitchfork bifurcations. The rungs do not possess discrete symmetries and together with the snaking branches form the snakes-and-ladders structure.
To investigate the importance of discrete symmetries for the snaking structure of the Swift-Hohenberg equation with qubic-quintic nonlinearity, Houghton & Knobloch (2011) introduced an additional quadratic term in the equation to break the odd symmetry of the system R 2 : x → −x, u → −u. When the amplitude of this symmetry-breaking term u 2 is increased from zero, the bifurcation structure changes. One pair of snaking branches breaks into disconnected pieces while the other pair splits into two distinct curves in an amplitude versus r bifurcation diagram. Those curves are connected by zand s-shaped non-symmetric solution branches that are composed of solutions that for = 0 form rungs and the snaking branches that disappear when the symmetry is broken. In the Swift-Hohenberg equation with cubic-quintic nonlinearity, breaking a specific discrete symmetry destroys the snakes-and-ladders bifurcation structure that solutions of Navier-Stokes in plane Couette geometry remarkably well resemble. To investigate the significance of discrete symmetries of the three-dimensional Navier-Stokes equations, we break a symmetry of plane Couette and study the structural stability of the snakes-andladders structure under controlled symmetry breaking. We specifically apply wall-normal suction to break the rotational symmetry of plane Couette flow. At small amplitude, suction modifies the snakes-and-ladders structure of plane Couette flow in the same way, a quadratic term modifies the snakes-and-ladders structure of the Swift-Hohenberg equation with the cubic-quintic nonlinearity. At higher values of suction velocity, the snaking branches that remain intact at small suction velocity also break down. The solutions form separated branches which involve solutions that can be followed to zero suction velocity but are not part of the snakes-and-ladders structure of plane Couette flow without suction. We thereby identify previously unknown localised solutions of plane Couette flow that exist in a much wider range of Reynolds numbers than the snaking solutions studied previously. The structural modifications of the bifurcation structure at large amplitude of the symmetry breaking terms is not observed in the Swift-Hohenberg system where increasing the amplitude of the quadratic term transforms the snakes-and-ladders structure into a modified snakes-and-ladders structure similar to the bifurcation structure of the localised solutions of the Swift-Hohenberg equation with quadratic-cubic nonlinearity. Thus, at small suction amplitudes, suction has a similar effect as the symmetry-breaking term within the cubic-quintic Swift-Hohenberg case while at larger suction amplitudes, the bifurcation structure of 3D Navier-Stokes solutions significantly differs from that of the analogous Swift-Hohenberg problem. Both in the one-dimensional model system and in the full three-dimensional Navier-Stokes problem, the initial breakdown of the snakes-and-ladders structure can be explained in terms of symmetry breaking.
This manuscript is organised as follows: In section 2, we introduce the plane Couette system with wall-normal suction and discuss its symmetry properties. Section 3 reviews the snaking solutions and the snakes-and-ladders structure of plane Couette flow at zero suction. Section 4 presents the key observations on how wall-normal suction modifies the snaking structure. In section 5, features of the modified snaking structure are discussed and related to symmetry properties including symmetry subspaces of the flow. In the final section 6, results are summarised and an outlook for future work is provided.

System and methodology
Plane Couette flow (PCF) is the flow of a Newtonian fluid between two parallel plates at a distance of 2H which move in opposite directions with a constant relative velocity of 2U w . The laminar solution of the Navier-Stokes equations in plane Couette flow is a linear profile. We investigate the effect of wall-normal suction on the snaking solutions of plane Couette flow. The wall-normal suction modifies the laminar flow. After nondimensionalisation of the velocities with respect to half of the relative velocity of the plates, U w , and the lengths with respect to half of the gap width between the plates, H, the laminar solution takes the form where the coordinate system (x, y, z) is aligned with the streamwise (ê x ), the wallnormal (ê y ), and the spanwise directions (ê z ); U(y) is the laminar solution; V s is the nondimensionalised wall-normal suction velocity; and Decomposition of the total velocity into the laminar solution, U(y) as a base flow, and the deviation from the laminar solution, u(x, y, z, t), yields the Navier-Stokes equations in perturbative form where the Reynolds number is defined as Re = U H/ν, with ν the kinematic viscosity of the fluid. The boundary conditions for u are periodic in streamwise and spanwise directions, where L x and L z are the length and width of the computational domain in streamwise and spanwise directions, respectively. At the top and bottom plate Dirichlet conditions are satisfied. Zero pressure gradient in both streamwise and spanwise directions is imposed. Following Gibson & Schneider (2016), we present bifurcation diagrams in terms of energy dissipation of the deviation from the laminar solution, normalised by the length of the channel L x . The considered solutions are localised and thus independent of the width L z of the domain. Consequently, we do not normalise by the width L z . The streamwise averaged energy dissipation for an equilibrium or travelling wave solution equals the energy input rate I and can be expressed as in units of U 2 w . D serves as a measure for the width of a localised solution that approaches laminar flow (u = 0) and only generates nonzero dissipation in the non-laminar part of the solution.
Plane Couette flow (V s = 0) is invariant under continuous translations in x and z directions, discrete reflection in z direction and a discrete rotation of 180 • around the z axis. The rotation is equivalent to two successive discrete reflections in x and y directions. Following the conventions of Gibson et al. (2008), symmetries of the governing equations and the boundary conditions of PCF are expressed as where δx and δz are arbitrary displacements in the streamwise and spanwise directions, respectively. All products of σ z , σ xy and τ compose the symmetry, or equivariance, group of plane Couette flow (V s = 0). The equivariance group contains the discrete symmetries the snaking solutions are invariant under. These are the inversion symmetry σ xyz and the shift-reflect symmetry τ x σ z : A nonzero wall normal suction velocity breaks the rotational symmetry of the system, σ xy . Therefore, for plane Couette flow with a nonzero wall-normal suction, the equivariance group consists of all products of the translational symmetry, τ , and the reflection symmetry, σ z . Consequently, for nonzero wall-normal suction, the inversion symmetry σ xyz of PCF is broken.
The significance of symmetries for the dynamics of the flow is twofold. First, all flow fields which are invariant under the action of a symmetry contained in the equivariance group form a symmetry subspace within the system's state space that the dynamics is invariant under. That means, if an initial condition u is located in a symmetry subspace, its evolution under the governing equations remains in the same symmetry subspace. Second, if an equilibrium or travelling wave solution is invariant under a symmetry and thus located in a symmetry subspace parametric continuation will not change the symmetry of the solution but the entire continuation branch remains in the same symmetry subspace. Symmetries can only change if the branch is created or terminates in a symmetry-breaking bifurcation such as a pitchfork bifurcation.
If an equilibrium or travelling wave solution is not invariant under a symmetry contained in the equivariance group of the system the action of that symmetry on the solution creates an additional solution. Two equilibria / travelling waves that are related by that symmetry have the same global integrated properties such as dissipation. Hence, symmetry-related invariant solution branches are represented by a single curve in bifurcation diagrams presenting a typical global integrated property as a function of a control parameter.
The invariant solutions presented in this study including equilibria and travelling waves are identified numerically using the Newton-Krylov-Hookstep root-finding tools contained in Channelflow 2.0 (www.channelflow.ch, Gibson et al. (2019)). To construct the required Krylov subspace, Channelflow employs successive time integrations of the Navier-Stokes equations to evaluate the action of Jacobian by finite differencing. For integrating the Navier-Stokes equations in time, a third-order accurate semi-implicit backward differentiation scheme is used together with a pseudo-spectral Fourier-Chebychev-Fourier discretisation. In the homogeneous streamwise, and spanwise directions, we dealiase the nonlinear term according to the 2/3 rule. The computational domain has a length of L x = 4π, a width of L z = 24π, and a height of 2H = 2 and is descretised by N x = 32, N y = 35, and N z = 432 collocation points. All of the solutions are localised in spanwise directions and thus independent of the width of the computational domain. Re z y x z Figure 1: (top left) A part of the snakes-and-ladders bifurcation structure of the localised invariant solutions in plane Couette flow in a box with a length of 4π and a width of 24π. The symmetry-related solutions are visualised at the points indicated on the bifurcation diagram: The flow fields at point (i) are two symmetry-related travelling waves; point (ii) corresponds to two symmetry-related equilibrium solutions; and point (iii) represents four symmetry-related rung states. All flow fields are visualised in terms of the midplane streamwise velocity (red/blue contours) and the streamwise averaged cross-flow velocity (vector plots). In addition, the contour line of zero average streamwise velocity is presented. The red/blue streaks indicate ±0.2U w . The entire computational domain is shown. Note the localisation of all solutions in spanwise direction.
characteristic of homoclinic snaking. Two snaking curves winding upward in dissipation while oscillating in Reynolds number represent equilibrium (EQ) and travelling wave (TW) solutions, respectively. Each of the curves, both for EQs and TWs represents two symmetry-related solution branches. Pairs of equilibrium solutions and pairs of travelling waves are related by the rotational symmetry of PCF, EQ1 = σ xy EQ2, Figure 2: Symmetry relations (arrows) between solution branches of plane Couette flow without suction. The pair of travelling waves T W 1/T W 2 is located within an invariant symmetry subspace (shaded orange ellipse). Likewise, the pair of equilibria EQ1/EQ2 shares a symmetry subspace (shaded cyan ellipse). All rungs R1, R2, R3 and R4 are non-symmetric. and T W 1 = σ xy T W 2. The operation of the discrete rotational symmetry leaves global solution measures such as the norm or the streamwise averaged energy dissipation unchanged. Consequently, both symmetry-related branches appear as a single curve in the dissipation D versus Re bifurcation diagram. In the snaking region between approximately 169 < Re < 177 the dissipation of the snaking solutions increases. This is linked to the solution growing in spanwise direction while undergoing sequences of saddle-node bifurcations that create additional streaks at the solution fronts. The internal structure of the solutions remains unchanged.
The equilibrium solutions are invariant under the inversion symmetry, σ xyz which is the product of the rotational symmetry, σ xy , and the reflection symmetry, σ z . A non-trivial invariant solution which is invariant under inversion is phase locked and can neither travel in the streamwise nor the spanwise directions so that the solution is an equilibrium and not a travelling wave. In contrast, the travelling wave solutions are invariant under a different symmetry, the shift-reflect symmetry, τ x σ z which locks only the z phase of the solutions. As a result the travelling waves are free to travel in the streamwise direction. Since the equilibrium solutions are invariant under the inversion symmetry (u = σ xyz u), the rotational symmetry relation between the two equilibria EQ1 = σ xy EQ2 reduces to EQ1 = σ xy σ xyz EQ2 = σ z EQ2. Consequently, the rotational symmetry relation between the equilibrium solutions is equivalent to a mirror z-reflection symmetry relation between them. The two symmetry-related equilibria are thus mirror images with respect to the spanwise direction. For the travelling waves, their rotational symmetry relation implies that T W 1 and T W 2 travel in opposite directions but at equal phase speed.
In addition to the snaking branches, there are non-symmetric localised solutions forming rungs connecting the equilibrium branches to the travelling wave branches. Rungs are created in pitchfork bifurcations close to the saddle-node bifurcations along the snaking branches. There are four rung solution branches corresponding to each curve in the bifurcation diagram. They each connect one of the two equilibrium branches to one of the two travelling wave branches. Every pitchfork bifurcation on a specific equilibrium branch, creates two rung states which each connect this specific equilibrium branch to one of the two travelling waves. Since the rung states are non-symmetric invariant solutions which are created in symmetry-breaking pitchfork bifurcations off symmetric invariant solutions, the two simultaneously created rungs are related to each other by the symmetry of the symmetric solution branch which is broken. This is either the inversion symmetry of the equilibria R1 = σ xyz R2 and R3 = σ xyz R4 or the shift-reflect symmetry of the travelling waves. The shift-reflect symmetry relation between the two rung state solutions connecting to one travelling wave can be interpreted as only a z-reflection symmetry relation σ z between the rungs R1 = σ z R3 and R2 = σ z R4, because the shift is absorbed in the continuous translation symmetry for a solution that is not phase locked but free to travel in the streamwise direction. Hence at any Reynolds number, the four rung state solutions have the same dissipation, and the four branches of the rung states in the bifurcation diagram are represented by a single curve. The symmetry relations between all branches of the snakes-and-ladders structure in PCF are schematically summarised in figure 2 together with all relevant invariant symmetry subspaces. The two equilibria are located in the inversion symmetry subspace and they are related to each other by a z-reflection symmetry. The two travelling waves are in the subspace of the shift-reflect symmetry, and are related by the rotational symmetry, σ xy . The rung state solutions are non-symmetric and live outside these two symmetry subspaces.

Modified snakes-and-ladders bifurcation structure for non-vanishing suction
Wall-normal suction breaks the inversion symmetry of PCF and leads to modifications of all discussed solution branches. In the following we first discuss modifications of the bifurcation diagram due to suction. Then, we present how solutions along the solution branches are modified. Figure 3 shows how the snakes-and-ladders bifurcation structure is modified when wall-normal suction is applied. The two travelling wave branches undergoing snaking are symmetry-related at zero suction and thus appear as a single curve. For non-zero suction, the symmetry relating both branches is not contained in the equivariance group. Consequently, the symmetry relation vanishes and both travelling wave branches split. The two travelling wave branches remain continuous and undergo snaking but the critical Re at which saddle-node bifurcations occur alternates. Following the terminology of Knobloch (2015), we refer to the piece of a snaking curve between two forward saddlenode bifurcations including a backward saddle-node bifurcation as an oscillation in the snaking branch. In the snaking region, both travelling wave branches are composed of two alternating oscillations with different spans in Reynolds number, a narrow oscillation and a wide oscillation. Where travelling wave branch T W 1 undergoes a narrow oscillation the other one T W 2 undergoes a wide oscillation. For the next pair of oscillations the situation reverses with, T W 1 undergoing a wide, and T W 2 a narrow oscillation.

Effect of suction on the bifurcation diagram
Equilibrium branches that undergo snaking in the zero-suction case are invariant under the inversion symmetry, which is broken by suction. As a result, for non-zero suction, the continuous equilibrium solution branches vanish. Instead the branch breaks into disconnected segments. The pitchfork bifurcations that -for the zero suction casecreate rung states bifurcating from the equilibrium branches, are broken by the wallnormal suction. At non-zero suction the disconnected remainders of the equilibrium branches together with remainders of rung states form new branches that connect to The suction velocity is increased from V s = 0 (left) to V s = 10 −4 (center) and V s = 2 · 10 −4 (right). Left panel: Travelling wave (magenta) and equilibrium (green) snaking branches are shown together with rungs (black). Center and right panels: The travelling wave branch splits into T W 1 (thick red) and T W 2 (thick blue). The snaking equilibrium branch has broken into disconnected segments and together with remainders of rungs forms returning states RS and connecting states CS. Two types of RS (thin red / blue) connect to T W 1 / T W 2, respectively. The CS states (thin black) connect T W 1 and T W 2. travelling wave branches. These branches fall into two groups. Some branches emerge in a pitchfork bifurcation on one of the travelling wave branches and terminate on the same travelling wave branch in another pitchfork bifurcation. Hereafter, these branches are called returning state, RS. There are two RS branches which are related by a z-reflection symmetry, RS1 = σ z RS2. All other non-symmetric branches connect T W 1 to T W 2. These branches are, hereafter, referred to as connecting state, CS. In the continuation diagram, each CS curve represents two symmetry-related branches, CS1 = σ z CS2, which bifurcate from the travelling wave branches in pitchfork bifurcations.
The RS branches as well as the CS branches form closed bifurcation loops, as every branch has a symmetry-related counterpart starting and terminating in the same pitchfork bifurcation on T W branches. In figure 4(a) a part of the bifurcation diagram enlarging one oscillation of the travelling waves is shown. The loops can be directly observed in figure 4(b) where solutions are shown in terms of dissipation versus spanwise wave speed. The spanwise wave speed differentiates between both z-reflection symmetryrelated branches. The spanwise wave speed of the two symmetry-related branches has the same magnitude but opposite sign. Due to their shift-reflect symmetry, the spanwise wave speed of the travelling wave branches is zero.

Evolution of the flow fields on the solution branches
Each saddle-node bifurcation along the snaking branches is associated with spanwise growth of the solution as evidenced by the increasing value of dissipation. The neutral eigenmode associated with each saddle-node is localised at the fronts of the localised solution where the spatially periodic internal pattern connects to unpatterned laminar flow (Schneider et al. 2010a). When the snaking branch undergoes a saddle-node bifurcation, additional structures in form of downstream streaks together with associated pairs of counter-rotating downstream vortices is added to the solution. The solution thereby grows while its interior structure remains essentially unchanged.
Both travelling wave snaking branches T W 1 and T W 2 are invariant under shift-reflect symmetry. This symmetry is preserved and not broken by the saddle-node bifurcations along the branch. Consequently, the saddle-node bifurcations simultaneously add structures symmetrically at both fronts of the travelling wave. Figure 5 shows the evolution of the travelling waves in the snaking region for one oscillation. In this range, T W 1 undergoes a wide oscillation and T W 2 a narrow oscillation. Although the symmetry relation between T W 1 and T W 2 is broken for non-zero suction, for small values of the suction velocity the velocity fields visually still appear as if they were related by rotational symmetry.  (left), the envelope of the localised solutions shifts to the left, grows symmetrically and shifts to the right. After the sequence, the contour line remains centered at a maximum. Along the CS branch the structure shifts left, grows symmetrically and shifts left again. As a result, the contour line in the center of the periodic pattern turns from a maximum into a minimum. Figure 6 visualises the flow fields of the returning state RS and the connecting state CS branches at the points indicated in figure 4(a). Each point along both branches corresponds to two velocity fields related by z-reflection symmetry σ z . We only visualise one of the two symmetry-related velocity fields. Neither the RS nor the CS solutions are invariant under any discrete symmetry. Consequently, the growing and shrinking of the solution along the branch is not symmetric. To visualise changes in the velocity field along the branch, we overlay two successive solutions shown in figure 6. Figure 7 shows the zero contour line of the average streamwise velocity of two consecutive flow fields. The differences are small and located at the fronts. To further highlight the small differences and reveal growing and shrinking of the solution, contours of the amplitude difference of the average streamwise velocity | u 2 x | − | u 1 x | are visualised.
The returning state branch RS (thin blue line in figure 4) bifurcates in a pitchfork bifurcation from T W 2 at (6a). Towards (6c), while decreasing in Re, the solution increases in strength at its left front and weakens on its right. One may describe the localised solutions as a superposition of a periodic pattern with an envelope that supports the internal core structure and damps the exterior towards laminar flow. The change of the solution along (6a) → (6c) corresponds to a left-shift of the envelope (see figure 7 (top left panel)). From (6c) → (6e), close to the next saddle-node bifurcation, the flow field gets stronger at both sides (see figure 7 (middle left panel)). This segment of the RS branch is the remainder of an equilibrium branch at zero suction. The simultaneous and almost symmetric growth of the solution at both fronts along this segment resembles the growth of the equilibrium solution along the equilibrium branch at zero suction. Finally, from (6e) → (6g), close to the pitchfork bifurcation off T W 2, the solution increases in strength on the right and weakens on the left front ( figure 7 (bottom left panel)). The evolution along this part of the RS branch corresponds to a right-shift of the envelope. Along a full RS branch, the envelope of the localised solutions first shifts to the left, then grows symmetrically, and finally shifts to the right. After this sequence, both fronts have grown equally so that the RS branch acquires the symmetry of T W 2 it bifurcated from. The symmetry related branch RS = σ z RS bifurcates and reconnects to T W 2 together with RS but the growth along the branch is inverted: A right-shift is followed by symmetric growth and a final left-shift.
The connecting state CS branch bifurcates from T W 2 in a pitchfork bifurcation at point (6b). From (6b) → (6d) the flow fields on the CS branch strengthens on the left and weakens on the right front, which corresponds to a left-shift of the envelope (figure 7 (top right panel)). From (6d) → (6f ) the flow fields grow almost symmetrically at both fronts ( figure 7 (middle right panel)). This part of the CS branch is the remainder of an equilibrium branch at zero suction. Finally, from (6f ) → (6h) the flow fields again strengthen on the left front and weaken at the right front, corresponding to a left-shift ( figure 7 (bottom right panel)). At point (6h) the CS branch connects to the T W 1 branch in a pitchfork bifurcation. Along the CS branch the envelope first moves to the left, then grows symmetrically, and finally moves to the left again. As a result of this growth sequence involving a net shift to the left, the centrally located low speed streak (a blue streak in figure 6) is replaced by the neighbouring high speed streak at the left (a red streak in figure 6), which now sits at the center of the localised solution. A centrally located high speed streak is characteristic of T W 1. Consequently, the CS branch starts on T W 2 and terminates on T W 1. The symmetry-related branch CS = σ z CS bifurcates from T W 2 and, as CS connects to T W 1. However, along CS the solution exhibits a net shift to the right until the solution is no longer centered at a low speed streak but on the next high-speed streak at the right. Since continuous shifts in spanwise direction are part of the system's equivariance group, despite net shifts in opposite directions, both connecting state branches CS and CS connect to the same T W 1 branch.

Discussion
Non-vanishing suction velocity V s causes splitting of the T W 1 and T W 2 branches and creates new CS and RS branches from the remainders of EQ and rung branches. In this section we relate these modifications of the snakes-and-ladders bifurcation structure due to the suction to the breaking of symmetries of PCF with zero suction. Moreover, we discuss additional modifications of the bifurcation structure observed at large amplitudes of the suction velocity V s .

Snaking solutions and symmetry subspaces in the presence of suction
Wall normal suction breaks the rotational symmetry of plane Couette flow. As a result, any symmetry subspace or symmetry relation originating from rotational symmetry is not present for non-zero suction velocity. Consequently, for non-zero suction, there is no inversion symmetry subspace, σ xyz and states that are inversion-symmetric equilibria at zero suction travel in the streamwise and the spanwise directions. Moreover, the rotational symmetry relation between the two travelling wave branches and the inversion symmetry relation between rung states vanishes. Figure 8 schematically indicates the configuration of symmetry subspaces and symmetry relations of states in the presence of wall-normal suction. Since the symmetry relation between the two travelling wave branches is broken by the wall-normal suction, we expect a separation of the travelling wave branches T W 1 and T W 2 in the bifurcation diagram showing dissipation versus Reynolds number. Both solution branches that are equilibria at zero suction remain symmetry related, EQ1 = σ z EQ2 which implies equal dissipation so that they remain represented by a single curve. At zero suction all four rung states are related by two different symmetry transformations. For V s = 0, the inversion symmetry relation σ xyz between R1 and R2 and between R3 and R4, the rotational symmetry relation σ xy between T W 1 and T W 2 and the inversion symmetry subspace, all present for V s = 0 (see figure 2) are broken.
One of those two symmetries is broken by suction so that the four rungs split into two groups of two symmetry-related branches each, R1 = σ z R3 and R2 = σ z R4. Thus two separate curves in the bifurcation diagram represent two symmetry-related solution branches each. Rungs and equilibrium branches at zero suction transform into connecting states CS and returning states RS at non-zero suction.

Front growth controls bulk velocity and oscillation width
The mean pressure gradient in both spanwise and downstream direction is imposed to be zero and each solution selects its bulk velocity U bulk , the y − z averaged streamwise velocity. Figure 9 shows the variation of the bulk velocity of the snaking solutions both for PCF with zero suction and for a suction velocity of V s = 10 −4 . At zero suction, the EQ branch is invariant under σ xyz , implying zero bulk velocity. The bulk velocity of the travelling wave solutions periodically varies around zero as the solution undergoes snaking. The magnitude of the U bulk oscillations is independent of the spatial extent of the solution which suggests the non-zero bulk velocity is generated by the fronts, at which high-and low-speed streaks are growing symmetrically Gibson & Schneider (2016). While low-speed streaks are created, the bulk velocity becomes negative, and when high-speed streaks grow at the fronts the bulk velocity becomes positive. When suction is applied, the oscillations in bulk velocity are overlaid by an additional negative component that is linearly proportional in the dissipation. The linear dependence on dissipation and thereby size of the solution suggests that the linearly growing component of the bulk velocity is generated by an unchanging internal structure of the solution while the oscillations are due to growth at the fronts.
The front-mediated oscillations in bulk velocity relative to the linear trend are correlated with the width of snaking oscillations of the travelling wave branches. Excursions to higher bulk velocity are observed when the branch undergoes a wide oscillation, i.e. a wider range of Reynolds numbers (see figure 4). Likewise, narrow oscillations are observed when the bulk velocity approaches local minima. Maxima of the bulk velocity are linked to high-speed streaks growing at the front while for minima in the bulk velocity, growing low-speed streaks are observed. Since T W 1 remains related to T W 2 by the approximate though broken rotational symmetry σ xy , a growth of high-speed streaks along the T W 1 branch implies the growth of low-speed streaks along T W 2 and vice versa. The growth of high-and low-speed streaks of approximately symmetry-related travelling wave branches thus explains why wide and narrow oscillations of T W 1 and T W 2 both alternate and moreover occur such that a wide oscillation of T W 1 coincides with a narrow one of T W 2 and vice versa.

Splitting of the travelling waves
In PCF without suction, both travelling wave solution branches, T W 1 and T W 2, are represented by a single curve in the bifurcation diagram in terms of D versus Re. Wallnormal suction breaks the rotational symmetry of the system σ xy and results in splitting of the travelling wave solution branches (see figure 3). A subset of the bifurcation diagram is enlarged in figure 10(a), where dissipation D as a function of Re is presented for T W 1 and T W 2. The branches are shown for three values of the wall-normal suction velocity, V s = 0, V s = 10 −4 and V s = 2·10 −4 . For non-zero suction, branches split so that T W 1 and T W 2 curves are symmetrically located around the zero-suction curve, with equal distance but on opposing sides. The distance between both split curves grows linearly with the suction velocity. To rationalise this splitting behaviour we consider modifications of the solution of the governing equations at leading, namely linear, order in V s . At zero wallnormal suction V s = 0, the laminar flow solution is a linear profile U(y) = yê x . Non-zero wall suction modifies the travelling wave solutions T W 1 and T W 2. When decomposing the velocity into a laminar flow and deviations, this modification affects both the laminar base flow U and the deviation from the modified laminar base flow u, as we will describe in this section.

Relation between travelling waves for inverted suction velocity
Solutions cannot only be followed from zero to positive suction velocity V s = +|V s | but also to negative values V s = −|V s |. The latter physically implies blowing through the wall. Since the rotational symmetry σ xy , relating T W 1 to T W 2 at V s = 0 and broken Consequently, the modification of the laminar solution with respect to the laminar solution of PCF for zero suction velocity U = U − U Vs=0 is whereŪ is a function of Reynolds number and the wall-normal coordinate y only.Ū is independent of the suction velocity so that the laminar flow modification U is linearly proportional to the suction velocity V s .

The modification of the deviation from the laminar solution
The travelling wave solutions T W 1 and T W 2 are equilibria in a moving frame of reference that translates at their specific wave speeds c x in theê x direction with respect to the lab frame. The travelling wave solutions thus satisfy the condition where c x is the wave speed associated with the total flow field u t = U(y)+u(x, y, z). The solution of this equation (u t , p, c x ) may be decomposed into three parts: the solution for zero suction, a modification of the laminar solution due to suction, and the modification of the deviation from the laminar solution: where (u 0 , p 0 , c x0 ) denotes the total solution for zero suction, U is the modification of the laminar solution, and (u , p , c x ) is the modification of the deviation from the laminar solution. For small suction velocity, quadratic interactions among U , u and c xêx resulting from the inertial term are neglected. At leading order, the condition for a travelling wave solution at small suction thus reads Rearranging using the base flow modifications at linear order in V s (Equation 5.2), yields where the suction velocity has been isolated. Only the primed variables depend on V s so that the entire left-hand side is proportional to V s . The equation determines (u , p , c x ) as a function of V s . Since the operators acting on the primed variables are linear (right-hand side), the solution (u , p , c x ) is linearly proportional to the suction velocity V s . Dissipation of T W 1 with wall-normal suction velocity V s can be expressed as where D Vs=0 is the dissipation of both T W 1 and T W 2 at zero suction velocity. According to relation 5.1, the dissipation of T W 2 with suction velocity V s is equal to the dissipation of T W 1 for inverted suction velocity (−V s ). As a result, the dissipation of T W 2 follows from the same expression when the suction velocity in Equation 5.4 is inverted (V s → −V s ). Equation 5.4 implies u (−V s ) = −u (V s ). Consequently, the dissipation of T W 2 with wall-normal suction velocity V s is given by Suction thus changes the dissipation of travelling wave solutions T W 1 and T W 2 relative to the zero suction case by amounts that are equal in magnitude but opposite in sign. Equation 5.4 implies that the entire solution (u , p , c x ) is linear in V s so that also modifications of the wave speed c x are proportional to the suction velocity. Equation 5.1 implies that the modification of the wave speed of T W 2 with wall-normal suction velocity V s is equal to minus the modification of the wave speed of T W 1 for inverted suction velocity (−V s ). Together with the linear dependence of wave speed modifications c x on V s , this implies that the modifications of the wave speed are equal for both T W 1 and T W 2. Figure 10 velocity, V s = 0, V s = 10 −4 and V s = 2 · 10 −4 . This data confirms that small values of the wall-normal suction velocity V s change the wave speeds of both T W 1 and T W 2 by the same amount. The amount is linearly proportional to the suction velocity.

Alternating bifurcations of CS and RS off travelling wave branches
Both travelling wave branches, T W 1 and T W 2 undergo alternating wide and narrow oscillations in the dissipation versus Reynolds number bifurcation diagram. Along a wide oscillation a pair of high-speed streaks is added at the fronts of the solution, while along a narrow oscillation a pair of low speed streaks is added. After an entire period consisting of a successive wide and narrow oscillation, the fronts connecting the internal periodic pattern of the travelling waves to laminar flow recover their original structure and the solution has grown by two additional pairs of high-and low-speed streaks. The solution thus only grows at the fronts while its internal structure remains unchanged. As a result of this front-driven growth process, the critical Reynolds numbers of corresponding saddlenode bifurcations along successive wide or narrow oscillations line up. The spatial growth of the solution indicated by increasing values of dissipation D continues periodically.
The growth of the travelling wave branches is driven by saddle-node bifurcations with neutral modes that are localised at the solution fronts (Schneider et al. 2010a). Likewise, the pitchfork bifurcations close to the saddle-nodes that create the RS and CS have neutral modes acting only at the fronts. As a result, the entire bifurcation structure, including the bifurcating RS and CS branches, repeats periodically after a pair of oscillations. Moreover, the front structure of T W 1 reflects that of the T W 2 branch shifted by half a period. Since all bifurcations are driven at the fronts, the equivalent front structure of both travelling wave branches implies that all branches bifurcating off T W 1 have corresponding branches bifurcating off T W 2. Consequently, the entire bifurcation structure of T W 1 including bifurcations off the branch corresponds to that of the T W 2 branch shifted by half a period.
As discussed in section 4.2, CS and RS branches are characterised by the sequence of directions in which the envelope of the localised solutions is shifted as one marches along the branch. Starting from the pitchfork bifurcations off a travelling wave at lower D where the branches emerge towards the pitchfork at higher D where the branches terminate on a travelling wave, both CS and RS encounter two envelope shifts. The connecting state CS branch is characterised by two shifts to the same directions either to the right (+ê z ) or to the left (−ê z ). Two shifts in opposing directions (left after right or vice versa) characterise the returning state RS branch. Figure 11 enlarges a part of the bifurcation diagram showing both T W 1 and T W 2 branches, two successive CS branches and one RS branch connected to T W 2. Note that each CS or RS curve in the bifurcation diagram represents a pair of z-reflection-related branches of CS or RS. In the figure, we indicate the shifts of the envelope for those CS and RS branches that inherit segments from EQ1 for V s = 0 (from cs 2 → cs 3 , from rs 2 → rs 3 and from cs 6 → cs 7 ). The symmetry-related counterparts of these CS and RS branches contain remainders from EQ2. Due to the symmetry relation, any shift of the branch shown in the figure corresponds to a shift in the opposite direction for the (not shown) symmetry-related branch. If one CS branch twice shifts to the left (from cs 1 → cs 4 ), the CS for the next higher dissipation (from cs 5 → cs 8 ) twice shifts in the opposite, here right, direction.
For V s = 0, the two rung state branches that bifurcate off an equilibrium branch in a pitchfork bifurcation are related by σ xyz . For non-zero but small values of V s , the remainders of these two rung state branches form segments of the CS and RS branches (in the figure, cs 3 → cs 4 / rs 1 → rs 2 and rs 3 → rs 4 / cs 5 → cs 6 ). The solutions on these remainders of the rung state branches remain visually almost symmetry-related at small V s . As a result of the approximate inversion symmetry, along these two segments of the CS and RS branches, the directions of the shifts of the envelope has to be identical. Consequently, the shifts along segments of the RS branches are slaved to the shifts of the CS branches along the corresponding segments by the broken inversion symmetry relation of the rung states. The non-symmetric solution branch between two CS branches in the bifurcation diagram must thus be characterised by shifts in two opposing directions and is therefore an RS branch. Symmetry thus implies that CS and RS branches alternate in the bifurcation diagram.

Snaking breakdown
The effect of suction on plane Couette snaking solutions was studied for small suction velocities up to V s = 2 · 10 −4 . At larger suction velocities, the bifurcation structure fundamentally changes. Figure 12(a) shows the bifurcation diagram of T W 1, T W 2, and CS branches for V s = 6 · 10 −4 . The range in Re for which RS branches exist shrinks with increasing V s (see figure 3) and at V s = 6 · 10 −4 , RS branches are not computed. At V s = 6·10 −4 , the snaking travelling wave solution branches, T W 1 and T W 2, have broken into multiple disconnected travelling wave branches that reach high Reynolds numbers far beyond the snaking range at V s = 0. The two symmetry-related CS branches no longer connect T W 1 and T W 2. Instead each CS branch emerging in pitchfork bifurcations off each disconnected travelling wave branch now extends to high Reynolds numbers but no longer appears to terminate in a second pitchfork bifurcation.
Parametric continuation down in V s starting from all branches identified at V s = 6 · 10 −4 including those reaching high Reynolds numbers yields previously unknown localised solution branches of PCF without suction. Figure 12(b) shows the bifurcation diagram with the characteristic snakes-and-ladders structure together with these newly-  figure 13. At V s = 6 · 10 −4 , the snaking branches of T W 1 and T W 2 are broken, and T W 1, T W 2 and CS branches are merged with solution branches that for V s = 0 are separated from the snakes-and-ladders bifurcation structure and extend to high Reynolds numbers.
found localised solution branches for V s = 0. The non-snaking localised solutions include branches that for V s = 6 · 10 −4 are merged with the T W 1 branch (hereafter referred to as T W 1-connected), the T W 2 branch (hereafter referred to as T W 2-connected), and the CS branches (hereafter referred to as CS-connected). The solutions of T W 1-and T W 2connected branches are invariant under the action of the shift-reflect symmetry, τ x σ z . The solutions of the CS-connected branches are not invariant under any discrete symmetry of the system. Instead there are pairs of CS-connected solution branches related by σ z . Each pair bifurcates from either a T W 1-or T W 2-connected branch in a pitchfork bifurcation. Snaking solutions in PCF at V s = 0 only exist in a small Reynolds number range of approximately 169 < Re < 177. The localised non-snaking solution branches, T W 1-, T W 2-and CS-connected, span a much wider range in Reynolds number (see figure 12(b)). Figure 13 visualises the flow fields of T W 1-, T W 2-and two successive CS-connected branches at the points indicated in figure 12(b). The internal periodic part and the front structures of the lower branch T W 1-and T W 2-connected solutions (panels a and b); the upper branch T W 1-and T W 2-connected solutions (panels e and f ); and the two successive CS-connected solutions (panels c and d, and their symmetry-related  figure  12(b). The flow fields in the left panels are located on the T W 1-connected branch (panels a and e) and on the CS-connected branch that bifurcates from the T W 1-connected branch (panel c). The right panels visualise the flow fields that are located on the T W 2connected branch (panels b and f ) and on the CS-connected branch that bifurcates from the T W 2-connected branch (panel d). The flow fields are visualised as in figure 1. counterparts) appear identical. The solutions mainly differ in the number of high-and low-speed streaks. From lower to upper branches of both T W 1-and T W 2-connected two high-speed streaks grow at the solution fronts. The T W 1-connected solution is centered at a high-speed streak while in the center of the T W 2-connected solution a low-speed streak is located. Application of the rotational symmetry σ xy to T W 1-, T W 2-and CSconnected solutions creates symmetry-related invariant solutions. However, the branches of these solutions do not merge with the snakes-and-ladders bifurcation structures in the presence of wall-normal suction.

Conclusion
We investigate the structural stability of the snakes-and-ladders bifurcation structures of invariant solutions of the three-dimensional Navier-Stokes equations in plane Couette flow. Salewski et al. (2019) show that adding a Coriolis force in PCF (V s = 0) that maintains the equivariance group, preserves the snakes-and-ladders bifurcation structure. Here, we show that applying a non-vanishing suction velocity that breaks the rotational symmetry of PCF modifies the snakes-and-ladders bifurcation structure. For non-zero but small suction velocity, the curve representing the branches of both travelling waves T W 1 and T W 2 = σ xy T W 1 (at V s = 0) in a bifurcation diagram showing dissipation D as a function of Re split in two different snaking curves with alternating span of the oscillations in Reynolds number. The both equilibrium branches EQ1 and EQ2 = σ z EQ1 (at V s = 0) break up. The pitchfork bifurcations of the equilibrium branches that create rungs at V s = 0 are broken. The result are new solution branches formed from remainders of both the the broken equilibrium branches and the rungs. The returning state branches RS connect one of the travelling wave branches to itself while the connecting state branches CS connect T W 1 and T W 2. Specific features of the bifurcation diagram including the symmetric splitting of the T W 1 and T W 2 branches and the ordering of RS and CS follow from symmetry arguments.
At small but non-vanishing suction velocity, the modifications of the snakes-andladders bifurcation structure of three-dimensional solutions of Navier-Stokes equations in plane Couette flow are analogous to modifications of the snakes-and-ladders bifurcation structure observed within the one-dimensional Swift-Hohenberg equation with the qubicquintic nonlinearity (SHE35), when a quadratic term is added. Introducing a quadratic term u 2 with amplitude in the SHE35 breaks the odd symmetry of this system, R 2 : x → −x, u → −u (Houghton & Knobloch 2011). Breaking R 2 in the SHE35 modifies the snakes-and-ladders bifurcation structure in the following way. For non-vanishing both solution branches that are invariant under the action of the even symmetry, R 1 : x → −x, u → u split into two distinct curves with alternating span of the oscillations in r in an amplitude versus control parameter r bifurcation diagram. This is analogous to the splitting of the T W 1 and T W 2 in the Navier-Stokes problem. As the EQ branches in Navier-Stokes, the continuous snaking branches of solutions that are invariant under R 2 at = 0 break up into disconnected segments for non-zero . Together with remainders of rung states at = 0 the remainders of these solutions form disconnected s-and z-shaped branches that connect to the split snaking solution branches that are invariant under R 1 . The s-and the z-shaped solution branches in the modified snakes-and-ladders bifurcation structure of the one-dimensional SHE35 with the symmetry-breaking quadratic term thus behave analogous to the CS and RS solution branches in the modified snakes-and-ladders bifurcation structure of Navier-Stokes equations in PCF in the presence of non-vanishing but small wall-normal suction velocity. When the amplitude of the quadratic term is increased towards large values, the snakes-and-ladders bifurcation structure in the SHE35 eventually transforms into a modified bifurcation structure that resembles the snakesand-ladders structure of the Swift-Hohenberg with the quadratic-cubic nonlinearity. This is not observed in the Navier-Stokes problem, where large suction causes the modified snakes-and-ladders bifurcation structure to disintegrate into disconnected branches of localised solutions. In conclusion, at small amplitudes, wall-normal suction in plane Couette flow and the symmetry-breaking quadratic term in the SHE35 yield analogous effects on the snakes-and-ladders bifurcation structure. At larger amplitudes, however, the effect of wall-normal suction in plane Couette flow and the symmetry-breaking term in SHE35 affect the bifurcation structure in significantly different ways.
Following the solution branches generated by suction back to regular plane Couette flow with V s = 0, additional disconnected branches of spatially localised invariant solutions of PCF are identified. Disconnected branches of solutions with varying width but matching front structure exists. The solutions thus share similarities with the snaking solutions but do not undergo homoclinic snaking themselves. While the snaking solutions only exist in a limited range of Reynolds number, the non-snaking solution branches span a much larger range of Reynolds numbers extending to those Reynolds numbers in which localised turbulent patterns are observed (Barkley & Tuckerman 2005). Consequently, the non-snaking localised invariant solutions identified here might be more relevant for supporting localised transitional turbulence than the previously known snaking branches. The dynamical relevance of the non-snaking localised solutions should be investigated in the future.