# Transition to turbulence in a shear flow

## I Introduction

In many flows the transition to turbulence proceeds via a sequence of bifurcations to flows of ever increasing spatial and temporal complexity. Analytical and experimental efforts in particular on layers of fluid heated from below  and fluids between rotating concentric cylinders  have lead to the identification and verification of several routes to turbulence, which typically involve a transition from a structureless laminar state to a stationary spatially modulated one and then to more complicated states in secondary and higher bifurcations.

Transitions in shear flows do not seem to follow this pattern . Typically, a transition to a turbulent state can be induced for sufficiently large Reynolds number with finite amplitude perturbations, just as in a subcritical bifurcation. However, in the most spectacular cases of plane Couette flow between parallel plates and Hagen-Poiseuille flow in a pipe , there is no linear instability of the laminar profile for any finite Reynolds number that could give rise to a subcritical bifurcation. The turbulent state seems to be high dimensional immediately, without clear temporal or spatial patterns (unlike the rolls in Rayleigh-Bénard flow). And the transition seems to depend sensitively on the initial conditions. Based on these characteristic features it has been argued that a novel kind of transition to turbulence different from the well-known three low-dimensional ones is at work .

Recent activity has focussed on three features of this transition: the non-normality of the linear eigenvalue problem , the occurence of new stationary states without instability of the linear profile  and the fractal properties of the lifetime landscape of perturbations as a function of amplitude and Reynolds number . The non-normality of the linear stability problem implies that even in the absence of exponentially growing eigenstates perturbations can first grow in amplitude before decaying since the eigenvectors are not orthogonal. During the decay other perturbations could be amplified, giving rise to a noise sustained turbulence . The amplification could also cause random fluctuations to grow to a size where the nonlinear terms can no longer be neglected . Then the dynamics including the nonlinear terms could belong to a new asymptotic state, different from the laminar profile, perhaps a turbulent attractor. Presumably, this attractor would be built around stationary or periodic solution. Here, the observation of tertiary structures  comes in since they could form the basis for the turbulent state. Finally, the observation of fractality in the lifetime distribution suggests that the turbulent state is not an attractor but rather a repeller: Infinite lifetimes occur only along the stable manifolds of the repeller, all other initial conditions will eventually decay. Permanent turbulence would thus correspond to noise induced excitations onto a repeller.

In plane Couette flow some of the features described above have been identified, but only with extensive numerical effort . The aim of the present work is to present a simple model that is based on the Navier-Stokes equation and captures the essential elements of the transition. It is motivated in part by the desire to obtain a numerically more accessible model which perhaps will provide as much insight into the transition as the Lorenz model  for the case of fluids heated from below (presumably at the price of similar shortcomings). The two and three degree of freedom models proposed by various groups (and reviewed in ) to study the effects of non-normality mock some features of the Navier-Stokes equations considered essential by their inventors but they are not derived in some systematic way from the Navier Stokes equation. The model used here differs from the one proposed by Waleffe  in the selection of modes.

Attempts to built models for shear flows using Fourier modes immediately reveal an intrinsic difficulty: In the case of fluids heated from below the nonlinearity arises from the coupling of the temperature gradient to the flow field so that two wave vectors, $`𝐤`$ and $`2𝐤`$, suffice to obtain nonlinear couplings. In shear flows, the nonlinearity has to come from the coupling of the flow field with itself through the advection term $`(𝐮)𝐮`$. This imposes rather strong constraints on the wave vectors. At least three wave vectors satisfying the triangle relation $`𝐤_1+𝐤_2+𝐤_3=0`$ are required to collect a contribution from the advection term. A minimal model thus has at least six complex variables. Three of these decay monotonically to zero, leaving three for a nontrivial dynamics. In the subspaces investigated (B.E., unpublished), the most complex behaviour found is a perturbed pitchfork bifurcation, which may be seen as a precursor of the observed dynamics: for Reynolds numbers below a critical value, there is only one stable state. Above that value a pair of stable and unstable states is born in a saddle-node bifurcation. The stable state can be excited through perturbations of sufficient amplitude. The basins of attraction of the two stable states are intermingled, but the boundaries are smooth.

