A

# Transient growth in linearly stable Taylor-Couette flows

## Abstract

Non-normal transient growth of disturbances is considered as an essential prerequisite for subcritical transition in shear flows, i.e. transition to turbulence despite linear stability of the laminar flow. In this work we present numerical and analytical computations of linear transient growth covering all linearly stable regimes of Taylor–Couette flow. Our numerical experiments reveal comparable energy amplifications in the different regimes. For high shear Reynolds numbers Re the optimal transient energy growth always follows a scaling, which allows for large amplifications even in regimes where the presence of turbulence remains debated. In co-rotating Rayleigh-stable flows the optimal perturbations become increasingly columnar in their structure, as the optimal axial wavenumber goes to zero. In this limit of axially invariant perturbations we show that linear stability and transient growth are independent of the cylinder rotation ratio and we derive a universal scaling of optimal energy growth using WentzelâKramersâBrillouin theory. Based on this, a semi-empirical formula for the estimation of linear transient growth valid in all regimes is obtained.

i
\checkfont

eurm10 \checkfontmsam10 Transient growth in Taylor–Couette flows]Transient growth in linearly stable Taylor–Couette flows S. Maretzke, B. Hof and M. Avila]S\lsI\lsM\lsO\lsN\nsM\lsA\lsR\lsE\lsT\lsZ\lsK\lsE1, B\lsJ\lsÖ\lsR\lsN\nsH\lsO\lsF and M\lsA\lsR\lsC\nsA\lsV\lsI\lsL\lsA

nstability, transition to turbulence

## 1 Introduction

The flow of viscous fluid between two coaxial independently and uniformly rotating cylinders, Taylor–Couette flow, is a paradigmatic system to study the stability and dynamics of rotating shear flows. For simplicity, we assume here that the system is infinite in the axial direction so that the annular geometry is uniquely determined by the dimensionless radius ratio of the inner and outer cylinders. A sketch of the Taylor–Couette system is shown in figure 1a.

The laminar Couette flow is determined by the inner and outer Reynolds numbers and , which are proportional to the rotation frequencies of the cylinders, and , respectively (see figure 1a). It is well known that the stability of Couette flow not only depends on the magnitudes of and , but also changes qualitatively with their ratio. In particular, Couette flow is stable to infinitesimal inviscid disturbances if and only if the fluid particles’ angular momentum increases in the radial direction. This result is known as Rayleigh’s criterion (Rayleigh, 1917). Consequently, inviscid instabilities solely depend on the ratio . Throughout this work, the term Rayleigh (un)stable is used to refer to the stability of Couette flow to inviscid disturbances. For viscous disturbances, there is a complex interplay of shear and centrifugal mechanisms determining the stability of the laminar flow (solid curve in figure 1b). Herein we use the expression linearly (un-) stable to refer to the viscous case.

In an attempt to separate the different effects that govern viscous stability, we adopt the parametrization introduced by Dubrulle et al. (2005), using shear Reynolds number Re and rotation number to parametrize the - plane (see figure 1b). As the name suggests, is a measure for the (absolute) shear in the flow, whereas depends solely on the ratio .

In the Rayleigh-stable regimes I and II in figure 1b the flow is also linearly stable. The remaining regimes III and IV are Rayleigh-unstable. Here viscosity has a stabilizing effect and the laminar flow first develops linear instabilities at finite non-zero Reynolds numbers. These already appear at moderate , except when approaching the boundaries of regimes I and II (Taylor, 1923). Indeed, the viscous linear stability boundary in figure 1b, determined from our numerical eigenvalue computations, shows that regime IV contains a relatively large range of moderate Reynolds numbers just above the axis, in which the viscous laminar flow is linearly stable. Note that defines the boundary to regime I, as the sign of is irrelevant here. In contrast, the region of linear stablity in III turns out to be negligibly small. Therefore, this regime is not studied in the present work, as the focus lies on linearly stable flows.

We subdivide the Rayleigh-stable regime according to the angular velocity profile of the base flow: The quasi-Keplerian regime II is characterized by , i.e. radially decreasing angular velocity, whereas regime I is defined by a positive gradient . In figure 1b these domains are separated by the solid-body line given by . Because of the absence of shear, for these configurations is constant, corresponding to a solid-body rotation flow profile. The transition from regime II to III defines the Rayleigh line where Rayleigh’s stability criterion ceases to be fulfilled and a centrifugal (linear) instability of the laminar flow emerges. In experiments, this results in the formation of a new stationary flow, characterized by the famous toroidal Taylor vortices (Taylor, 1923). Similar instabilities and associated patterns occur in the counter-rotating regime IV above the linear stability boundary plotted in figure 1b. For moderate Re very good agreement between this theoretical curve and experimentally observed instabilities has been achieved.

