## Abstract

We report the computation of a family of travelling wave solutions of pipe flow up to *Re*=75 000. As in all lower branch solutions, streaks and rolls feature prominently in these solutions. For large *Re*, these solutions develop a critical layer away from the wall. Although the solutions are linearly unstable, the two unstable eigenvalues approach 0 as *Re*→∞ at rates given by *Re*^{−0.41} and *Re*^{−0.87}; surprisingly, the solutions become more stable as the flow becomes less viscous. The formation of the critical layer and other aspects of the *Re*→∞ limit could be universal to lower branch solutions of shear flows. We give implementation details of the GMRES-hookstep and Arnoldi iterations used for computing these solutions and their spectra, while pointing out the new aspects of our method.

## 1. Introduction

In this paper, we look at a lower branch travelling wave solution in the *Re*→∞ limit. The travelling wave that we chose to compute has an asymmetric arrangement of streaks, with two fast streaks located preferentially on one side of the pipe. Schneider *et al*. (2007*b*) found that states with such an asymmetry arise in direct numerical simulations of transition to turbulence. Pringle & Kerswell (2007) computed such a travelling wave using a bifurcation point of a mirror-symmetric family around *Re*=1000. Our computations of the same travelling wave go up to *Re*=75 000 and help elucidate aspects of the *Re*→∞ asymptotic limit.

The fast streaks near the wall are the most prominent and stable structures in lower branch travelling wave solutions of pipe flow (Faisst & Eckhardt 2003; Wedin & Kerswell 2004). The fast streaks are regions in a circular section where the streamwise velocity significantly exceeds the laminar value. The fast and slow streaks can form different patterns. The pattern that characterizes some of the computed solutions is an invariance with respect to rotation about the pipe axis by 2*π*/*m*, where *m*=2, 3, 4, …. The rolls, which correspond to positive and negative streamwise vorticity, form complementary patterns. Although the computed solutions use periodic boundary condition in the axial direction and very short pipes, they do pick structures that transitional pipe flow tends to develop (Hof *et al*. 2004; Willis & Kerswell 2008). The data analysis techniques used to extract these patterns are set up to pick patterns with rotational symmetry (Eckhardt *et al*. 2007; Schneider *et al*. 2007*a*; Willis & Kerswell 2008). The streak pattern of the asymmetric travelling wave does not have any *m*-fold rotational symmetry as evident from figure 1.

Wang *et al*. (2007; also see Waleffe 2003) showed that the *Re*→∞ limit of a symmetric lower branch solution of plane Couette flow is characterized by a number of features. The streaks remain *O*(1), but the magnitude of the rolls and of the fundamental and higher streamwise modes decrease algebraically with *Re*. The scaling exponents for the rolls and the fundamental streamwise mode of the asymmetric travelling wave are −1.08 and −0.97, which may be compared with −1 and −0.9 for the symmetric solution of plane Couette flow. Higher streamwise modes decrease even faster.

The most important consequence of these scalings is the development of a critical layer away from the circular boundary of the pipe. The theory of Wang *et al*. (2007) successfully identifies the critical curve as given by *w*_{0}(*r*, *θ*)=*c*_{z}, where *w*_{0} is the *z*-averaged streamwise velocity and *c*_{z} is the wave speed in the *z* direction. The fundamental component of the radial velocity is concentrated in a region around the critical curve and drops off to zero away from that region. We find that the size of the region decreases at the rate *Re*^{−0.32} as *Re* increases, which compares well with the rate of *Re*^{−1/3} derived by Wang *et al*. (2007) using formal arguments. The exponents for the rates at which the sizes of the regions decrease with *Re* are different for the fundamental mode of the streamwise velocity and the mean streamwise vorticity. These are found to be −0.26 and −0.23 in §4. These exponents present a challenge to asymptotic theory.

At the end of §4, we suggest that it might be useful to calculate the analogue of the critical curve for puffs. Puffs have a well-defined extent and travel down the pipe with a well-defined speed. The analogue of the critical curve would be a surface, embedded inside the puff, on all points of which the streamwise velocity equals the speed of the puff. Such a surface could be helpful in elucidating the structure of the puff.

The Newton equations for solving a nonlinear system can sometimes be solved efficiently in a Krylov subspace (Brown & Saad 1990; Sancheź *et al*. 2004). We point out two new aspects of the extensions to the Newton–Krylov procedure introduced by Viswanath (2007). The first novelty is the formulation of the Newton equations. In the case of pipe flow, the formulation allows for translation of the velocity field along the pipe axis or rotation around the pipe axis. The second novelty is the GMRES-hookstep combination explained in §5.

For large *Re*, the lower branch asymmetric travelling wave looks very different from both the laminar solution of pipe flow and the sort of turbulence that is typically observed at such *Re*. Unlike the laminar solution, the travelling wave develops streaks, for instance. Unlike fully developed turbulence, there is no rapid decay of correlations. The form of the asymmetric travelling wave is nearly independent of the *z* direction at high *Re*. Thus, one may ask if the lower branch solutions are relevant for high *Re* turbulence and if they can be realized in the laboratory. The answer to the first question is probably no. The second question is a difficult challenge to experiment. That the computations are restricted to small pipes is less of an issue for high *Re* owing to the scaling of the streamwise modes mentioned above and discussed in §3.

## 2. Preliminary data

The asymmetric travelling waves were computed at a number of values of *Re* in the range 1500≤*Re*≤75 000. Some basic data are summarized in table 1. The choice of units and boundary conditions follows that of Faisst & Eckhardt (2004). The pipe radius is chosen as the unit of length. The unit of velocity is equal to the centreline velocity of the Hagen–Poiseuille laminar flow. The Reynolds number is based on the pipe radius, centreline velocity of the Hagen–Poiseuille laminar flow and the kinematic viscosity. The boundary condition is no slip at the pipe wall and periodic in the axial direction. The mass flux of the flow, which is fixed at 0.5, drives the flow. The pipe length or period is 2*πΛ*. We took *Λ*=1/1.44, but this choice has no special significance in the *Re*→∞ limit.

The quantities *L*, *M*, *N* listed in table 1 parametrize the spatial grid used to represent the velocity field. The spatial coordinate system *r*, *θ*, *z* was cylindrical, with *u*, *v*, *w* being the three components of velocity, respectively. The three components of vorticity are denoted as *ξ*, *η* and *ζ*. The radial component of the velocity field ** u** is represented as(2.1)with the discretization using 2

*M*and 2

*N*Fourier points along

*θ*and

*z*, respectively. The coefficients are even in

*r*for

*m*odd, and odd for

*m*even. Each is represented using its values at the Chebyshev points

*r*=cos(

*πi*/