Thus more wave vectors are needed and they have to couple in a nontrivial manner to sustain permanent dynamics. The specific set of modes used is discussed in section II. It is motivated by boundary conditions for the laminar profile and the observation that wave vectors pointing to the vortices of hexagons satisfy the triangle conditions in a most symmetrical manner. Other than that the selected vectors are a matter of trial and error. In the end we arrive at a model with 19 real amplitudes, two force terms and 212 quadratic couplings. Without driving and damping the dynamics is energy conserving, as would be the corresponding Euler equation (suitably truncated). Moreover, the perturbation amplitudes can be put together to give complete flow fields. Thus the model has a somewhat larger number of degrees of freedom, but the dynamics should provide a realistic approximation to shear flows.

The outline of the paper is as follows. In section II we present the model, in particular the selected wave vectors, the equations of motion and a discussion of symmetries. In section III we focus on the dynamical properties of initial perturbations as a function of amplitude and Reynolds number. In section IV we discuss the stationary states, their bifurcations and their stability properties. We conclude in section V with a summary and a few final remarks.

## II The model shear flow

Ideal parallel shear flows have infinite lateral extension. Both in experiment and theory this cannot be realized. We therefore follow the numerical tradition and chose periodic boundary conditions in the flow and neutral direction. The flow is confined by parallel walls a distance $`d`$ apart. A convenient way to built a low dimensional model is to use a Galerkin approximation. Solid boundaries would require the vanishing of all velocity components and complicated Galerkin functions where all the couplings can only be calculated numerically. However, under the assumption that here as well as in many other situations the details of the boundary conditions effect the results only quantitatively but not qualitatively, we can adopt free-free boundary conditions on the walls and use simple trigonometric functions as basis for the Galerkin expansion. Similarly, the nature of the driving (pressure, boundary conditions or volume force) should not be essential so that we take a volume force proportional to some basis function (or a linear combination thereof). This still leaves plenty of free parameters to be fixed below.

### A Galerkin approximation

We expand the velocity field in Fourier modes,

$$𝐮(𝐱,t)=\underset{𝐤}{}𝐮(𝐤,t)e^{i𝐤𝐱}.$$
(1)   
Incompressibility demands

$$𝐮(𝐤,t)𝐤=0.$$
(2)   
The Navier-Stokes equation for the amplitudes $`𝐮(𝐤,t)`$ becomes

$`_t𝐮(𝐤,t)=`$ $``$ $`ip_𝐤𝐤i{\displaystyle \underset{𝐩+𝐪=𝐤}{}}\left(𝐮(𝐩,t)𝐪\right)𝐮(𝐪,t)`$ (3)   
$``$ $`\nu 𝐤^2𝐮(𝐤,t)+f_𝐤`$ (4)   
where $`p_𝐤`$ are the Fourier components of the pressure (divided by the density), $`\nu `$ is the kinematic viscosity and $`f_𝐤`$ are the Fourier components of the volume force sustaining the laminar profile.

There are three constraints on the components $`𝐮(𝐤)`$: incompressibility (2), reality of the velocity field,

$$𝐮(𝐤)=𝐮(𝐤)^{}$$
(5)   
and the boundary conditions that the flow is limited by two parallel, impenetrable plates. The ensuing requirement $`u_z(x,y,z)=0`$ at $`z=0`$ and $`z=d`$ (where $`d`$ is the separation between plates) is most easily implemented through periodicity in $`z`$ and the mirror symmetry

$$\left(\begin{array}{c}u_x\\ u_y\\ u_z\end{array}\right)(x,y,z)=\left(\begin{array}{c}u_x\\ u_y\\ u_z\end{array}\right)(x,y,z),$$
(6)   
which in Fourier space requires

$$\left(\begin{array}{c}u_x\\ u_y\\ u_z\end{array}\right)(k_x,k_y,k_z)=\left(\begin{array}{c}u_x^{}\\ u_y^{}\\ u_z^{}\end{array}\right)(k_x,k_y,k_z).$$
(7)   
This is not sufficient to fix the coefficients: the dynamics also has to stay in the relevant subspace, and thus the time derivatives have to satisfy similar requirements.