However, similarly to plane Couette and Poiseuille flow (cf. Romanov, 1973; Davey, 1973; Drazin & Reid, 1981), certain Taylor–Couette flows may undergo subcritical transition to turbulence in the absence of unstable eigenvalues. This phenomenon has been observed both by Coles (1965) in the Rayleigh-unstable counter-rotating regime IV as well as by Wendt (1933) and Taylor (1936) for a stationary inner cylinder (i.e. at the lower boundary of the Rayleigh-stable regime I: see figure 1b). Recent studies by Borrero-Echeverry et al. (2010) have confirmed the rapid lifetime increase of turbulent spots with the Reynolds number in the latter setting. Hence, we may infer the existence of subcritical turbulence within regime I in spite of the lack of experimental and numerical data for such flows.

On the other hand, the existence of turbulence remains controversial in the equally Rayleigh-stable quasi-Keplerian regime II (Yecko, 2004; Ji et al., 2006; Paoletti & Lathrop, 2011; Balbus, 2011). As the name suggests, these flows are of great importance in modelling astrophysical objects with Keplerian velocity profiles, such as accretion disks (for details, see Pringle, 1981). However, endcap effects render this regime difficult to explore experimentally. In fact, Avila (2012) has shown state-of-the-art Taylor–Couette apparatus to be possibly unsuited to adequately produce the respective flow fields at the required Reynolds numbers. Based on Re bounds derived from a variational formulation of the stability problem, Busse (2007) conjectured that turbulence cannot exist in the quasi-Keplerian regime. Yet, this result is predicated on the hypothesis that the extremizing vector fields are independent of the streamwise coordinate. To the best of our knowledge there is no general proof ruling out the existence of turbulence in the literature.

Whether linear or nonlinear, stability analysis boils down to the evolution of initial perturbations to a stationary state. For stationary flows, the development of the perturbation energy is given by the Reynolds-Orr equation, which is valid for both fully nonlinear and linearized dynamics (Schmid & Henningson, 2001). Remarkably, this implies that nonlinear instabilities may exist only if the linearized Navier-Stokes equations have solutions that grow in energy, i.e. transition requires linear growth.

At first glance, this theory seems contradictory to subcritical transition being a manifestation of nonlinear instability despite linear stability. However, the apparent paradox is resolved by the non-normality of the linearized Navier-Stokes operator, i.e. the non-orthogonality of its eigenmodes (Kato, 1995). This potentially allows for transient growth of infinitesimal perturbations (Boberg & Brosa, 1988; Trefethen et al., 1993), i.e. temporary energy growth even in the case of linear stability (as illustrated, for example, by Grossmann, 2000). As in other flow geometries the non-normality of the Taylor–Couette operator grows with the shear Reynolds number Re so that the maximum energy amplification, , may reach several orders of magnitude at sufficiently large Re (Reddy & Henningson, 1993). For instance, numerical simulations by Yecko (2004) of the rotating plane Couette geometry showed an asymptotic scaling of for one quasi- Keplerian flow configuration in the limit .

Hristova et al. (2002) and Meseguer (2002) were the first to study transient growth in the Taylor–Couette system. Both studies investigate counter-rotating flows. The former focuses on the growth behaviour of a single axisymmetric mode, whereas the latter computes optimal linear energy amplifications at the subcritical stability boundary measured by Coles (1965). Most prominently, Meseguer (2002) partly observes a strong correlation and finds a sharp threshold value for relaminarization in the experiments. These results reinforce the potential significance of non-normal growth in subcritical transition.

This article is concerned with transient growth in all regimes of linearly stable Taylor–Couette flows, identifying universal properties, especially in the limit of high Reynolds numbers. After briefly presenting the governing equations of the Taylor–Couette problem and our numerical formulation in 2 and 3, we discuss some tests of the method and numerical issues of transient growth computations in 4. In 5 the main numerical results for the asymptotic scaling of optimal transient growth and the corresponding optimal perturbations are presented. Furthermore, a semi-empirical formula for the estimation of by Re and the cylinder radius ratio is obtained. The latter is revealed to be universal by the analytical results for axially independent perturbations derived in 6. For such disturbances we further verify the characteristic scaling via a WentzelâKramersâBrillouin (WKB) approximation to the linearized evolution equations in 7. In the final section 8 we discuss our results and draw some conclusions concerning subcritical instability.