*L*), 0≤

*i*≤(

*L*−1)/2. Note that

*L*is always odd. The vorticity component

*ξ*has an analogous representation. As the velocity field has zero divergence, the entire velocity field can be recovered using

*u*,

*ξ*, and , where and are averages of

*v*and

*w*with respect to both

*θ*and

*z*. After setting the modes with |

*m*|=

*M*or |

*n*|=

*N*to zero, we are left with (

*L*−2)+((2

*N*−1)(2

*M*−1)−1)(

*L*−3)/2 independent degrees of freedom.

In terms of the modes, the boundary conditions become and . The constant mass flux condition implies a pressure gradient along *z* that can change from instant to instant for an evolving flow.

The wave speed of the travelling wave is given by *c*_{z}. To find each travelling wave, one solves for a velocity field *u*_{0} such that ** u**(

*r*,

*θ*,

*z*,

*t*)=

*u*_{0}(

*r*,

*θ*,

*z*−

*c*

_{z}

*t*) is a solution of the Navier–Stokes equation. The artificial parameter

*T*, which occurs in table 1, arises in the solution procedure and its meaning is explained in §5.

The rate of energy dissipation per unit mass is given by 2*D*/*Re*, where *D* is the integral ofover the volume of the pipe. The rate of energy input per unit mass is given by 2*I*/*Re*, wherewith *p* being pressure and with the integral being over the volume of the pipe. The friction coefficient (Wedin & Kerswell 2004) is the same as *I*, but with a different normalization. *D* and *I* are normalized to be 1 for the laminar flow. From table 1, we see that *D*=*I* for all the travelling waves in agreement with energy conservation. Kinetic energy, denoted KE in table 1, is also normalized to evaluate to 1 for laminar flow.