### B The wave vectors

The choice of wave vectors is motivated by the geometry of the flow and the aim to include nonlinear couplings. The basic flow shall be a flow in $`y`$-direction, neutral in the $`x`$-direction and sheared in the $`z`$-direction. Thus we take the first three wave vectors in $`z`$-direction,

$$𝐤_1=\left(\begin{array}{c}0\\ 0\\ 1\end{array}\right),𝐤_2=\left(\begin{array}{c}0\\ 0\\ 2\end{array}\right),𝐤_3=\left(\begin{array}{c}0\\ 0\\ 3\end{array}\right).$$
(8)   
The negative vectors $`𝐤_i`$ also belong to the set but will not be numbered explicitely. In these units, the periodicity in the $`z`$-direction is $`2\pi `$, so that the separation between the plates is $`d=\pi `$ because of the mirror symmetry (6). The amplitude $`𝐮(𝐤_1)`$ will carry the laminar profile and $`𝐮(𝐤_3)`$ can be excited as a modification to the laminar profile. $`𝐤_2`$ is needed to provide couplings through the nonlinear term. These three vectors satisfy a triangle identity $`𝐤_1+𝐤_2𝐤_3=0`$, but the nonlinear term vanishes since they are parallel.

The next set of wave vectors contains modulations in the flow and neutral direction,

$$𝐤_4=\left(\begin{array}{c}1\\ 0\\ 0\end{array}\right),𝐤_5=\left(\begin{array}{c}1/2\\ \sqrt{3}/2\\ 0\end{array}\right),𝐤_6=\left(\begin{array}{c}1/2\\ \sqrt{3}/2\\ 0\end{array}\right).$$
(9)   
Together with $`𝐤_i`$ they form a regular hexagon, so that they provide nontrivial couplings in the nonlinear term. The periodicity in flow direction is $`4\pi /\sqrt{3}`$, in the neutral direction it is $`4\pi `$.

Finally, this hexagon is lifted upwards with $`𝐤_1`$ and $`𝐤_2`$ to form the remaining 12 vectors,

$`𝐤_7`$ $`=`$ $`𝐤_1+𝐤_4𝐤_8=𝐤_1+𝐤_5𝐤_9=𝐤_1+𝐤_6`$ (10)   
$`𝐤_{10}`$ $`=`$ $`𝐤_1𝐤_4𝐤_{11}=𝐤_1𝐤_5𝐤_{12}=𝐤_1𝐤_6`$ (11)   
$`𝐤_{13}`$ $`=`$ $`𝐤_2+𝐤_4𝐤_{14}=𝐤_2+𝐤_5𝐤_{15}=𝐤_2+𝐤_6`$ (12)   
$`𝐤_{16}`$ $`=`$ $`𝐤_2𝐤_4𝐤_{17}=𝐤_2𝐤_5𝐤_{18}=𝐤_2𝐤_6.`$ (13)   
The full set $`𝐤_i`$, $`i=1\mathrm{}18`$ is shown in Fig. 1.

The Fourier amplitudes $`𝐮(𝐤_i)`$ have to be orthogonal to $`𝐤_i`$ because of incompressibility (2). If they are expanded in basis vectors perpendicular to $`𝐤_i`$, the pressure drops out of the equations and need not be calculated. We therefore chose normalized basis vectors

$`𝐧(𝐤_i)`$ $`=`$ $`({\displaystyle \frac{k_xk_z}{k_x^2+k_y^2}},{\displaystyle \frac{k_yk_z}{k_x^2+k_y^2}},1)^T/\sqrt{1+k_z^2/(k_x^2+k_y^2)}`$ (14)   
$`𝐦(𝐤_i)`$ $`=`$ $`(k_y,k_x,0)^T/\sqrt{k_x^2+k_y^2}`$ (15)   
so that $`𝐧`$, $`𝐦`$ and $`𝐤`$ form an orthogonal set of basis vectors. For the negative vectors $`𝐤_i`$ we chose the basis vectors $`𝐧(𝐤_i)=𝐧(𝐤_i)`$ and $`𝐦(𝐤_i)=𝐦(𝐤_i)`$. If the $`x`$ and $`y`$ components of $`𝐤`$ vanish, the above definitions are singular and replaced by