## 2 The linearized Taylor–Couette problem

### 2.1 Principal equations

We consider an incompressible Newtonian fluid with kinematic viscosity confined between two coaxial independently rotating cylinders with radii that are infinite in the axial direction. The annular geometry and its governing parameters are visualized in figure 1a. Non-dimensionalized with the gap width as length scale, viscous time and the pressure scale , the system is governed by the dimensionless incompressible Navier–Stokes equations and continuity equation

 ∂tv=−(v⋅∇)v−∇~p+Δv (1a) ∇⋅v=0 (1b)

where is the reduced pressure and the velocity field of the fluid.

The independent variables are the viscous time and the spatial vector parametrized in cylindrical coordinates (see figure 1a). The dimensionless geometry parameters are given by , and the radius ratio . Let and be the (constant) angular velocities of the inner and outer cylinder, respectively. Defining the inner and outer Reynolds numbers and the no-slip boundary condition at the inner and outer cylinder walls read

 v|r=ri={Re}ieφandv|r=ro={Re}oeφ (2)

where , and denote the orthonormal radial, azimuthal and axial unit vectors. The unusual appearance of the Reynolds number in the boundary conditions is due to the non-dimensionalization with the viscous timescale .

A well-known solution of the boundary value problem (1) and (2) is laminar Couette flow , given by

 vB=(Ar+Br)eφand~pB=12A2r2+2ABln(r)−B22r2 (3a) A:={Re}o−η{Re}i1+ηandB:=η({Re}i−η{Re}o)(1−η)(1−η2). (3b)

In order to investigate its stability, equations (1) are linearized about yielding the linearized Navier–Stokesequations for the evolution of an infinitesimally small perturbation :

 ∂t~u=−(vB⋅∇)~u−(~u⋅∇)vB−∇~q+Δ~u (4a) ∇⋅~u=0 (4b) ~u|r=ri=~u|r=ro=0 (4c)

By a Fourier ansatz in the azimuthal and axial coordinates , for , the evolution equation can be written as

 ∂tu=Lu−∇cq. (5)

Herein a subscript for an operator denotes the conjugate with , i.e. . The operator is given by (Meseguer, 2002)

 Lu=−(vB⋅∇)cu−(u⋅∇)vB+Δcu=:⎛⎜⎝LrrLrφ0LφrLφφ000Lzz⎞⎟⎠⎛⎜⎝uruφuz⎞⎟⎠ (6a) Lrr=Lφφ=Lzz−1r2 = D+D−n2+1r2−k2−inrvBφ Lrφ = 2rvBφ−2inr2 Lφr = 2inr2−D+vBφ (6b)

with the abbreviations and . The domain of admissible velocity fields in equation (5) is the twice continuously differentiable subspace

 V:={v∈H3∩C2((ri;ro)):v(ri)=v(ro)=0,∇c⋅v=0} (7)

of the Hilbert space . Here we define with the inner product

 ⟨⋅,⋅⟩:H×H→C;(q1,q2)↦∫roriq∗1q2rdr. (8)

where the superscript denotes the conjugate transpose of a scalar, vector or matrix. For simplicity, we likewise denote the canonical inner product in , , by . The induced norm squared is proportional to the total kinetic energy of a perturbation and is therefore denoted as the energy norm.

A modal ansatz in the time coordinate , i.e. and for yields the eigenvalue problem

 λuλ=Luλ−∇cqλ,(uλ,qλ)∈V×H. (9)

For the axisymmetric case , DiPrima & Habetler (1969) have shown the discreteness of the eigenvalues and completeness of the corresponding generalized eigenfunctions in . If we assume that this remains true for , then the laminar Couette flow (3) is linearly stable if and only if all eigenvalues of (9) have negative real parts.

### 2.2 The parameter space for transient growth

In addition to the experimental parameters , and , the evolution problem (5) depends on the azimuthal and axial wavenumbers and . Owing to the cylindrical symmetry of the Taylor–Couette geometry the parametric analysis may be confined to (for details, see Meseguer & Marques (2000)). The parameter determines the curvature of the system and thus the rotational influence. The limit corresponds to plane Couette flow as demonstrated with respect to transient growth by Hristova et al. (2002), whereas implies infinite curvature at the inner cylinder wall.