The Navier–Stokes equation for pipe flow, with periodic boundary along *z*, is unchanged by the shift-reflect transformation. The shift-reflect transformation reflects the velocity field about the plane *θ*=0 or *θ*=*π*, and shifts it along *z* by half a pipe length. All the asymmetric travelling waves have only two unstable eigenvalues in the shift-reflection symmetric subspace. Those are given as *λ*_{1} and *λ*_{2} in table 1. Section 6 has a discussion of the spectrum of the travelling waves.

## 3. Scaling of modal kinetic energies

Figure 2*a* shows the variation of the KE in various modes as a function of *Re*. To find the KE for the *n*=1 streamwise mode, we form Fourier expansions of type (2.1) for *v* and *w* as well. The volume integral for KE is computed by setting all modes with *n*≠±1 equal to zero. The kinetic energies of the other streamwise modes are computed in a similar manner.

The KE of the rolls is obtained using *n*=0 mode only, but the *w* component is set to zero. Retaining only the *n*=0 modes is equivalent to averaging the velocity field with respect to *z*. The *z*-averaged *w* corresponds to streaks.

As evident from figure 2*a*, the magnitudes of the modes decrease with *Re* algebraically and are proportional to *Re*^{e} for high *Re* and a suitable exponent *e*. The exponents for the rolls, *n*=1, *n*=2 and *n*=3 obtained using *Re*≥8000 were −1.08, −0.97, −1.35 and −1.92, respectively. For comparison, the exponents for rolls and the *n*=1 mode are −1 and −0.9 for the symmetric lower branch solution of plane Couette flow (Wang *et al*. 2007).

Figure 3*b* shows that the wavespeed *c*_{z} increases with *Re*. An application of Wynn's *ρ*-algorithm (Wynn 1956) shows the limit of *c*_{z} as *Re*→∞ to be 0.88. The speed of the asymmetric travelling wave is nearly twice the speed of puffs in transitional pipe flow. In our units, the speed of the puff is approximately 0.45 around *Re*=2000 (Peixinho & Mullin 2006).

From figure 3, we conclude that the streaks converge as *Re*→∞ and that the plots in that figure are a good approximation to the limit. Those plots differ quite a bit from the plot at *Re*=3000 in figure 2, with the position of the two high-speed streaks being much more to the left of the pipe at *Re*=3000.

## 4. The critical layer

The Fourier expansion of *u* (2.1) can be rewritten aswhere the asterisk denotes complex conjugation. Similar expansions can be formed for *v*, *w* and the vorticity components. To illustrate the critical layer, we will begin by looking at |*u*_{1}|.

Wang *et al*. (2007) derived the equation *w*_{0}(*r*, *θ*)=*c*_{z} for the critical curve. The critical curve is shown as a thick black curve in figure 4*a*. It is closer to the centre of the pipe than to the pipe wall. The contour lines of |*u*_{1}| are all nestled around the critical curve. In particular, the contour lines occur as two groups near the indentation at the left of the critical curve. This compares well with fig. 3 of Wang *et al*. (2007). Figure 4*b* shows that |*u*_{1}| takes its maximum value on or very close to the critical curve and falls off rapidly away from the critical curve. The first two plots of figure 4 give a good idea of how |*u*_{1}| varies inside the unit circle. The critical region is a band around the critical curve where most of the variation of |*u*_{1}| and certain other quantities is concentrated. The band need not be of uniform width.

Figure 4*c* shows contour plots of *ζ*_{0}. The regions where *ζ*_{0} is positive or negative agree very well with the position of the rolls. Counter-rotating vortices are a well-known feature of lower branch solutions and of small perturbations of the laminar flow that trigger turbulence. Like rolls and streamwise modes, the scaling of whose magnitudes with *Re* is shown in figure 2, the magnitude of *ζ*_{0} also decreases with *Re*.

From figure 4*c*, it is evident that most of the variation of |*u*_{1}| is in a region around the critical curve. Similar plots can be produced for |*w*_{1}| or *ζ*_{0}. In such plots the peaks become notably sharper as *Re* increases.

The purpose of figure 5 is to estimate the rate at which the contour curves, such as those in figure 4*a*,*c*, approach the critical curve as *Re*→∞. For each value of *Re*, a specific contour curve is picked. For |*u*_{1}|, |*w*_{1}| and *ζ*_{0}, the chosen contour curve is for half their maximums. We pick the point on the contour curve that is farthest from the critical curve and plot its distance against *Re*. Such plots are a good way to measure the thickness of the critical region. They follow the convention in which the width of a Gaussian density function is measured at half its maximum.