$$𝐧=(1,0,0)^T𝐦=(0,1,0)^T.$$
(16)   
The amplitudes of the velocity amplitude are now expanded as

$$𝐮(𝐤_i,t)=\alpha (𝐤_i,t)𝐧(𝐤_i)+\beta (𝐤_i,t)𝐦(𝐤_i).$$
(17)   
The impenetrable plates impose further constraints on the $`\alpha (𝐤_i)`$ and $`\beta (𝐤_i)`$. For $`i=1`$, $`2`$ and $`3`$ the wave vector has no components in the $`x`$\- and $`y`$-directions, so that $`\alpha `$ and $`\beta `$ have to be real. For $`i=4`$, $`5`$ and $`6`$ the velocity field cannot have any components in the $`z`$-direction, hence $`\alpha =0`$. The remaining wave vectors $`𝐤_i`$ and $`𝐤_i`$ with $`i=7,\mathrm{},18`$, a total of 24, divide up into six groups of 4 vectors each,

$$𝐤=(k_x,k_y,k_z),𝐤^{}=(k_x,k_y,k_z),𝐤\text{and}𝐤^{}.$$
(18)   
The groups are formed by the vectors and their negatives in the pairs with indices (7,10), (8,11), (9,12), (13,16), (14,17) and (15,18). The amplitudes of the vectors in the sets are related by

$`\alpha (𝐤)`$ $`=`$ $`\alpha (𝐤)^{}=\alpha (𝐤^{})^{}`$ (19)   
$`\beta (𝐤)`$ $`=`$ $`\beta (𝐤)^{}=\beta (𝐤^{})^{}.`$ (20)   
Thus the full model has $`6+6+6\times 4=36`$ real amplitudes. Restricting the flow by a point symmetry around $`𝐱_0=(0,0,\pi /2)^T`$ eliminates the contributions from $`𝐤_2`$ and some other components, resulting in a 19-dimensional subspace with nontrivial dynamics and the following amplitudes:

$`\alpha (𝐤_1)`$ $`=`$ $`y_1\beta (𝐤_1)=y_2`$ (21)   
$`\alpha (𝐤_3)`$ $`=`$ $`y_3\beta (𝐤_3)=y_4`$ (22)   
$`\beta (𝐤_4)`$ $`=`$ $`iy_5\beta (𝐤_5)=iy_6\beta (𝐤_6)=iy_7`$ (23)   
$`\alpha (𝐤_7)`$ $`=`$ $`y_8\beta (𝐤_7)=y_9`$ (24)   
$`\alpha (𝐤_8)`$ $`=`$ $`y_{10}\beta (𝐤_8)=y_{11}`$ (25)   
$`\alpha (𝐤_9)`$ $`=`$ $`y_{12}\beta (𝐤_9)=y_{13}`$ (26)   
$`\alpha (𝐤_{13})`$ $`=`$ $`iy_{14}\beta (𝐤_{13})=iy_{15}`$ (27)   
$`\alpha (𝐤_{14})`$ $`=`$ $`iy_{16}\beta (𝐤_{14})=iy_{17}`$ (28)   
$`\alpha (𝐤_{15})`$ $`=`$ $`iy_{18}\beta (𝐤_{15})=iy_{19};`$ (29)   
components not listed vanish or are related to the given ones by the boundary conditions (20). A complete listing of the flow fields $`𝐮_i`$ associated with the coefficients $`y_i`$ such that $`𝐮=_iy_i𝐮_i`$ as well as of the equations of motion are available from the authors.

### C The equations of motion

In this 19-dimensional subspace $`y_1\mathrm{}y_{19}`$ the equations of motion are of the form