For reasons discussed in 1 we introduce the shear Reynolds number Re and the rotation number . Assuming and , the mapping is one-to-one so that the flow parameters and can be expressed via Re and :

 Re:=2|ηReo−Rei|1+ηandRΩ:=(1−η)(Rei+Reo)ηReo−Rei (10a) A=sgn(RΩ)Re2(RΩ+1)andB=−sgn(RΩ)ηRe2(1−η)2. (10b)

The parametrization of the different Taylor–Couette flow regimes in figure 1b by is summarized in table 1. As can be seen from (10b), only the parameter , which governs the solid-body rotation part of the base flow (see (3a)) depends on the rotation number . The shear term is independent of modulo sign, which will be essential for the results of 6 and 7. On the other hand, we have , so that the shear Reynolds number determines the overall magnitude of the flow.

Computing the commutator of the operator given by (6) with its adjoint , , reveals its non-normality, scaling with the shear Reynolds number. The eigenspaces are therefore non-orthogonal to one another (Kato, 1995), which potentially allows for significant transient growth at large Re. Detailed discussions of this mechanism can be found in Grossmann (2000) and Schmid & Henningson (2001, pp. 99-101).

As a consequence, initial perturbations may temporarily grow in energy before they ultimately decay – even if has only stable eigenvalues with . The maximum transient growth at time is given by . The evolution of may be written as a linear equation of the form (where, strictly speaking, due to the remaining pressure dependence in (5)). Thus can be expressed using the operator norm (Trefethen et al., 1993):

 G(t):=sup∥u(0)∥=1∥u(t)∥2=sup∥u(0)∥=1∥exp(~Lt)u(0)∥2=∥exp(~Lt)∥2 (11)

If denotes the energy norm, is equal to the greatest kinetic energy amplification that an initial perturbation can attain at time .

For a Taylor–Couette flow configuration given by the parameters , and , the optimal transient growth is defined by . A perturbation with is called optimal if for some . Note that is finite if and only if all eigenvalues of are stable.

## 3 Numerical formulation and implementation

### 3.1 The Galerkin method

The eigenvalue problem (9) is numerically solved using a Galerkin method. The implementation is similar to the Petrov-Galerkin method described by Meseguer & Marques (2000) and Meseguer et al. (2007), but based on Legendre rather than Chebyshev polynomials so that trial and projection basis are identical.

The basis choice is , where and are defined according to table 2 for different wavenumbers and . The functions and are given by

 hm(r):=r(1−x2)Lm(x)andgm(r):=r(1−x2)2Lm(x)forr∈[ri;ro]. (12)

where is the Legendre polynomial of degree and . Then every satisfies both the continuity condition and the boundary conditions by definition, since we have by construction

 hm(ri)=hm(ro)=gm(ri)=gm(ro)=g′m(ri)=g′m(ro)=0. (13)

The problem is discretized by truncating at the polynomial resolution , i.e. defining , and expanding possible solutions to the eigenvalue problem (9) in terms of , . Plugging this ansatz into equation (9) and projecting on some yields

 λ∑m

The pressure terms vanish due to the boundary and divergence conditions.

Thus, equations (14) for all , , can be written in the form of a generalized eigenvalue problem

 Missing or unrecognized delimiter for \right (15)

for the coefficient vector where and are -matrices being Hermitian positive definite (Meseguer & Marques, 2000).

### 3.2 Computation of transient growth

Now let be the eigenfunctions corresponding to the eigenvalues and eigen(-coefficient-)vectors solving the generalized eigenvalue problem (15). Consider some perturbation expanded in , i.e. where denotes the time-dependent coefficient vector. Since the are (approximate) solutions to the eigenvalue problem (9), it follows that

 b(t)=exp(diag(λ)t)b(0) (16)

where denotes the diagonal matrix constructed from and exp is the matrix exponential. Thus the evolution of the perturbations kinetic energy reads

 ∥u∥2=⟨u,u⟩=2N∑i,j=1b∗ibj⟨qi,qj⟩=b∗\mathsfbiMb=∥\mathsfbiFb∥22=∥\mathsfbiFexp(diag(λ)t)b(0)∥22 (17)