Fits using *Re*≥8000 show that the thickness scales as *Re*^{−0.32}, *Re*^{−0.26} and *Re*^{−0.23} for |*u*_{1}|, |*w*_{1}| and *ζ*_{0}, respectively. The exponents do not change appreciably if fits are made by dropping the data points with smaller *Re*.

Perhaps the main achievement of Wang *et al*. (2007) is to give a formula for the critical curve. In the context of pipe flow, the critical curve is the set of all points (*r*, *θ*) such that *w*_{0}(*r*, *θ*)=*c*_{z}. We have used that formula throughout this section. Their calculations apply directly to |*u*_{1}| and |*v*_{1}|, and predict that the contour curves of those quantities will approach the critical curve at a rate given by *Re*^{−1/3}. The exponent that we found for |*u*_{1}|, which came in at −0.32, is in excellent agreement with that prediction. The exponents for |*w*_{1}| and *ζ*_{0} indicate that the contour curves of those quantities concentrate more slowly on the critical curve than those of |*u*_{1}|. A more refined theory is probably needed to explain those exponents.

The thickness of the critical layer is highly unlikely to be uniform around the critical curve. The manner in which the thickness varies along the critical curve appears worthy of investigation. It appears that the variation of the thickness could be related to the structure of the rolls. Even at low *Re*, such as *Re*=1500, contour plots still show that structures tend to develop around the critical curve. This motivates a suggestion that will end this section.

Puffs are structures observed in transitional pipe flow that have a well-defined extent. They travel down the pipe with a well-defined speed. It could be interesting to calculate the surface formed by all points of the puff whose streamwise velocity equals the speed at which the puff moves down the pipe. Such a surface would be the analogue of the critical curve for a lower branch travelling wave.

## 5. Implementation of GMRES-hookstep and Arnoldi iterations

In §2, we pointed out that the velocity field for pipe flow with suitable boundary conditions can be recovered from , , *u* and *ξ*. If we pack the information in those variables into a single column vector ** x** with real components, it is possible to recover the entire velocity field given

*x*.

*X*(

*t*;

*x*) is the column vector that results from allowing the flow to evolve for time

*t*. To compute

*X*(

*t*;

*x*), a velocity field is constructed starting from

*x*and then allowed to evolve for time

*t*using a direct numerical simulation code.

*X*(

*t*;

*x*) is constructed from the final velocity field. We have generally used Runge–Kutta methods with constant step sizes (except for the last step) to compute

*X*(

*t*;

*x*). The reason is that the discretized flow is then a dynamical system that is smooth and close to the Navier–Stokes flow. Adaptive time stepping strategies introduce non-smoothness and imply that the discretized flow is no longer a dynamical system.

The methods for computing travelling waves and other solutions that will be described depend upon the shear flow mainly in the computation of *X*(*t*; *x*). The other dependence is in the definition of the translation operators. Given the Fourier representation (2.1) of *u*(*r*, *θ*, *z*), the representation after a translation along the axis and a rotation about the axis is given by(5.1)We use linear operators defined by(5.2)to effect the translation and the rotation in (5.1). In particular,The definition of the linear operators depends upon the shear flow. The definition of the linear operators for plane Couette flow is identical to that for pipe Poiseuille flow (Viswanath 2007). The operators can be made to act on a vector ** x** that encodes a velocity field in an obvious way, by making them act on each component of the velocity field. Then encodes a translated and rotated velocity field. Expressing the translation and rotation of a velocity field using makes it possible to differentiate with respect to

*s*

_{θ}and

*s*

_{z}while deriving the Newton equations.

Given the ability to compute *X*(*x*; *t*) and the linear operators of (5.2), the numerical methods described in this section need to know nothing more about the shear flow. Determining the exact dimension of the vector ** x** can be a little tricky because one needs to eliminate Fourier coefficients that are conjugates of certain others and so on (Viswanath 2007). It is unlikely that one may leave out some essential components as this error will become manifest when trying to construct the velocity field from

*x*. It is more likely that

*x*may end up having duplicates. In principle, that would make some of the matrices that occur later singular. In practice, the effect of having duplicates in

*x*is probably to introduce some error without being disastrous.