$$\dot{y}_i=\underset{j,k}{}A_{ijk}y_jy_k\nu K_iy_i+f_i.$$
(30)   
Of the driving force all components but $`f_2`$ and $`f_4`$ vanish. Moreover, if the $`f`$’s are taken to be proportional to $`\nu `$, the resulting laminar profile has an amplitude independent of viscosity (and thus Reynolds number). These components give rise to a laminar profile that is a superposition of a $`\mathrm{cos}(z)`$ profile (from $`f_2`$) and a $`\mathrm{cos}(3z)`$ profile (from $`f_4`$). This allows us to approximate the first two terms of the Fourier expansion of a linear profile with velocity $`u_y=\pm 1`$ at the walls,

$$𝐮_0=\frac{8}{\pi ^2}(\mathrm{cos}z+\frac{1}{9}\mathrm{cos}3z)𝐞_y.$$
(31)   
that can be obtained with a driving $`f_2=4\nu /\pi ^2`$ and $`f_4=4\nu /9\pi ^2`$ (see Fig. 2).

The nonlinear interactions in the Navier-Stokes equation conserve the energy $`E=\frac{1}{2}𝑑V𝐮^2`$. In the 19-dimensional subspace, the corresponding quadratic form is

$$E=V\left(\underset{i=1}{\overset{7}{}}y_i^2+2\underset{i=8}{\overset{19}{}}y_i^2\right).$$
(32)   
The above equations conserve this form without driving and dissipation. With dissipation but still without driving, the time derivative is negative definite, indicating a monotonic decay of energy to zero.

Finally, we define the Reynolds number using the wall velocity of the linear profile, $`u_0=1`$, the half width of the gap, $`D=d/2=\pi /2`$ and the viscosity $`\nu `$,

$$Re=u_0D/\nu =\pi /2\nu .$$
(33)   
The other geometric parameters are a period $`4\pi /\sqrt{3}`$ in flow direction and $`4\pi `$ perpendicular to it.

### D Symmetries

We achieved the impenetrability of the plates by requiring the mirror symmetry:

$$\left(\begin{array}{c}u_x\\ u_y\\ u_z\end{array}\right)(x,y,z)=\left(\begin{array}{c}u_x\\ u_y\\ u_z\end{array}\right)(x,y,z).$$
(34)   
The reduction from 36 to 19 modes was achieved by restricting the dynamics to a subspace where the flow has the point symmetry around $`𝐱_0=(0,0,\pi /2)^T`$, a point in the middle of the shear layer,

$$\left(\begin{array}{c}u_x\\ u_y\\ u_z\end{array}\right)(x,y,z+\pi /2)=\left(\begin{array}{c}u_x\\ u_y\\ u_z\end{array}\right)(x,y,z+\pi /2).$$
(35)   
In addition, there are further symmetries that can be used to reduce the phase space. There is a reflection on the $`y`$-$`z`$-plane,

$$T_1:\left(\begin{array}{c}u_x\\ u_y\\ u_z\end{array}\right)(x,y,z)\left(\begin{array}{c}u_x\\ u_y\\ u_z\end{array}\right)(x,y,z).$$
(36)   
and two shifts by half a lattice spacing,

$`T_2:`$ $`\left(\begin{array}{c}u_x\\ u_y\\ u_z\end{array}\right)(x,y,z)\left(\begin{array}{c}u_x\\ u_y\\ u_z\end{array}\right)(x+2\pi ,y,z)`$ (37)   
$`T_3:`$ $`\left(\begin{array}{c}u_x\\ u_y\\ u_z\end{array}\right)(x,y,z)\left(\begin{array}{c}u_x\\ u_y\\ u_z\end{array}\right)(x+\pi ,y+\pi /\sqrt{3},z).`$ (38)   
When applied to the flow these transformations induce changes in the variables $`y_i`$ (typically exchanges or sign changes), but the equations of motion are invariant under these transformations. Thus, if a certain flow has this symmetry, it leads to constraints on the variables $`y_i`$, and if it does not have this symmetry immediately a new flowfield can be obtained by applying this symmetry transformation. We do not attempt to analyze the full symmetry structure here and confine our discussion to two illustrative examples which are relevant for the stationary states discussed below. Demanding invariance of the flow field to the reflection symmetry $`T_1`$ leads to the following constraints on the variables $`y_i`$:

$`y_1`$ $`=`$ $`y_3=y_5=y_8=y_{15}=0`$ (39)   
$`y_6`$ $`=`$ $`y_7y_{10}=y_{12}`$ (40)   
$`y_{11}`$ $`=`$ $`y_{13}y_{16}=y_{18}y_{17}=y_{19}.`$ (41)   
The non vanishing components, $`y_2`$, $`y_4`$, $`y_6=y_7`$, $`y_9`$, $`y_{10}=y_{12}`$, $`y_{11}=y_{13}`$, $`y_{14}`$, $`y_{16}=y_{18}`$, $`y_{17}=y_{19}`$ thus define a 9 dimensional subspace.

For the combined symmetry $`T_1T_2`$ we find the constraints

$`y_1`$ $`=`$ $`y_3=y_5=y_8=y_{15}=0`$ (42)   
$`y_6`$ $`=`$ $`y_7y_{10}=y_{12}`$ (43)   
$`y_{11}`$ $`=`$ $`y_{13}y_{16}=y_{18}y_{17}=y_{19}`$ (44)   
and again a 9 dimensional subspace with non vanishing components $`y_2`$, $`y_4`$, $`y_6=y_7`$, $`y_9`$, $`y_{10}=y_{12}`$, $`y_{11}=y_{13}`$, $`y_{14}`$, $`y_{16}=y_{18}`$, $`y_{17}=y_{19}`$. The dimensions of the invariant spaces vary from a minimum of 6 (for each a $`T_1T_3`$ and $`T_1T_2T_3`$ invariance) and a maximum of 10 (for $`T_2T_3`$-invariance).

As mentioned, one can classify flows according to their symmetries. The most asymmetric flows are eightfold degenerate as the application of the eight combinations of the symmetries give eight distinct flows. The laminar flow profile is invariant under all the linear transformations and is the only member of the class with highest symmetry. The other stationary states discussed below fall into equivalence classes with eight members or four members if they are invariant under $`T_1`$ or $`T_1T_2`$.

## III Dynamics of perturbations

A stability analysis shows that the laminar flow profile is linearly stable for all Reynolds numbers. The matrix of the linearization is non-normal with a block structure along the diagonal. To bring this structure out more clearly, it is best to order the equations in the sequence 1, 2, 3, 4, 5, 7, 15, 8, 9, 14, 13, 19, 12, 18, 6, 11, 17, 10, 16. The matrix of the linearization then is upper diagonal, with a clear block structure: there are 10 eigenvalues isolated on the diagonal, three $`2\times 2`$ blocks and one $`3\times 3`$ block as well as several couplings between them in the upper right corner. While some eigenvalues can be complex, all of them have negative real part as shown in Fig. 3. For vanishing viscosity, the eigenvalues become zero or purely imaginary.

Large amplitude perturbations, however, need not decay. Already in the linear regime the non-orthogonality of the eigenvectors can give rise to intermediate amplifications into a regime where the nonlinear terms become important . In a related study on plane Couette flow  we used the lifetime of perturbations to get information on the dynamics in a high-dimensional phase space. As in that case, the amplitude of the velocity field in the $`z`$-direction indicates the survival strength of a perturbation. Linearizing the equations of motion around the base flow $`𝐮_0`$ gives for the perturbation $`𝐮^{}`$ the equation

$$_t𝐮^{}=(𝐮_0)𝐮^{}(𝐮^{})𝐮_0p^{}+\nu \mathrm{\Delta }𝐮^{}.$$
(45)   
The second term on the right hand side describes the energy source for the perturbation, and depends, because of $`𝐮_0=u_0(z)𝐞_y`$ and thus

$$(𝐮^{})𝐮_0=u_z^{}_zu_0(z)𝐞_y$$
(46)   
in an essential way on the $`z`$-components of the perturbation. Thus, if the amplitudes $`y_8`$, $`y_{10}`$, $`y_{12}`$, $`y_{14}`$, $`y_{16}`$ and $`y_{18}`$ become too small, the decay of the perturbation cannot be stopped any more. These modes account also for most of the off-diagonal block-couplings. A model for sustainable shear flow turbulence has to include some of these modes.