Here is the Hermitian positive definite Gramian matrix , a Cholesky decomposition and denotes the standard 2-norm on . Hence, the maximum transient growth at time is given by (see Meseguer, 2002)

 G(t) = sup∥u(0)∥=1∥u(t)∥2=sup∥\mathsfbiFb∥2=1∥\mathsfbiFexp(diag(λ)t)b∥22 (18) \lx@stackrelv=\mathsfbiFb= sup∥v∥2=1∥\mathsfbiFexp(diag% (λ)t)\mathsfbiF−1v∥22=∥\mathsfbiFexp(diag(λ)t)\mathsfbiF−1∥22.

So is equal to the squared maximum singular value of . Moreover, if denotes the corresponding right-singular vector, is the initial -coefficient vector of a perturbation that attains optimal transient growth at time . By means of singular value decomposition, we thus obtain both maximum transient growth and corresponding perturbations in the finite-dimensional subspace spanned by . This yields a lower bound to the maximum attainable by arbitrary initial conditions in . As discussed in 4.2 we find convergence of this estimate to the total maximum.

### 3.3 Outline of the code

By definition of only integrals over polynomial functions have to be evaluated in order to calculate and . Hence, these are computed exactly using Gauss–Legendre quadrature with Gauss-Lobatto collocation points of degree , where (see Canuto et al., 2006, pp. 69 ff.). Moreover, the derivatives in the operator are implemented by means of the corresponding differentiation matrices given in Canuto et al. (2006, p. 76).

The code used in this work is based on the scientific computing package Scipy for the interactive language Python. The linear algebra algorithms are provided by the package Scipy.Linalg based on the standard ATLAS, LAPACK and BLAS implementations.

The optimization of in the time coordinate and in the continuous wavenumber are performed via the Scipy.Optimize implementation of Brent’s method (for details see Press et al., 2007, sec. 9.3). With respect to the discrete wavenumber , is optimized by brute force. If the optimal transient growth is found at the upper boundary of the considered domains, i.e. for , or the respective intervals are enlarged in subsequent steps until a local maximum is located in their interior.

## 4 Numerical issues

In this section, the performance of the numerical implementation presented in section 3 is tested by comparison to results in the literature. Furthermore, eigenvalue and transient growth convergence are studied for test cases in order to justify the choice of polynomial resolution used to obtain the numerical results in 5. We find that the optimal transient growth may converge without the Y-shaped spectrum being properly resolved. This observation suggests a minor significance of the spectrum in transient growth computations, disagreeing with the conclusions drawn by Reddy & Henningson (1993) for channel flows.

### 4.1 Eigenvalue decomposition

Our discretization of the eigenvalue problem (9) has been tested against the results on eigenvalue-critical Reynolds numbers presented in Krueger et al. (1966, table 2) as well as by replication of the plotted spectra given by Gebhardt & Grossmann (1993, fig. 3a-d). Agreement within the respective accuracies has been found. Additionally, we have compared our Galerkin method to the Petrov-Galerkin scheme of Meseguer et al. (2007). No significant deviations are found between the converged spectra.

For these methods we study the convergence of the approximated least stable eigenvalue as the number of Legendre or Chebyshev polynomials is increased. In figure 2 the relative errors compared to (converged) reference values are plotted against . The test parameters are , , and at shear Reynolds numbers for figure a and for figure b.

The plots in figure 2 show plateaus of non-convergence for low , which are due to the difficulty of identifying the respective eigenvalue in a non-converged spectrum. For moderate resolutions ( in figure a and in b) spectral accuracy, i.e. exponential convergence rates, is observed for both methods. Notably however, the convergence turns out to be significantly quicker for the Legendre-polynomial-based Galerkin method presented in this work: spectral accuracy is attained using significantly fewer polynomials and the limiting machine precision is reached already for () and () compared to and , respectively, in the case of the Petrov-Galerkin scheme (see figure 2).

The required resolution for convergence grows with the shear Reynolds number Re and – much more significantly – as soon as subsequent, more stable eigenvalues are considered. In fact, it turns out to be numerically impossible to resolve significant parts of the eigenvalue spectrum for . This also affects the computation of transient growth discussed in the next subsection.

### 4.2 Computation of Transient Growth

In table 3 our results concerning the optimal transient growth for and the corresponding optimal wavenumbers are compared to the numerical data of Meseguer (2002, table 1). The values of and differ by less than .