A big part of the numerical method for computing travelling waves, relative periodic orbits and other solutions that will now be described, are the well-known GMRES and Arnoldi iterations. Trefethen & Bau (1997) give a lucid account of these methods and more importantly their convergence properties. Pointers to the original literature can be found in the end notes of their book or in many other well-known textbooks of numerical linear algebra.

### (a) GMRES-hookstep iteration

A relative periodic orbit is a solution of the Navier–Stokes equation where the initial velocity field evolves for time *T*, which is the period, to reach a certain final state. In the case of pipe flow, it must be possible to translate the final velocity field along the axis and then rotate it to get back the initial velocity. If *x*_{0} encodes the initial velocity field,(5.3)where *s*_{θ} and *s*_{z} are shifts in the azimuthal and streamwise directions, respectively. To find a relative periodic orbit, one must solve for *x*_{0}, *s*_{θ}, *s*_{z} and the period *T* such that the nonlinear equation (5.3) is satisfied.

A relative periodic orbit is the most general object that our method can find. Periodic orbits are a special case where *s*_{θ}=*s*_{z}=0. Travelling waves are a special case where *T* is fixed to be a small but not too small number. A travelling wave will satisfy (5.3) for any *T*>0 and suitably chosen *s*_{θ}, *s*_{z}. In general, the solution of (5.3) could be a relative periodic orbit that is not a travelling wave. *T* is chosen small enough to make it likely that the solution of (5.3) is a travelling wave, although it is not important to have a small *T* if we already know that the initial guess for *x*_{0} is near a travelling wave. An equilibrium or steady solution is also a special case of a relative periodic orbit. The reason for treating travelling waves as special cases of relative periodic orbits is explained at the end of this section.

Suppose , *s*_{x}, *s*_{z}, *T* is an initial guess to a solution of (5.3) and that(5.4)Linearizing one gets the following Newton equations (Viswanath 2007):(5.5)In the above system, *I* is the identity matrix whose dimension equals that of ; and *f*(*x*) is such that d*x*/d*t*=*f*(*x*) is the spatially discretized Navier–Stokes equation written in terms of the vector ** x** that encodes the discretized velocity field. The code for evaluating

*f*(

*x*) can be extracted from a direct numerical simulation code with a little work. One can also approximate

*f*(

*x*) as (

*X*(

*h*,

*x*)−

*x*)/

*h*, where

*h*is small. We have not tried approximating

*f*(

*x*) using differences, but it is probably fine to do so. The last three rows of the linear system (5.5) correspond to phase conditions (Viswanath 2007).

To find a relative periodic orbit, one step of the Newton iteration would be to solve (5.5) for the *δ*s and add those corrections to the initial guess. To find a travelling wave, (5.5) must be modified by dropping the last row and the last column because *T* is fixed. If the travelling wave has the shift-reflect symmetry, as the travelling wave family studied in this paper does, then *s*_{θ}=0, because rotation around the pipe axis breaks that symmetry. In such a case, we must drop the first and the third of the last three columns, and likewise with the rows. To find an equilibrium solution, we must drop the last three columns and rows. All the special cases of a relative periodic orbit mentioned above can be dealt with in this manner. In each case, we denote the resulting linear system as *AΔ*=*b*.

To solve such a linear system using a Krylov subspace method like GMRES, it is not necessary to invert *A* nor is it even necessary to form *A* explicitly. It is enough if *A* can be applied to vectors. The only difficulty in applying *A* to a vector arises in calculatingwhere ** c** is a column vector of the same dimension as . That quantity can be calculated using differences as(5.6)where

*ϵ*is chosen such that . The choice of the norm will be discussed shortly. Even when is nearly equal to

*y*

_{0}, which is defined by (5.4), it is important not to substitute for

*y*

_{0}in (5.6).

The GMRES iteration for solving *AΔ*=*b* finds an orthonormal matrix *Q*_{k} at the *k*th stage such that (Trefethen & Bau 1997). In implementing this step, it may be best to use the square root of the KE of the vector field that *x* encodes as the norm over *x*. At the *k*th stage GMRES would solve the least-squares problem , where *y*∈*R*^{k} and *e*_{1} is the *k*+1 dimensional vector with a 1 at the top followed by 0s. The approximation to *Δ* at that stage would be *Δ*_{k}=*Q*_{k}*y*. We do not attempt to solve the Newton equation this way, however. The Newton equation is useful only if the solution *Δ* is tiny enough that the linearization that led to the Newton equation is valid. That is often not the case because the initial guesses are typically not so accurate. A well-known solution is to minimize subject to the constraint ‖*Δ*‖≤*δ*, where *δ* has to be chosen small enough that the linearization within that radius is valid (Dennis & Schnabel 1996). The resulting step is called the hookstep (Dennis & Schnabel 1996).