We chose a fixed initial flow field with a random selection of amplitudes $`y_1,\mathrm{},y_{19}`$, scaled it by an amplitude parameter $`A`$ and measured the lifetime as a function of $`A`$ and Reynolds number Re. Fig. 4 shows the time evolution of such a perturbation at $`Re=400`$ with one mode driven and for different amplitudes. For small $`A`$ there is an essentially exponential decay, whereas for larger amplitudes the perturbation swings up to large amplitude and shows no sign of a decline at all. The results for many amplitudes and Reynolds numbers are collected in Fig. 5 in a landscape plot. For small Reynolds number and/or small amplitude the lifetimes of perturbations are short, indicated by the light areas. For Reynolds numbers around 100 isolated black spots appear, indicating the occurence of lifetimes larger than the integration time (which increases with $`Re`$ so that $`t_{max}/Re=4\pi `$). The spottiness for $`Re`$ between about 100 and 1000 is due to rapid changes in lifetimes from pixel to pixel. For Re above 1000 the long lifetimes dominate. These results are in good agreement with what has been observed in plane Couette flow. Fig. 5b shows a similar plot for the case with two modes driven; it is qualitatively similar, but quantitatively shifted to higher Reynolds numbers.

In connection with the non-normality of the linearized eigenvalue problem it has been argued that the upper limit on the size of perturbations for which the non-linear terms in the dynamics can be neglected decreases algebraically like $`Re^\alpha `$. Different exponents have been proposed, ranging from 1 to 3 . It seems that for large $`Re`$ (where the model is actually less reliable because of the limited spatial resolution) the envelope of the long lived states in the fractal life time plot decays like $`Re^1`$.

The sensitive dependence of lifetimes on initial conditions and parameters is further highlighted in Fig. 6 and 7. The first shows the lifetime in the plane of the amplitudes $`y_{16}`$ and $`y_{17}`$ at Reynolds number $`Re=400`$ with all other components fixed. There is considerable structure on many scales. One notes ’valleys’ of short lifetimes between ’plateaus’ of longer lifetimes and a granular structure within both. The striations are reminiscent of features seen near fractal basin boundaries . Fig. 7 shows successive magnifications of lifetime versus amplitude plots at Re=200. Even after a magnification by $`10^7`$ there is no indication of a continuous and smooth variation of lifetime with amplitude.

## IV Stationary states

Motivated by the observation of new stationary structures in plane Couette flow for sufficiently high Reynolds number  we searched for non-trivial stationary solutions and studied their generation, evolution and symmetries.

We computed the stationary states with the help of a Monte Carlo algorithm. The initial conditions for the $`y_i`$’s were chosen randomly out of the interval $`[1/2,1/2]`$ and the Reynolds number was chosen randomly matrixwith an exponential bias for small $`Re`$ in the interval $`[10,10000]`$. With these initial conditions we entered a Newton algorithm. If the Newton algorithm converged, we followed the fixed point in Reynolds number as far as possible. We included about 200000 attempts in the Monte Carlo search.

The stationary states found for a single driven mode are collected in Fig. 8. No stationary states (besides the laminar profile) were found for Reynolds numbers below about 190. Between 190 and about 500 there are eight stationary states which divide into two groups of four symmetry related states each. With increasing Reynolds number more and more stationary states are found and they reach down to smaller and smaller amplitude. The envelope of all states reflects the $`Re^1`$ behaviour found for the borderline where nonlinearity becomes important. For two driven modes (Fig. 9) the situation is similar.

The appearance of the branches of the stationary states and in particular their coalescence near $`Re=190`$ suggests that the states are born out of a saddle-node bifurcation. And indeed, the eigenvalues as a function of $`Re`$ show two eigenvalues moving closer together and collapsing at zero for $`Re=190.41`$ (Fig. 10). However, these eigenvalues are not the leading ones, so that one set of states has three unstable eigenvalues, the other two unstable ones. It is thus a ‘saddle-node’ bifurcation into unstable states.