The convergence of the maximum transient growth shows remarkable characteristics which partly contradict the significance of the linearized Navier-Stokes operator’s spectrum for such computations claimed, for example, by Reddy & Henningson (1993).

These features are discussed with reference to the example displayed in figure 3: for three different resolutions (corresponding to figures a, b and c), the eigenvalues (top), the evolution of the maximum transient growth (middle) and the moduli of the components , and of the corresponding optimal perturbation are plotted for comparison. The example parameters are , , , and .

A few aspects are noteworthy. Around its maximum is already surprisingly well approximated by only Legendre polynomials, whereas the optimal perturbation is far from its actual shape (see figure a). For (figure b) the curve is converged within an error while its maximum is even approximated up to . Likewise, the optimal perturbation is practically converged. At the same time the characteristic Y-like structure of the eigenvalue spectrum (cf. Gebhardt & Grossmann, 1993) is by no means well resolved for not to mention (top). In fact, it takes as many as polynomials for convergence of the two meeting branches (see figure c). However, this does not seem to affect the transient growth quantities – even though the converged spectrum in figure c (top) is even much more stable as a whole than its approximation for (figure b).

In contrast to these observations Reddy & Henningson (1993) stress the significance of the two eigenvalue branches and especially their meeting point for transient growth in channel flows. As for Taylor–Couette flow, this is only confirmed if the Y structure is resolved in the first place. This turns out not to be necessary, which is a lucky circumstance in two respects. On the one hand, the two branches consist of discrete eigenvalues for , rendering their convergence numerically infeasible for . This agglomeration of eigenvalues can be explained by the dominance of the convective multiplicative terms in the linearized Taylor–Couette operator over the viscous differential contributions in the limit (see equations (6a) and (6b)). The asymptotic degeneracy into a pure multiplication operator, i.e. a mapping with a continuous, matrix-valued function , corresponds to a transition from discrete eigenvalues to a continuous spectrum.

On the other hand the standard Cholesky decomposition of the matrix (see 3.2) tends to fail at large Re if the eigenvalue spectrum is over-resolved. In the example shown in figure 3 this happens for – just as the crucial meeting point is resolved. Accordingly, one might expect to miss a sudden jump in the maximum transient growth if the method breaks down precisely at this point. Note, however, that no such discontinuity is observed in those cases where the intersection can still be resolved, i.e. for smaller Re.

We may thus conclude that the transient growth of the linearized Taylor–Couette operator is already converged while its approximated spectrum is still far from its natural shape. Startling at first glance, this is yet another manifestation of transient growth’s non-modal nature: the non-eigendirections are those of significance.

Nevertheless, numerical artifacts in the form of spurious unstable eigenvalues have to be avoided by choosing sufficiently high resolutions . However, must not be too large either in order to keep the Cholesky decomposition stable (although preconditioning or more stable algorithms such as the one presented by Ogita & Oishi (2012) might be another alternative). For a given set of parameters , Re, it turns out that resolving the transient growth peak for optimal wavenumbers , tends to require the highest resolutions. Moreover, the necessary are mostly independent of and at least of the same magnitude for different . Here greater curvature, i.e. , results in slower convergence. Consequently, for practical computations, suitable resolutions are determined for different ranges of Re by the convergence of (computationally challenging) test cases, more precisely less than deviation in the optimal transient growth for and .

For greater lower resolutions may be sufficient and greater Reynolds numbers than might be resolvable. However, universal convergence for any parameters , , and within about may be assumed if is chosen according to the resulting canonical resolutions given in table 4. They are found to approximately follow a power law of the form with and .

Starting from these, is temporarily reduced in subsequent steps whenever the Cholesky decomposition fails and temporarily increased if unstable eigenvalues occur in order to identify possible numerical artifacts. In the case of converged unstable eigenvalues the computation of the matrix and thus of the transient growth is confined to the stable eigenmodes in in agreement with the analysis of Meseguer (2002).

These computation guidelines have been applied to obtain the numerical results presented in 5.

## 5 Numerical results

In this section the numerical results concerning stability and transient growth in Taylor–Couette flows are presented.

### 5.1 Optimal transient growth in various regimes

According to the numerical strategy discussed in 3 and 4.2, the optimal transient growth is computed for logarithmically equidistant shear Reynolds numbers and .

By studying test cases we find that with and is a suitable choice to determine the transient growth maximum in time for all considered parameter regimes. Optimization in the wavenumbers is carried out by default in the range and . Additionally, as discussed in 3.3 the ranges for , and are enlarged whenever the optimization terminates near one of the upper bounds.