We approximate the hookstep using GMRES as follows. To find *Δ*_{δ,k} that approximates the true hookstep *Δ*_{δ} we solve the minimization problem(5.7)subject to the constraint ‖*y*‖≤*δ*. That minimization can be solved using the singular value decomposition (Dennis & Schnabel 1996; Golub & van Loan 1996). Let be a reduced singular value decomposition (*V*′ is the transpose of the real unitary matrix *V*). Let . If the diagonal entries of the diagonal matrix *D* are *d*_{i}, is found using , 1≤*i*≤*k*, where either *μ*>0 is such that ‖*q*‖=*δ* or *μ*=0 if that allows ‖*q*‖≤*δ*. Finding *μ* is an easy one-dimensional root-finding problem. The solution of (5.7) is *y*=*Vq* and the GMRES hookstep is *Δ*_{δ,k}=*Q*_{k}*y*.

To complete the description of the GMRES-hookstep method, we have to describe the choice of *k*, or the stopping criterion for finding a *Δ*_{δ,k} that approximates *Δ*_{δ}, and also describe how *δ* is updated every time a new Newton system (5.5) is formed. There is a natural stopping criterion for GMRES without the constraint ‖*y*‖≤*δ*. That is because the relative residual error at the end of *k* iterations can be easily found as . For GMRES hookstep, we have no practical way of knowing how close is to . Thus, there is no way to assess the quality of *Δ*_{δ,k}. The stopping criterion in our implementation is to pick a *k* that is large enough to ensure *r*_{k}≤0.01. In other words, we stop when the GMRES iterate *Δ*_{k} is an acceptable substitute for the true solution of *AΔ*=*b* believing then that the Krylov subspace matrix *Q*_{k} has enough column vectors to ensure that *Δ*_{δ,k} is an acceptable substitute for *Δ*_{δ}. There is no theoretical support for this stopping criterion, but it works very well in practice.

The choice of *δ* follows standard trust-region prescriptions (Dennis & Schnabel 1996). The choice for *δ* for the very first GMRES-hookstep iteration can be anything that looks reasonable. To assess the quality of a *δ*, we take as the error in the initial guess. Once *Δ*_{δ,k} is computed, we update to , where dim is the dimension of and the subscripting of *Δ*_{δ,k} follows Matlab notation. The quantities *s*_{θ}, *s*_{z} and *T* are also updated, if applicable. The linearization used to find *Δ*_{δ,k} predicts that the reduction in error in going from to should be about . If the prediction is very good *δ* can be increased, and if it is bad *δ* must be decreased and a new GMRES hookstep must be computed. This completes the description of the GMRES-hookstep method for solving (5.3), each iteration of which begins with a guess for *x*_{0} and for the shifts and the period, forms the Newton system (5.5), uses that Newton system to find *Δ*_{δ,k}, checks if *δ* is acceptably small and then uses *Δ*_{δ,k} to form a better guess. The iterations can be stopped if the error as measured by is less than the relative error due to spatial discretization of the velocity field.

It is surprising that the method for computing *Δ*_{δ,k} is a new contribution considering it is quite a natural thing to do. In an early paper on the use of Krylov subspaces for globally convergent modifications of Newton's method, Brown & Saad (1990) formulated a minimization problem ((4.2) of their paper) and called it the model trust-region problem. The solution to that problem is theoretically equivalent to *Δ*_{δ,k}. The equivalence is similar to that between GMRES and ORTHODIR, which pre-dated GMRES, with our formulation being more direct. We have described a practical method for finding *Δ*_{δ,k} with a criterion for choosing *k*. We were not able to find implementations of GMRES hookstep in the literature, although one may exist that we were not able to track down.