With increasing $`Re`$ more and more stationary states appear, partly through secondary bifurcations, partly through additional saddle-node bifurcations. Their number increases rapidly with Reynolds number (Fig. 11) and this increase goes in parallel with the increase in density of long lived states, Fig. 5. The detailed structure of the bifurcation diagram is rather complex and has not yet been fully explored. We note here that the various stationary states may be grouped according to their symmetries introduced in section II D and that we found only stationary states which belong to equivalent classes with four or eight members. The stationary states of the classes with four members are invariant under the transformation $`T_1`$ or $`T_1T_2`$. In addition, there are forward directed bifurcations generating two new branches with the same symmetry properties (eight or four member class) and inverse bifurcations of two branches belonging to eight member equivalent classes. We also found a backward directed bifurcation generating branches of an eight member class, which is born out of a four member class branch. The scenarios described above are marked in the bifurcation diagram Fig. 8.

## V Concluding remarks

The few degrees of freedom shear model introduced here lies halfway between the simplest models of non-normality and full simulations. Its dynamics has turned out to be surprisingly rich. There are a multitude of bifurcations introducing new stationary states besides the laminar profile, there are secondary bifurcations, and the distribution of life times shows fractal structures on amazingly small scales. It seems that as one goes from the low-dimensional models via the present one to full simulations one notes not only an increase in numerical complexity but also the appearence of qualitatively new features .

The simplest models with very few degrees of freedom focus on the non-normality of the linearized Navier-Stokes problem and emphasize the amplification of small perturbations. If the non-linearity is included a transition to another kind of dynamics, sometimes as simple as relaxation to a stationary point, is found .

Next in complexity are models like the one presented here that share with the few degree of freedom models the amplification and the transition but the additional degrees of freedom allow for chaos. When nonlinearities become important the dynamics does not settle to a fixed point or a limit cycle but continues irregularly for an essentially unpredictable time. The time is unpredictable because of the fractal life time distribution which seems to persist down to amazingly small scale: tiny variations in Reynolds number or amplitudes of the perturbation can cause major variations in life times. This fractal behaviour is the new quality introduced by the additional degrees of freedom. Indications for this behaviour are seen in the experiments by Mullin on pipe flow . It is interesting to ask just how few degrees of freedom are necessary to obtain this behaviour. Reducing our model to the $`T_1`$ subspace gives one with just 9 degrees of freedom (comparable in number and flow behaviour to the ones of Waleffe ) that still shows this fractal life time distributions. Further reduction, as in the four mode model of , seems to eliminate them.

The full, spatially extended shear flows share essential features with the model but add new problems. Spatially resolved simulations of the present model  as well as plane Couette flow with rigid-rigid boundary conditions  show the occurence of additional stationary states at sufficiently high Reynolds number that are unstable. A novel and as yet unexplained feature in spatially extended plane Couette flow, which we believe to be connected to the high dimensionality of phase space, is the difference between Reynolds numbers where the first stationary states are born (about 125 in units of half the gap width and half the velocity difference) and the ones where experiments begin to see long lived states (about 300–350) .

The fractal life time distributions have obvious similiarities to chaotic scattering . Drawing on this analogy one would like to identify permanent structures in phase space away from the laminar profile that could sustain turbulence. This has partly been achieved by the search for stationary states. Many have been found but irritatingly only for Reynolds numbers above about 190 while long lived states seem to appear much earlier. The solution to this puzzle must be periodic states and indeed we have found a few periodic states in a symmetry reduced model at lower Reynolds numbers, close to the occurence of the first long lived states. This suggests that the dynamical system picture that long lived states have to be connected to persistent structures in phase space is tenable.

There are several features of the model that can be studied further. In particular, quantitative characterizations of the fractal life time distribution, visualizations of the flow field, a detailed analysis of the primary and secondary bifurcation, an investigation of the dependence on the aspect ratio of the periodicity cell are required and look promising. We expect the lessons to be learned from this simple model to be useful in understanding the dynamics of full plane Couette and other shear flows. Work along these directions continues.