By the choice of rotation numbers the linearly stable regimes I and II are parametrized considering (I) and (II) (see table 1). Furthermore, transient growth is studied in the counter-rotating regime IV near the line in the - left quadrant by choosing . For a global overview the results for are presented. The lines in the - plane defined by this choice for are visualized in figure 4 for orientation, along with the numerically computed (viscous) linear stability boundary. Figure 5 shows the optimized transient growth and the corresponding optimal axial wavenumber , respectively, against Re for the considered parameter sets.

The most prominent feature in the double-logarithmic plots of figure a are the nearly identical asymptotic slopes of the lines in the linearly stable regimes for showing a characteristic power law with (compare dashed line in figure a). Notably, even in the Rayleigh-unstable counter-rotating regime IV (circles in figure a), seems to approach this scaling as long as the computation is not destabilized by dominant linear instability. In fact, for constant Re the energy amplifications in the different regimes differ only by factors and not – as possibly expected – by orders of magnitude. Within the linearly stable regimes I and II these deviations are most distinct in the vicinity of the Rayleigh line and the boundary to regime IV where larger amplifications occur.

Hence, the numerical results suggest that optimal transient growth in linearly stable Taylor–Couette flows roughly follows a common scaling for . Note that this scaling result is in perfect agreement with those by Yecko (2004) obtained for Keplerian flows at fixed in rotating plane Couette geometry.

### 5.2 Optimal axial wavenumber

Beyond the magnitude of transient growth studied in 5.1 the spatial structure of the optimal perturbations is of great interest. The latter is determined by the optimal axial and azimuthal wavenumbers and which attain the optimal transient growth shown in figure a. The are plotted in figure b with logarithmic horizontal axes. Note the discontinuities of the curves whenever the (discrete) optimal azimuthal wavenumbers changes.

The plots reveal a characteristic quasi-two-dimensional, columnar structure of the optimal perturbations in regime II (also observed by Yecko (2004)): for the optimal transient growth is consistently attained by axially independent perturbations, i.e. (compare and in figure b). The transition to typically occurs already for Reynolds numbers as small as . Only near the Rayleigh line – that is, for – is a sharp divergence of for is observed. Here holds up to the greatest studied shear Reynolds numbers .

While is only obtained in the quasi-Keplerian regime (II), in regime I (represented by ) seems to decay (slowly) to zero for . At least weak axial dependence is observed for in these flows. Once again, the asymptotic decay is most distinct near the solid-body line and is lost near the transition to counter-rotation at . Here an almost constant optimal wavenumber is observed.

For further illustration figures 6 and 7 show contour plots of in the regimes II and I, respectively. The boundary lines have been computed by a bisecting algorithm with relative accuracy . The extent of the shaded regions in figure 6 emphasizes the dominance of axially independent, columnar perturbations for quasi-Keplerian flows. \sidecaptionvposfigurec

In the counter-rotating regime (IV) we observe a growing with Re. This difference might be explained by emerging linear instabilities which first appear for in this regime and thus render fully three-dimensional perturbations less dissipative.

The behaviour of the optimal azimuthal wavenumber is not discussed in detail here. Notably however, axisymmetric perturbations (corresponding to ) never attain significant energy growth up to high Reynolds numbers except for a small neighbourhood of where the dominant Taylor-vortex-related instability of regime III emerges. On the other hand, usually transient growth of the same order is attained for different . Numerical results indeed suggest that for sufficiently large shear Reynolds numbers, depends more on the geometrical parameter rather than on Re or which parametrize the base flow. In general, an azimuthal wavenumber seems to be optimal if the associated wavelength is , i.e. twice the gap width, leading to vortices, that are of about the same radial and streamwise dimension (see e.g. figure a, centre right).