Like the work of Brown & Saad (1990) much of the later literature deals with the dog-leg and other strategies; for instance see Luksan & Vlcek (1997). The dog-leg is an approximation to the hookstep that is made up of only the gradient direction and the Newton step (Dennis & Schnabel 1996). It is preferred over the hookstep mainly because its computation does not require the singular value decomposition. Since the hookstep moves away from the Newton step smoothly, one may suggest that the Krylov subspace approximates the hookstep better than the gradient. The dog-leg is also much more complicated to implement within a Krylov subspace than the computation of *Δ*_{δ,k} described here. Having to compute the singular value decomposition is not a problem because the way the Newton system (5.5) is set up means that *k* is small (being approx. 150 at most but more typically approx. 50). Since the dog-leg is only an approximation to the hookstep, and is in fact harder to implement within a Krylov subspace, we see no reason to prefer it over the GMRES-hookstep method.

### (b) Arnoldi iteration

Ignoring spatial discretization errors, the eigenvalues *μ*_{i} of the matrix(5.8)are the eigenvalues of the corresponding relative periodic or periodic solution. If *x*_{0} encodes the velocity field of a travelling wave or a relative periodic solution, then *μ*_{i}=exp(*λ*_{i}*T*) where *λ*_{i} are the eigenvalues of the travelling wave or the equilibrium solution.

The matrix (5.8) will be dense and large, but it can be applied to vectors as in (5.6). The Arnoldi iteration forms *Q*_{k}, *Q*_{k+1} and *H*_{k+1,k} like GMRES, with the one difference being that the starting vector ** b** is arbitrary. We usually take

*x*

_{0}as the starting vector but either rotate and translate it or add some noise to ensure that it does not have the shift-reflect symmetry. In the case of both pipe and channel flows, the laminar solution must be subtracted from

*x*

_{0}to get the right boundary conditions. If

*H*

_{k}is the matrix obtained by dropping the last row of

*H*

_{k+1,k}, and

*H*

_{k}

*y*=

*μy*, then

*μ*is an approximation for an eigenvalue of (5.8) with

*Q*

_{ky}being an approximation for the corresponding eigenvector.

The approximations *μ* and *y* must be checked for correctness. If *μ* is real, one only has to apply the matrix (5.8) with *Q*_{k}*y* and verify if the resulting vector has the right amplitude and direction; if *μ* is complex one has to apply the matrix to the real part of *Q*_{k}*y*. In figures 6 and 7, we accept an eigenvalue if the result of applying the matrix has an error in direction that is less than 1 degree and the error in amplitude is less than 1 per cent. Most eigenvalues and eigenvectors are much more accurate than that, and it is reasonable to expect the eigenvalues to be more accurate than the eigenvectors.