In contrast, the dominant axial wavenumbers in regime I correspond to wavelengths of rather than gap widths. The axial dependence of the optimal perturbations is thus indeed weak compared to azimuthal (and radial) variations. For comparison, recall that one observes axial symmetry and order-one axial wavelengths for the usual Taylor vortices corresponding to and . Moreover, we observe that, the stronger the rotational influence on the fluid’s stability expressed by smaller and/or larger , the smaller are the attained for (figures b). The observed columnwise preference of the optimal perturbations is thus in good agreement with the Taylor–Proudman theorem, stating that a rapidly rotating inviscid fluid is (preferably) uniform along its rotational axis. On the other hand, this preference does not seem to be manifested in the dominant least stable eigenmodes observed in quasi-Keplerian flows: numerical optimization of the principal eigenvalue’s real part over and in the test cases , and (data not plotted) indeed shows significantly non-columnar modes with to be least dissipative for . The principal zero mode , on the other hand, is found to decay about one order of magnitude more slowly than the optimal non-axisymmetric ones in the considered parameter range. Note, furthermore, that eigenvalues corresponding to perturbations with predominantly streamwise (i.e. azimuthal) or spanwise (axial) flow, respectively, alternate along the real axis in the least stable parts of all studied spectra, where the spanwise modes even turn out to be slightly less stable. The study thus demonstrates that the structure of optimal non-modal perturbations may be entirely different from that of the dominant eigenmodes.

Changing does not seem to have any further qualitative effects on transient growth according to the results in figure 5, as long as none of the limits is considered. A further study of this parameter is therefore omitted in the following.

### 5.3 Evolution of optimal perturbations

In the sequel, three different optimal perturbations , and are considered at a constant shear Reynolds number and . The rotation numbers are given by , and corresponding to regimes II, I and IV. The optimal wavenumbers are given by , and and . The time evolution of these modes is computed by eigenmode decomposition at a polynomial resolution .

In figure 8 the resulting real parts of , and are shown at a sequence of snapshots throughout the transient growth evolution. The flow fields are plotted in radial–azimuthal projection (top) and radial-axial projection (bottom) with on the horizontal axis except for where the latter is omitted due to the axial independence. The radial-axial plots have been rescaled so that exactly one axial wavelength is displayed. Arrow lengths scale with the absolute flow velocities although different scalings are applied in figures a, b and c. The colour map from yellow to red marks regions with relatively low or high energy densities in the current fields.

The perturbations’ total kinetic energy evolution in relation to the transient growth maxima are plotted in figure 9 with time scale renormalized by . The considered in figure 8 are identified by markers. \sidecaptionvposfigurec

The radial–azimuthal projections in figure 8 reveal essentially similar transient growth mechanisms of the considered modes: the optimal initial perturbations have a spiral-like structure of streamwise elongated vortices. Recalling the different angular velocities and of the driving inner and outer cylinders(i.e. for , for and in the counter-rotating case , respectively), we find that the initial spiral orientations are always misfit to the base flow. This “misfit” character is a manifestation of the perturbations’ non-modal nature and thus typical of transient growth as emphasized by Grossmann (2000). The spiral velocity fields are tilted by the base flow and thereby gain energy (compare figures 8 centre-left and figure 9). As for the axially independent perturbation in a the energy maximum then occurs exactly at the turning point of the spiral orientation whereas in cases 2 and 3 it is attained shortly after this point (centre-right in figure 8). Subsequently, the perturbation is further deformed into a “fit” flow direction, i.e. an eigendirection, and meanwhile decays.

This shear-induced detilting dynamics of perturbations, with initial vorticity leaning against the background shear profile, essentially represents a Taylor–Couette analogue of the so-called Orr mechanism. The latter has been identified, e.g. in the early two-dimensional studies of optimal transient growth by Farrell (1988), as an important ingredient of linear non-modal growth in plane channel flows, providing a potential explanation for the emergence of finite-amplitude disturbances required for nonlinear instabilities. Notably, in the cases studied here, this mechanism leads to transient spiral structures that resemble those of the linearly unstable, axially independent eigenmodes reported by Gallet et al. (2010) – compare our figure a center-left with 4(d) in Gallet et al. (2010). The latter arise in the case of an additionally imposed radial inflow through the outer cylinder, which seemingly stabilizes the misfit tilt of the vortices, rendering the transient energy growth, observed in the present work, sustained.

Especially for the columnar axially independent perturbation , the energy growth and decay is rather sudden, leading to the sharp peak depicted in figure 9. This phenomenon is similar to the dynamics observed in plane channel flows and is possibly due to the rapid flow near the inner cylinder wall (see figure a, centre-right), which leads to high dissipation around the energy maximum. On the other hand, the optimal perturbations and in the regimes I and IV seem to be stabilized in this respect by their axial dependence, leading to and larger growth than that attained for and slower decay in figure 9. This interpretation is supported by the fact that, in spite of the small wavenumber