If *x*_{0} is the initial velocity field of a travelling wave, its wave speeds are given by *c*_{θ}=(*s*_{θ}+2*πp*)/*T* and *c*_{z}=(*s*_{z}+2*πΛq*)/*T*, where *p* and *q* are integers. The values of *p* and *q* are found by advancing the initial velocity field by an amount of time that is not too large, and then translating and rotating the final velocity field to see which values of *p*, *q* imply the best match to the initial velocity field. In the case of the asymmetric travelling wave, *c*_{θ}=*s*_{θ}=0 owing to symmetry and care is needed for determining *c*_{z} at high *Re* because there is very little energy in the streamwise modes with *n*≠0.

In the case of travelling waves, there is a delicate numerical point that arises in passing from a complex eigenvalue *μ* of (5.8) to an eigenvalue *λ*=log(*μ*)/*T* of the travelling wave. Figure 6*c* shows the *λ*s that correspond to the *μ*s in figure 6*a*. The imaginary part of the complex log is not unique, and to determine it for the *λ*s one has to in effect determine the rate of rotation of the real part of the eigenvector in the space spanned by the real and imaginary parts. If the column *c* is the real part of the eigenvector, the matrix–vector productfor *t* not too large will give the correct rate of rotation. To find that matrix–vector product, we can again use differences as in (5.6) but there are two mathematically equivalent ways to do so. The first way is to use the quotient difference(5.9)where is determined using the same direct numerical simulation code and the same time step used to compute , and the second way is to use(5.10)We must use (5.9), although (5.10) involves less work. The numerical errors in using the quotient difference (5.10) will be intolerably high.

The eigenvalues in figure 6*a*,*b* are mostly inside the unit circle and stable. Most of the eigenvalues of the matrix (5.8) are stable owing to the dissipation term in the Navier–Stokes equation. For a demonstration of the effect of the dissipation term, note that the stable eigenvalues for *Re*=1500 are closer to the circle than those of *Re*=15 000, even though the computation at *Re*=15 000 uses a larger *T* (see table 1) which brings the stable eigenvalues closer to the centre.

Setting up the eigenvalue problem for travelling waves using direct numerical simulation and the matrix (5.8) may seem contrived owing to the need to choose an artificial parameter *T* and the need to use direct numerical simulation. Contrived it may be, but the contrivance does serve a purpose. Without it we will have a spectrum that will look like the one in figure 6*c*, but with a lot of eigenvalues with very large and negative real parts not shown there. For a matrix with such a spectrum, the Arnoldi iteration will not work well because it will be forced to chase the eigenvalues with large and negative real parts. With matrix (5.8), those eigenvalues move very close to 0, and the extremal part of the spectrum that is approximated well is also the interesting part of the spectrum for stability considerations.

## 6. Spectrum of lower branch travelling waves as *Re*→∞

The Arnoldi iterations for travelling waves at various *Re* were carried out using *k*=150. For *Re*=1500, 122 out of 150 eigenvalues of *H*_{k} turned out to be correct. For *Re*=15 000 as well, 122 out of the 150 eigenvalues were correct, but the time of integration was higher with *T*=15.

At *Re*=1500, the asymmetric travelling wave has two real unstable eigenvalues, whose eigenvectors are invariant under shift reflection. Those two eigenvalues persist as *Re*→∞. Surprisingly, those two eigenvalues approach 0 as *Re*→∞. Figure 7*b* shows that the rate of decrease of those eigenvalues is algebraic. The most unstable eigenvalue approaches 0 at the rate *Re*^{−0.41}. The other eigenvalue approaches 0 at the faster rate *Re*^{−0.87}. For the symmetric lower branch solution of plane Couette flow, there is just one unstable eigenvalue and that decreases at the rate *Re*^{−0.46} or *Re*^{−0.48} (Wang *et al*. 2007; Viswanath 2008). Figure 7*a*,*b* shows that the spectrum as a whole approaches the imaginary axis as *Re* increases.

In addition to the two real unstable eigenvalues with eigenvectors in the symmetric subspace, there is an unstable complex pair at *Re*=1500 which can be seen in figure 6*a*. That pair moves inside the circle as *Re* increases. At *Re*=3000 and 5000, there is a third real and weakly unstable eigenvalue. For *Re*≥8000, there seem to be only two unstable eigenvalues, and both of those have eigenvectors that are invariant under shift reflection.

## 7. Conclusion

We have demonstrated the existence of a critical layer in the *Re*→∞ limit for a family of lower branch travelling waves. The theory of Wang *et al*. (2007) gives the right formula for the critical curve. The scaling of the size of the critical region for |*u*_{1}| is in excellent agreement with their theory. Further development of the asymptotic theory appears necessary to explain the scaling of the size of the critical regions for |*w*_{1}| and *ζ*_{0}. Comparison with a family of lower branch equilibrium solutions of plane Couette flow suggests that the formation of the critical layer and many of its properties could be universal to all lower branch solutions of shear flows as *Re*→∞.

Certain parts of puffs, which are structures observed in transitional pipe flow, are characterized by streaks and rolls (Hof *et al*. 2004; Willis & Kerswell 2008). We have suggested that the critical surface of a puff could be helpful in visualizing its structure. In particular, the arrangement of rolls and streaks could be correlated with the shape of the critical surface.

In §5, we have given a detailed account of the GMRES-hookstep iteration for computing relative periodic solutions, travelling waves, periodic solutions and equilibria for shear flows. Our account emphasizes the implementation aspects of GMRES-hookstep and of the Arnoldi iteration, which is used for finding eigenvalues. Together with the derivation of the Newton equations (Viswanath 2007), this account is sufficiently detailed to enable implementation of these iterations.

## Acknowledgments

The author thanks the mathematics department of the Indian Institute of Science, Bangalore, for its hospitality and support. The author thanks F. Waleffe and J. F. Gibson for their helpful discussions, and W. R. Morrow for catching a bad typo. This work was partially supported by NSF grants DMS-0407110 and DMS-0715510.

## Footnotes

One contribution of 10 to a Theme Issue ‘Turbulence transition in pipe flow: 125th anniversary of the publication of Reynolds' paper’.

- © 2008 The Royal Society