Transient growth in linearly stable Taylor-Couette flows

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 2/3-scaling with Re, 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 cylinders' rotation-ratio and we derive a universal 2/3-scaling of optimal energy growth with Re using WKB-theory. Based on this, a semi-empirical formula for the estimation of linear transient growth valid in all regimes is obtained.


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 1(a).
The laminar Couette flow is determined by the inner and outer Reynolds numbers Re i and Re o , which are proportional to the rotation frequencies of the cylinders, Ω i and Ω o , respectively (see figure 1a). It is well known that the stability of Couette flow not only depends on the magnitudes of Re i and Re o , but also changes qualitatively with their ratio. In particular, Couette flow is stable to infinitesimal † Email address for correspondence: simon.maretzke@googlemail.com 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 Re i /Re o . 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 R Ω to parametrize the Re i -Re o plane (see figure 1b). As the name suggests, Re ∼ Ω o − Ω i is a measure for the (absolute) shear in the flow, whereas R Ω depends solely on the ratio Ω o /Ω i . In the Rayleigh-stable regimes I and II in figure 1(b), 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 Re = O(10 2 -10 3 ), except when approaching the boundaries of regimes I and II (Taylor 1923). Indeed, the viscous linear stability boundary in figure 1(b), determined from our numerical eigenvalue computations, shows that regime IV contains a relatively large range of moderate Reynolds numbers (Re o , Re i ) just above the Re o axis, in which the viscous laminar flow is linearly stable. Note that Re i = 0 defines the boundary to regime I, as the sign of Re o 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 Ω B of the base flow. The quasi-Keplerian regime II is characterized by ∂ r Ω B < 0, i.e. radially decreasing angular velocity, whereas regime I is defined by a positive gradient, ∂ r Ω B > 0. In figure 1(b) these domains are separated by the solid-body line given by Re i = ηRe o . Because of the absence of shear, for these configurations Ω B 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 1(b). 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. Davey 1973;Romanov 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 & Schatz (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;Balbus 2011;Paoletti & Lathrop 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, G max , 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 G max ∼ Re 2/3 for one quasi-Keplerian flow configuration in the limit Re → ∞. 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 Re T (R Ω ) measured by Coles (1965). Most prominently, Meseguer (2002) partly observes a strong correlation and finds a sharp threshold value G max,T = 71.58 ± 0.16 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 G max ∼ Re α of optimal transient growth and the corresponding optimal perturbations are presented. Furthermore, a semi-empirical formula for the estimation of G max 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 G max ∼ Re 2/3 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.

Principal equations
We consider an incompressible Newtonian fluid with kinematic viscosity ν confined between two coaxial independently rotating cylinders with radii r i < r o that are infinite in the axial direction. The annular geometry and its governing parameters are visualized in figure 1(a). Non-dimensionalized with the gap width d := r o − r i as length scale, viscous time ν −1 d 2 and the pressure scale ν −2 d 2 , the system is governed by the dimensionless incompressible Navier-Stokes equations and continuity equation wherep is the reduced pressure and v is the velocity field of the fluid. The independent variables are the viscous time t and the spatial vector x parametrized in cylindrical coordinates x =: (r, ϕ, z) T (see figure 1a). The dimensionless geometry parameters are given by r i := r i d −1 , r o := r o d −1 and the radius ratio η := r i r −1 o . Let Ω i and Ω o be the (constant) angular velocities of the inner and outer cylinder, respectively. Defining the inner and outer Reynolds numbers Re i := (d/ν)r i Ω i and where e r =: (1, 0, 0) T , e ϕ =: (0, 1, 0) T and e z =: (0, 0, 1) T 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 time scale d 2 /ν.
In order to investigate its stability, equations (2.1) are linearized about (v B ,p B ), yielding the linearized Navier-Stokes equations for the evolution of an infinitesimally small perturbation (ũ,q): (2.4c) By a Fourier ansatz in the azimuthal and axial coordinatesũ(r, ϕ, z) := u(r)e i(nϕ+kz) , q(r, ϕ, z) := q(r)e i(nϕ+kz) for k ∈ R, n ∈ Z, the evolution equation can be written as Herein a subscript c for an operator T denotes the conjugate with e i(nϕ+kz) , i.e. T c := e −i(nϕ+kz) T e i(nϕ+kz) . The operator L is given by (Meseguer 2002) with the abbreviations D := ∂ r and D + := ∂ r + 1/r. The domain of admissible velocity fields u = (u r , u ϕ , u z ) T in (2.5) is the twice continuously differentiable subspace of the Hilbert space H 3 . Here we define H := L 2 ((r i ; r o )) with the inner product where the superscript * denotes the conjugate transpose of a scalar, vector or matrix. For simplicity, we likewise denote the canonical inner product in H 3 , (u 1 , u 2 ) → u 1,r , u 2,r + u 1,ϕ , u 2,ϕ + u 1,z , u 2,z , by ·, · . The induced norm squared u 2 := u, u is proportional to the total kinetic energy of a perturbation u and is therefore denoted as the energy norm.

Regime I
Regime II Regime III Regime IV Solid-body line Rayleigh line  figure 1(b); lines of constant R Ω are axes meeting in the origin.
A modal ansatz in the time coordinate t, i.e. u := u λ e λt and q := q λ e λt for λ ∈ C yields the eigenvalue problem (2.9) For the axisymmetric case n = 0, DiPrima & Habetler (1969) have shown the discreteness of the eigenvalues {λ} and completeness of the corresponding generalized eigenfunctions in V. If we assume that this remains true for n = 0, then the laminar Couette flow (2.3) is linearly stable if and only if all eigenvalues of (2.9) have negative real parts.
2.2. The parameter space for transient growth In addition to the experimental parameters Re i , Re o and η, the evolution problem (2.5) depends on the azimuthal and axial wavenumbers n and k. Owing to the cylindrical symmetry of the Taylor-Couette geometry, the parametric analysis may be confined to Re i , n, k 0 (for details, see Meseguer & Marques (2000)). The parameter η ∈ (0, 1) determines the curvature of the system and thus the rotational influence. The limit η → 1 corresponds to plane Couette flow as demonstrated with respect to transient growth by Hristova et al. (2002), whereas η → 0 implies infinite curvature at the inner cylinder wall.
For reasons discussed in §1, we introduce the shear Reynolds number Re and the rotation number R Ω . Assuming Re i 0 and Re i = ηRe o , the mapping (Re i , Re o ) → (Re, R Ω ) is one-to-one so that the flow parameters A and B can be expressed via Re and R Ω : The parametrization of the different Taylor-Couette flow regimes in figure 1(b) by R Ω is summarized in table 1. As can be seen from (2.10b), only the parameter A, which governs the solid-body rotation part ∝ Ar of the base flow (see (2.3a)) depends on the rotation number R Ω . The shear term ∝ B/r is independent of R Ω modulo sign, which will be essential for the results of § § 6 and 7. On the other hand, we have A, B ∝ Re, so that the shear Reynolds number determines the overall magnitude of the flow. Computing the commutator of the operator L given by (2.6) with its adjoint L * , [L * , L ] = O(Re 2 ), 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 Spectral basis functions for m ∈ N used for the discretization of the eigenvalue problem (2.9) via (3.1) and (3.3) according to Meseguer et al. (2007).
this mechanism can be found in Grossmann (2000) and Schmid & Henningson (2001, pp. 99-101). As a consequence, initial perturbations u(0) may temporarily grow in energy before they ultimately decay -even if L has only stable eigenvalues λ ∈ C with Re(λ) < 0. The maximum transient growth at time t 0 is given by G(t) := sup u(0) =1 u(t) 2 . The evolution of u may be written as a linear equation of the form ∂ t u =L u (where, strictly speaking,L = L due to the remaining pressure dependence in (2.5)). Thus G can be expressed using the operator norm (Trefethen et al. 1993): (2.11) If · denotes the energy norm, G(t) is equal to the greatest kinetic energy amplification that an initial perturbation u(0) ∈ V can attain at time t 0.
For a Taylor-Couette flow configuration given by the parameters Re i , Re o and η, the optimal transient growth is defined by G max := sup t,n,k G(t). A perturbation u with u(0) = 1 is called optimal if u(t) 2 = G max for some t 0. Note that G max is finite if and only if all eigenvalues ofL are stable.

The Galerkin method
The eigenvalue problem (2.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 U := {u j m } j=1,2 m∈N 0 , where u 1 m and u 2 m are defined according to table 2 for different wavenumbers n and k. The functions h m and g m are given by where L m is the Legendre polynomial of degree m and x := 2r − (1 + η)(1 − η) −1 . Then every u j m satisfies both the continuity condition ∇ c · u j m = 0 and the boundary conditions by definition, since we have by construction Transient growth in Taylor-Couette flows

261
The problem is discretized by truncating U at the polynomial resolution N ∈ N, i.e. defining U N := {u j m } j=1,2 m<N , and expanding possible solutions u λ to the eigenvalue problem (2.9) in terms of U N , u λ := m<N, j=1,2 a j m u j m . Plugging this ansatz into (2.9) and projecting on some u i l yields The pressure terms vanish due to the boundary and divergence conditions. Thus, equations (3.3) for all l < N, i = 1, 2, can be written in the form of a 2N × 2N generalized eigenvalue problem for the coefficient vector a := (a 1 0 , . . . , a 1 N−1 , a 2 0 , . . . , a 2 N−1 ) T , where G and H are 2N × 2N matrices, with G being Hermitian positive definite (Meseguer & Marques 2000).

Computation of transient growth
where diag(λ) denotes the diagonal matrix constructed from λ and exp is the matrix exponential. Thus the evolution of the perturbations kinetic energy reads Here M is the Hermitian positive definite Gramian matrix M := ( q i , q j ), M = F * F is a Cholesky decomposition and · 2 denotes the standard 2-norm on C 2N . Hence, the maximum transient growth at time t 0 is given by (see Meseguer 2002) So G(t) is equal to the squared maximum singular value σ 2 0 of F exp(diag(λt))F −1 . Moreover, if v 0 denotes the corresponding right-singular vector, F v 0 is the initial Q-coefficient vector of a perturbation that attains optimal transient growth at time t. By means of singular value decomposition, we thus obtain both maximum transient growth G(t) and corresponding perturbations in the finite-dimensional subspace spanned by Q. This yields a lower bound to the maximum attainable by arbitrary initial conditions in V. As discussed in § 4.2 we find convergence of this estimate to the total maximum.

3.3.
Outline of the code By definition of U N only integrals over polynomial functions have to be evaluated in order to calculate G and H. Hence, these are computed exactly using Gauss-Legendre quadrature with Gauss-Lobatto collocation points of degree M, where M N + 6 (see Canuto et al. 2006, pp. 69ff.). Moreover, the derivatives in the operator L 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 G in the time coordinate t ∈ [0; t cut ] and in the continuous wavenumber k ∈ [0; k cut ] 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 n ∈ {0, 1, . . . , n cut }, G is optimized by brute force. If the optimal transient growth is found at the upper boundary of the considered domains, i.e. for t = t cut , k = k cut or n = n cut , the respective intervals are enlarged in subsequent steps until a local maximum is located in their interior.

Numerical issues
In this section, the performance of the numerical implementation presented in § 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 N 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.

Eigenvalue decomposition
Our discretization of the eigenvalue problem (2.9) has been tested against the results on eigenvalue-critical Reynolds numbers presented in Krueger, Gross & DiPrima (1966, table 2) as well as by replication of the plotted spectra given by Gebhardt & Grossmann (1993, figure 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 λ N 1 as the number N of Legendre or Chebyshev polynomials is increased.
are plotted against N. The test parameters are R Ω = −2, η = 0.5, n = 5 and k = 1 at shear Reynolds numbers Re = 8000 for figure 2(a) and Re = 128 000 for figure 2 The plots in figure 2(a,b) show plateaus of non-convergence for low N, which are due to the difficulty of identifying the respective eigenvalue in a non-converged spectrum. For moderate resolutions (N ∈ [20; 45] in figure 2(a) and N ∈ [45; 90] in 2(b)) spectral accuracy, i.e. exponential convergence rates, is observed for both methods. Notably, however, the convergence turns out to be significantly quicker   TABLE 3. Optimal transient growth G max := sup n,k,t G(t) according to Meseguer (2002, table 1) and present results; parameters are η = 0.881 and N = 50; n max and k max denote the azimuthal and axial wavenumbers that attain optimal transient growth G max .
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 for N = 43 (Re = 8000) and N = 83 (Re = 128 000) compared to N = 62 and N = 104, respectively, in the case of the Petrov-Galerkin scheme (see figure 2). The required resolution N 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 Re O(10 5 ). This also affects the computation of transient growth discussed in the next subsection.

Computation of transient growth
In table 3 our results concerning the optimal transient growth G max := sup n,k,t G(t) for η = 0.881 and the corresponding optimal wavenumbers n max are compared to the numerical data of Meseguer (2002, table 1). The values of k max and G max differ by less than 0.3 %.
The convergence of the maximum transient growth G 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 N ∈ {5, 15, 50} (corresponding to figure 3(a,b,c)), the eigenvalues (top), the evolution of the maximum transient growth G(t) (middle) and the moduli of the components |u r |, |u ϕ | and |u z | of the corresponding optimal perturbation u(0) (bottom) are plotted for comparison. The example parameters are Re = 10 000, R Ω = −2.0, η = 0.8, n = 5 and k = 1.
A few aspects are noteworthy. Around its maximum G is already surprisingly well approximated by only N = 5 Legendre polynomials, whereas the optimal perturbation is far from its actual shape (see figure 3a). For N = 15 (figure 3b) the curve {(t, G(t))} is converged within an error 1 % while its maximum is even approximated up to ≈0.01%. Likewise, the optimal perturbation u(0) 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 N = 15 not to mention N = 5 (top). In fact, it takes as many as N = 50 polynomials for convergence of the two meeting branches (see figure 3c). However, this does not seem to affect the transient growth quantities -even though the converged spectrum in figure 3(c) (top) is even much more stable as a whole than its approximation for N = 15 (figure 3b).
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 O(Re α ) discrete eigenvalues for α ≈ 1 2 , rendering their convergence numerically infeasible for Re O(10 5 ). This agglomeration of eigenvalues can be explained by the dominance of the O(Re) convective multiplicative terms in the linearized Taylor-Couette operator L over the viscous differential contributions in the limit Re → ∞ (see equations (2.6a) and (2.6b)). The asymptotic degeneracy into a pure multiplication operator, i.e. a mapping u → Au with a continuous, matrix-valued function A : [r i ; r o ] → R 3×3 , corresponds to a transition from discrete eigenvalues to a continuous spectrum.
On the other hand, the standard Cholesky decomposition of the matrix M (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 N 51 -just as the crucial meeting point is resolved. Accordingly, one might expect to miss a sudden jump in the maximum transient growth G 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 L 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 N. However, N 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, R Ω , it turns out that resolving the transient growth peak for optimal wavenumbers n = n max , k = k max tends to require the highest resolutions. Moreover, the necessary N are mostly independent of R Ω and at least of the same magnitude for different η. Here greater curvature, i.e. η → 0, results in slower convergence. Consequently, for practical computations, suitable resolutions N Re are determined for different ranges of Re by the convergence of (computationally challenging) test cases, more precisely less than 0.3% deviation in the optimal transient growth for η = 0.2 and N ∈ [N Re − 3; N Re ].
For greater η, lower resolutions N may be sufficient and greater Reynolds numbers than Re = 2 048 000 might be resolvable. However, universal convergence for any parameters R Ω , η, n and k within about 1% may be assumed if N is chosen according to the resulting canonical resolutions N Re given in table 4. They are found to approximately follow a power law of the form N Re = N 0 Re α with N 0 = 2.28 ± 0.06 and α = 0.293 ± 0.002. Starting from these, N 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 M and thus of the transient growth is confined to the stable eigenmodes in Q in agreement with the analysis of Meseguer (2002).
These computation guidelines have been applied to obtain the numerical results presented in § 5.

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

Optimal transient growth in various regimes
According to the numerical strategy discussed in § § 4.2 and 3, the optimal transient growth G max = sup n,k,t G(t) is computed for logarithmically equidistant shear Reynolds numbers 250 Re 2 × 10 6 and η ∈ {0.2, 0.5, 0.8}.
By studying test cases, we find that t ∈ [0; τ 0 ] with τ 0 = 2π/[Re α (1 − η)] and α = 0.85 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 n ∈ {0, 1, . . . , 8} and 0 k 5. Additionally, as discussed in § 3.3, the ranges for t, n and k are enlarged whenever the optimization terminates near one of the upper bounds.
By  figure 4 for orientation, along with the numerically computed (viscous) linear stability boundary. Figure 5(a,b) shows the optimized transient growth G max and the corresponding optimal axial wavenumber k max , respectively, against Re for the considered parameter sets.
The most prominent feature in the double-logarithmic plots of figure 5(a) are the nearly identical asymptotic slopes of the lines in the linearly stable regimes for R Ω ∈ {−3, −1.2, 1.2, 3} showing a characteristic power law G max ∼ Re α with α ≈ 2 3 ± 7 % (compare dashed line in figure 5a). Notably, even in the Rayleigh-unstable counter-rotating regime IV (circles in figure 5a), G max 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 G max (Re) in the different regimes differ only by O(1) 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 G max ∼ Re 2/3 for Re → ∞. Note that this scaling result is in perfect agreement with those by Yecko (2004) obtained for Keplerian flows at fixed R Ω = 1.5 in rotating plane Couette geometry.  figure 5(a,b), respectively. The quasi-Keplerian regime II is shaded for orientation.

Optimal axial wavenumber
Beyond the magnitude of transient growth studied in § 5.1, the spatial structure of the optimal perturbations u max is of great interest. The latter is determined by the optimal axial and azimuthal wavenumbers k max and n max , which attain the optimal transient growth G max shown in figure 5(a). The k max are plotted in figure 5(b) with logarithmic horizontal axes. Note the discontinuities of the curves whenever the (discrete) optimal azimuthal wavenumber n max changes.
The plots reveal a characteristic quasi-two-dimensional columnar structure of the optimal perturbations in regime II (also observed by Yecko (2004)): for Re > Re 0 the optimal transient growth G max (Re) is consistently attained by axially independent perturbations, i.e. k max = 0 (compare R Ω = −3 and R Ω = −1.2 in figure 5b). The transition to k max = 0 typically occurs already for Reynolds numbers as small as Re 0 = O(10 3 ). Only near the Rayleigh line -that is, for −1.2 R Ω < −1 -is a sharp divergence of Re 0 for R Ω → −1 observed. Here k max ≈ 1 holds up to the greatest studied shear Reynolds numbers Re = O(10 6 ).
While k max = 0 is only obtained in the quasi-Keplerian regime (II), in regime I (represented by R Ω = 3(1 − η)/η and R Ω = 1.2(1 − η)/η) k max seems to decay (slowly) to zero for Re → ∞. At least weak axial dependence k max 1 is observed for Re O(10 4 ) in these flows. Once again, the asymptotic decay k max → 0 is most distinct near the solid-body line R Ω → ∞ and is lost near the transition to counter-rotation at Re i = 0. Here an almost constant optimal wavenumber k max = O(1) is observed.
For further illustration, figures 6 and 7 show contour plots of k max in the regimes II and I, respectively. The boundary lines have been computed by a bisecting algorithm with relative accuracy = 10 −2 . The extent of the shaded regions in figure 6 emphasizes the dominance of axially independent columnar perturbations for quasi-Keplerian flows.
In the counter-rotating regime (IV), we observe a growing k max with Re. This difference might be explained by emerging linear instabilities, which first appear for k > 1 in this regime and thus render fully three-dimensional perturbations less dissipative. indeed suggest that, for sufficiently large shear Reynolds numbers, n max depends more on the geometrical parameter η than on Re or R Ω , which parametrize the base flow. In general, an azimuthal wavenumber n seems to be optimal if the associated wavelength is λ ≈ 2πη/(n(1 − η)) ≈ 2, i.e. twice the gap width, leading to vortices that are of about the same radial and streamwise dimension (see e.g. figure 8a, centre right). In contrast, the dominant axial wavenumbers k < 1 in regime I correspond to wavelengths of O(10) rather than O(1) 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 n = 0 and k = π. Moreover, we observe that, the stronger the rotational influence on the fluid's stability expressed by smaller η and/or larger |R Ω |, the smaller are the k max attained for Re → ∞ (figure 5b). 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 n and k in the test cases η = 0.5, R Ω = −2.0 and Re = 1000, 2000, . . . , 128 000 (data not plotted) indeed shows significantly non-columnar modes with k ∼ 5 to be least dissipative for n 1. The principal zero mode n = k = 0, 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 η → {0, 1} is considered. A further study of this parameter is therefore omitted in the following.

Evolution of optimal perturbations
In the sequel, three different optimal perturbations u max,1 , u max,2 and u max,3 are considered at a constant shear Reynolds number Re = 8000 and η = 0.5. The rotation numbers are given by R Ω,1 = −2.0, R Ω,2 = 2.0 and R Ω,3 = 0.8 corresponding to regimes II, I and IV. The optimal wavenumbers are given by k max,1 = 0, k max,2 ≈ 0.464 and k max,3 ≈ 1.200 and n max,1 = n max,2 = n max,3 = 3. The time evolution of these modes is computed by eigenmode decomposition at a polynomial resolution N = 50. Radial-azimuthal and radial-axial projections are shown. In the latter the plots are scaled to show exactly one axial wavelength along the horizontal axis; and, to aid visualization, a unit length (d = gap width) in the axial direction is indicated by the dashed line. The various parts show subsequent snapshots at times t = t j during the transient growth evolution; the t j are also marked in the energy evolution curves plotted in figure 9. Arrow lengths are scaled with the flow velocities whereas their shading from lighter to stronger colours reflects energy densities |u max,i | 2 . The relative rotation of the inner and outer cylinders in the different settings is indicated by arrows visualizing the frequencies Ω o and Ω i , respectively. The corresponding optimal axial wavenumbers are k max,1 = 0, k max,2 ≈ 0.464 and k max,3 ≈ 1.200. Taylor In figure 8 the resulting real parts of u max,1 , u max,2 and u max,3 are shown at a sequence of snapshots t j throughout the transient growth evolution. The flow fields are plotted in radial-azimuthal projection (top) and radial-axial projection (bottom) with z on the horizontal axis except for u max,1 , 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 figure 8(a-c). The colour map from lighter to stronger colours marks regions with relatively low or high energy densities |u max,i | 2 in the current fields.

Transient growth in
The perturbations' total kinetic energy evolution u max,i (t) 2 in relation to the transient growth maxima G max,i are plotted in figure 9 with time scale renormalized by τ 0 = 2π/[Re 0.85 (1 − η)]. The t j considered in figure 8 are identified by markers.
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 2n streamwise elongated vortices. Recalling the different angular velocities Ω i and Ω o of the driving inner and outer cylinders (i.e. Ω i > Ω o > 0 for R Ω = −2.0, Ω o > Ω i > 0 for R Ω = 2.0 and Ω i > 0 > Ω o in the counter-rotating case R Ω = 0.8, 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 figure 8 centre-left and figure 9). As for the axially independent perturbation in figure 8(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 twodimensional 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, Doering & Spiegel (2010) -compare our figure 8(a), centre-left, with figure 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 u max,1 , 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 8a, centre-right), which leads to high dissipation around the energy maximum. On the other hand, the optimal perturbations u max,2 and u max,3 in the regimes I and IV seem to be stabilized in this respect by their axial dependence, leading to 40 % and 115 % larger growth than that attained for R Ω = −2.0 and slower decay in figure 9. This interpretation is supported by the fact that, in spite of the small wavenumber k max,2 ≈ 0.464 of u max,2 , up to 86 % of the kinetic energy is transferred into the axial component around the transient growth maximum. These three-dimensional effects go beyond the classical Orr mechanism.
The characteristic structure of deforming elongated vortices is also reflected in the radial-axial projections in figure 8(b,c). A unique feature of the counter-rotating flow (R Ω,3 = 0.8) is the localization of the optimal perturbation near the inner cylinder walls, where the base flow is locally Rayleigh-unstable. This localization has also been observed in the spiral eigenvectors and in the saturated spiral instability (e.g. Langford et al. 1988). Hence, although the flow remains eigenvalue stable for the chosen parameters, emerging instabilities already seem to interact with non-modal growth mechanisms. This possibly explains the greater energy amplifications in regime IV. 5.4. Transient growth scaling for k = 0 The previous numerical results, especially those for the quasi-Keplerian regime II, motivate the transient growth analysis of axially independent perturbations with k = 0. Moreover, it will be shown in § 6.2 that G k=0 max , i.e. the optimal transient growth of k = 0 perturbations, is indeed independent of the rotation number R Ω .
The parallel slopes for high Reynolds numbers Re O(10 4 ) show a common scaling G k=0 max ∼ Re γ , where the proportionality factor may depend only on η. In order to estimate this, G k=0 max (Re) is computed for logarithmically equidistant Re ∈ [10 5 ; 4 × 10 6 ] for η ∈ {0.05, 0.1, 0.15, . . . , 0.95}. Fits of the form G k=0 max (Re) = a(η)(Re) γ (η) for each η yield exponents γ (η) ≈ 2 3 within an error 0.5 % except for γ (η = 0.2) ≈ 0.657. Hence, a common exponent γ = 2 3 is assumed to be universal and the factor a(η) Transient growth in Taylor Numerically computed optimal transient growth G k=0 max for axially independent perturbations plotted against the shear Reynolds number Re (log-log plot). The curves are independent of R Ω and parallel for Re O(10 4 ), corresponding to a common scaling G k=0 max = a(η)Re 2/3 . (b) Fitted scaling coefficients a(η) for the respective η and high shear Reynolds numbers Re O(10 4 ). The error bars represent data determined by fitting of numerical results for G k=0 max (Re) with errors by mean square deviation. The full line is a third-degree polynomial fit according to (5.1) and (5.2). is independently determined by another fit. The results are plotted in figure 10(b), where the error bars have been determined by the mean square deviation from the data.
In order to obtain an analytical formula for G k=0 max (Re) a third-degree polynomial is fitted to the data in figure 10b taking into account the extremum of a at η = 1, which is due to the system's symmetry with respect to exchanging r i and r o . The result is a 0 ≈ 9.218 × 10 −3 , a 1 ≈ 0.1198 and a 2 ≈ −9.072 × 10 −2 , (5.2) and the corresponding curve is also shown in figure 10(b). Good agreement between fit and data is found, especially for η 0.5, possibly due to the lesser impact of the azimuthal wavenumber's discreteness on the attainable optimal transient growth compared to η < 0.5. For arbitrary η, test cases give less than 7 % error if the analytical formula is applied for Re ∈ [10 4 ; 8 × 10 6 ] and less than 5 % in the interval [10 5 ; 2 × 10 6 ].
The maximum amplification of axially independent perturbations G k=0 max = a(η)Re 2/3 defines a lower bound for the total (k = 0) transient growth G max (Re) in every flow regime. Moreover, the estimate can be expected to hold within a factor of O(1) and is exact in the shaded regions of the quasi-Keplerian regime II in figure 6.

Analytical results for axially independent perturbations
The prominent role played by columnar axially independent perturbations together with their geometrical simplicity motivates an analytical study of their properties, which is pursued in this section. We begin by applying the conjugated curl operator to the linearized Navier-Stokes equation (2.5). This eliminates the pressure gradient terms, yielding  (6.2) For axially independent perturbations (k = 0), the azimuthal velocity u ϕ is determined from u r via the divergence condition and the evolution equations for u r and u z decouple (Gebhardt & Grossmann 1993). Using L rr = L ϕϕ the resulting equations read (in/r)L zz u z −DL zz u z [−(in/r)L rr + D + L rr (ir/n)D + + D + L ϕr + L rϕ D + ]u r   . (6.4) The first and second equations, which are equivalent, determine the evolution of u z : Transient growth in Taylor-Couette flows Using the results D + L ϕr + L rϕ D + = (2B/r 2 )D + − 4in/r 3 and [L rr , r∂ r ] = 2L rr + 2inA =: 2L 0 rr obtained in § A.1, the evolution equation for u r becomes Further using ∂ r (2/r)(rD + rD + − n 2 ) = 2L 0 rr rD + + irn[(2B/r 2 )D + − 4in/r 3 ] (see § A.1) yields This fourth-order partial differential equation is supplemented with the boundary conditions u r (r i ) = u r (r o ) = ∂ r u r (r i ) = ∂ r u r (r o ) = 0, which correspond to the no-slip boundary conditions at the cylinders u r (r i ) = u r (r o ) = u ϕ (r i ) = u ϕ (r o ) = 0.
6.1. Advection of perturbations by the basic flow and universal stability properties A remarkable property of (6.5) and (6.7) is revealed by considering the transformatioñ u r := e inAt u r andũ z := e inAt u z . The derivatives then read ∂ rũ * = e inAt ∂ r u * and ∂ tũ * = e inAt (∂ t + inA)u * , so substituting into (6.5) and (6.7) yields wheref r := (rD + rD + − n 2 )ũ r and f r := (rD + rD + − n 2 )u r . Asũ z andũ r satisfy (6.5) and (6.7) with A = 0, the A dependence of the perturbation's evolution u is entirely described by the factor e −inAt . This factor corresponds to a pure advection of the perturbation with the shear-free uniformly rotating part of the basic flow v B and thus it is locally and globally energy-conserving (|u r | 2 = e inAt e −inAt |ũ r | 2 = |ũ r | 2 ). Although these conclusions might seem obvious at first glance, note that they are not true in the general three-dimensional case k = 0. The minor importance of the parameter A has crucial consequences. Without loss of generality, A = 0 can be assumed when analysing the stability of Taylor-Couette flow to axially independent perturbations. The remaining parameter B characterizing the base flow v B depends only on the shear Reynolds number Re and not on the rotation number R Ω (see (2.10b)), which parametrizes the flow regime. Hence the linear stability of Taylor-Couette flow to axially independent perturbations is independent of R Ω and thus is identical in all regimes. Furthermore, the optimal transient growth G k=0 max for k = 0 provides a lower bound for the absolute maximum G max , which is universal in the sense that it depends only on η and Re. We note that these results can be expected to apply approximately also for weakly axially dependent perturbations in the vicinity of k = 0. 6.2. Global analysis of the evolution equations First, consider the evolution of u z described by (6.5). Operator L zz is the sum of a selfadjoint negative definite operator and a skew hermitian one. As these do not commute, L zz is an example of a non-normal operator, which nonetheless does not allow for transient growth (see § A.2 for a proof). The evolution (6.7) for the radial component u r may be split into two independent problems, ∂ t f r = L rr − 2∂ r 1 r f r and (rD + rD + − n 2 )u r = f r , (6.10) where the second is of Sturm-Liouville type (the solution is given in § A.4) and the first resembles (6.5). Using this factorization, it might be possible to construct an exact analytical solution of the evolution problem (6.7) by incorporating the boundary conditions via an influence matrix method. However, the outer problem in (6.10) remains cumbersome to solve, and, even if one were to write down an expression for an exact solution of (6.7), this would most likely turn out to be too involved to interpret the underlying physics. In the following, the analysis of the evolution equation (6.7) is therefore confined to the limit of asymptotically large Reynolds numbers Re → ∞ and is studied by means of scale analysis. In order to identify and motivate the scales to be studied quantitatively in the WKB analysis of § 7, we consider the energy evolution of a perturbation u = u r e r + u ϕ e ϕ , ∂ t u 2 = 2 Re u, L u = −2Re u, (u · ∇)v B + 2 Re u, r u , (6.11) where the pressure and convective terms drop out as in the derivation of the Reynolds-Orr equation. Using u ϕ = (ir/n)D + u r , the non-normal term in (6.11) becomes N(u) := −2 Re u, (u · ∇)v B = − 4B n Im u r , ∂ r u r , (6.12) while the self-adjoint dissipative summand reads D(u) := 2 Re u, r u =−2{n −2 ( DrD + u r 2 + (n 2 + 1) D + u r 2 ) + Du r 2 − 4 Re u r , r −1 D + u r + (n 2 + 1) r −1 u r 2 }. (6.13) Assume that u r varies on a typical length scale of order O((nRe) −α ) with α > 0. In the limit Re → ∞, the highest-order r derivative dominates in each term of (6.11). As B ∼ Re and u 2 = u r 2 + n −2 rD + u r 2 , we obtain from (6.12) and (6.13) (6.14) According to (6.13), D(u) is strictly negative, so by virtue of (6.14) dissipation always dominates for α > 1 3 . On the other hand, the non-normal term N(u) may be positive, so that, for α 1 3 , growth rates ∂ t ln u 2 = O((nRe) 1−α ) are possible.
Transient growth in Taylor-Couette flows   279 The question remains how long such growth may last. Let us consider a Fouriertype ansatz u r ∼ e imr with wavenumber m = O((nRe) α ). Note that locally this is valid because in the limit Re → ∞ boundary effects are confined to thin layers near the cylinder walls. Then N(u) is of optimal order in (6.14) and N(u) > 0 if and only if n −1 Bm < 0 by virtue of (6.12).
The total velocity field isũ = e i(nϕ+kz) u ∼ e i(nϕ+mr) , so the curves of constant phase (characteristics) are (locally) given by ϕ(r) = ϕ(r i ) − n −1 m(r − r i ). Starting at the inner cylinder the set of these lines form streamwise elongated spiral structures like the vortices in figure 8. To attain growth they have to be oriented according to the sign sgn(∂ r ϕ) = −sgn(n −1 m) = sgn(B) = −sgn (∂ r Ω).
(6.15) Thus, the characteristics of the perturbations have to be misfit to the base flow's angular velocity profile Ω B = r −1 v B ϕ , as observed in the numerical computations of § 5.3. Therefore, energy amplification may only occur transiently until the perturbation has been sheared into the 'fit' orientation by advection, analogously to the perturbation dynamics associated with the Orr mechanism in channel flows (see e.g Farrell 1988). Within the advective time scale T = O(Re −1 ), i.e. a cylinder rotation period, the shear uniformly distorts the flow profile between the inner and outer cylinders by a length of order O(1). Consequently, as the initial streamwise elongation of the characteristics is O(n −1 m) and m = O((nRe) α ), the time t 0 for the perturbation to be tilted into the fit direction is Viscosity prevents transient growth if α > 1 3 . Now assume u is an optimal perturbation for α < 1 3 . Then we can evolve this mode backwards until times of order O((nRe) −2/3 ) before its energy maximum, introduce the result as a new initial condition and thereby attain additional growth. Thus, optimal perturbations must vary on length scales O(nRe) −1/3 and t 0 = O((nRe) −2/3 ) is the natural time scale for transient growth.
Our numerical computations are in perfect agreement with these scaling results. However, we cannot obtain an analytical estimate for the optimal transient growth with this section's zeroth-order approach. Therefore, in the next section we introduce the time scale t 0 into the evolution (6.7) and analyse it by means of a first-order WKB approximation. Our analysis closely follows the work of Chapman (2002, pp. 47-53) for oblique modes in channel flows.

WKB ansatz
We make a WKB ansatz with amplitudeã and rapidly oscillating phase δ −1 φ, where bothã and φ depend ont and r. Together with the divergence condition, this yields u ϕ = (ir/n)D + u r = O(∂ r u r ) = O(δ −1 u r ). Hence the scalingã = δa, with a = O(1) and φ = O(1), is required in order that initial perturbations u = u r e r + u ϕ e ϕ have unit energy norm ( u(0) 2 = 1). We now substitute the WKB ansatz (7.2) into the evolution equation (7.1). Because a, φ = O(1), the evolution equation needs to be independently satisfied at each order in δ. At leading order O(δ −1 ) the equation reduces to (see § A.3) By using this solution to eliminate δ −1 terms in (7.1), we obtain Defining τ := i(∂ r φ) and ∂t = i(∂t∂ r φ)∂ τ = −(2B 0 /r 3 )∂ τ (Chapman 2002, p. 49) yields According to this solution, a becomes singular for τ → 0, which may raise doubts about its physical correctness. However, in this limit, the underlying separation of orders in the WKB approximation breaks down so that O(δ 1 ) terms in (7.4) or even in the leading-order equation have to be considered. These bound the blow-up, leading to an overall nearly singular amplitude behaviour in the complete linearized dynamics given by (6.7). In numerical simulations, this manifests itself in increasingly sharp peaks of the optimal perturbation's energy for Re → ∞, as visualized in figure 11. The larger Re, the longer the blow-up seems to follow the singular WKB solution (7.6) before the energy growth is capped near the maximum blow-up timet 0 . Most prominently, for Re = 1 024 000 it is only in a neighbourhood (1 ± 0.05)t 0 about the maximum that the singularity is smoothed out by additional terms, resulting in the sharpest peak in figure 11.

Construction of optimal perturbations
Assume that the amplitude's growth according to equation (7.6) is capped as soon as the next-order terms become relevant. Then the optimal energy growth is attained if: (i) the blow-up occurs at a common timet 0 over the whole radial domain r ∈ (r i , r o ); and (ii) the O(δ 1 ) terms in (7.4) are of the highest attainable order in τ . Condition (i) ensures that no averaging effects of the spatial integral evaluated for the computation of u(t) 2 limit the global energy maximum in time. It is equivalent to ∂ r φ(r,t 0 ) = ∂ r φ 0 (r) + (2iB 0 /r 3 )t 0 = 0, so that φ 0 = (iB 0 /r 2 )t 0 + c and without loss of generality φ = −(iB 0 /r 2 )(t −t 0 ).
(7.8) Equation (7.8) is also satisfied for φ 0 = (iB 0 /r 2 )t 0 + c. Hence, this is indeed the optimal initial phase giving the optimal perturbation according to WKB theory, u r = δa exp φ δ eq. (7.6) δr 2 (t −t 0 ) . (7.9)  by Chapman (2002, p. 51ff.). By equation (7.9) this should yield δr 2 B 0 Im(ln u r (r, 0)) =t 0 + O(δ). (7.13) Owing to the non-uniqueness of the complex logarithm, relation (7.13) needs to be assumed to be satisfied for some r 0 ∈ (r i , r o ). We choose r 0 = 1 2 (r i + r o ). In figure 12 the blow-up timet 0 computed from (7.13) is plotted against the radial coordinate r (solid curve). This WKB prediction is compared for Reynolds numbers Re ∈ {10 3 , 10 4 , 10 5 , 10 6 }, corresponding to δ ∈ {0.069, 0.032, 0.015, 0.007}, to the optimal time determined numerically from the full equations (dashed line). The expected error ranges are denoted by [t 0 + δ;t 0 − δ] (dash-dotted lines). Excellent agreement between the numerical results and WKB solution within the predicted error of order δ and convergence for Re → ∞ is found. Significant deviations are confined to a O(δ) neighbourhood of the cylinder walls in which growth is prevented a priori by the boundary conditions. Hence, the initial phase's behaviour as a key property of the derived WKB approximation has been numerically verified.

Discussion
Rayleigh-stable Taylor-Couette flows with the outer cylinder rotating faster than the inner one tend to become turbulent at moderate Reynolds numbers Re = O(1000) (Taylor 1936;Borrero-Echeverry & Schatz 2010;Burin & Czarnocki 2012). In the case of the quasi-Keplerian regime II, where the inner cylinder rotates faster than the outer one, the existence of turbulence remains debated (Ji et al. 2006;Paoletti & Lathrop 2011). At the same time, Rayleigh-unstable but linearly (eigenvalue) stable counter-rotating Taylor-Couette flows are known to undergo subcritical transition (Coles 1965).
In this work, the optimal linear transient growth G max , i.e. the maximum nonnormal energy amplification of infinitesimal perturbations, has been investigated. Our analysis covers the whole parameter regime of Taylor-Couette flow, spanned by the shear Reynolds number Re, the cylinder radius ratio η and the rotation number R Ω . We find that accurate transient growth computations are numerically feasible up to Re = O(10 6 ), even though the characteristic Y-shaped eigenvalue spectrum of the linearized Navier-Stokes operator cannot be resolved for such Reynolds numbers. This is in contrast to previous studies of channel flow (e.g. Reddy & Henningson 1993), which suggest that resolving the Y shape of the spectrum is necessary to accurately compute transient growth. For Taylor-Couette flow the transient growth maximum G max is well converged for resolutions where the approximated spectrum is still far from its natural shape. This allows us to examine the optimal transient growth for large Re. Our numerical computations show an asymptotic scaling G max ∼ Re α for Re O(10 4 ) with α ≈ 2 3 for all geometries considered, η ∈ {0.2, 0.5, 0.8}, and all linearly stable flows. This reveals energy growth of the same order in all regimes and allows for arbitrary transient amplifications if Re is sufficiently large. Moreover, the dynamics discussed in § 5.3 suggest that the underlying growth mechanisms (interpreted here as a curved analogue of the Orr mechanism) are essentially the same in the studied regimes I, II and IV. In the counter-rotating regime IV, there are additional amplifying effects of the Rayleigh instability. Notably, the observed spiral-shaped structures resemble those of the unstable eigenmodes emerging in the case of an imposed radial inflow at the rotating outer cylinder, reported by Gallet et al. (2010). A distinction between the regimes is found in the optimal axial wavenumber k max , which reflects the axial dependence of the optimal perturbations attaining maximum energy amplification. Although columnar structures, representing axially invariant modes, dominate within practically the whole quasi-Keplerian regime (II) above Re = O(1000) corresponding to k max = 0, weakly three-dimensional optimal perturbations 0 < k max < 1 are found in the likewise Rayleigh-stable regime I for Re → ∞. The reason why a weak axial structure enhances transient growth in the latter, but not in the former, remains open. For counter-rotating flows, greater k max = O(1) turn out to attain higher energy maxima.
Our numerical results reveal an important role of axially invariant perturbations for transient growth in linearly stable Taylor-Couette flow. Hence, the corresponding linearized Navier-Stokes equations have been studied analytically in § § 6 and 7. Firstly, the analysis has revealed that transient growth and linear stability are indeed independent of R Ω in the case k = 0. Then we have shown that optimal perturbations blow up and decay by the Orr mechanism within the time scale t 0 = O((nRe) −2/3 ). By introducing this scale in the linearized evolution equations, an optimal transient growth scaling G k=0 max (Re) = a(η)Re 2/3 for axially independent perturbations has been derived in the limit Re → ∞, following the channel flow WKB analysis of Chapman (2002). The results apply for all R Ω and thus in all flow regimes. For the coefficient a(η) the semi-empirical formula, given by (5.1) and (5.2) has been obtained by a cubic fit to the numerical data.
The expression G k=0 max (Re) = a(η)Re 2/3 provides a universal lower bound for the optimal transient growth of general three-dimensional perturbations. This bound attains the optimum in most of regime II according to the numerical results. However, while quasi-Keplerian flows thus indeed have the smallest possible energy amplification potential, the growth is nevertheless of the same order as in the other regimes. Temporary amplifications of disturbances may promote nonlinear instability if growing modes are consistently fed by nonlinear energy redistribution. Hence, by our scaling results, such a transient growth-mediated instability is as likely to exist in quasi-Keplerian flows as in any other regime. However, axially independent perturbations are possibly not equally fit to feed nonlinear instabilities as three-dimensional ones, e.g. because of their sharper growth and decay. In the future this question could be addressed by studying nonlinear generalizations of transient growth, such as applied for instance by Pringle & Kerswell (2010), Monokrousos et al. (2011) andPringle, Willis &Kerswell (2012). On the other hand, such investigations are computationally expensive and beyond the scope of the present work. Meseguer (2002) found a strong correlation between the experimentally observed nonlinear stability boundary (Coles 1965) and optimal transient growth G max in counter-rotating flows. Following these ideas, we estimate the threshold shear Reynolds number Re T for subcritical transition in quasi-Keplerian flows using our universal scaling result. To this end G max was computed numerically at the subcritical stability boundary of Taylor-Couette flow (results not shown) according to measurements by Mallock (1896), Wendt (1933), Taylor (1936), Coles (1965), Borrero-Echeverry & Schatz (2010), Burin & Czarnocki (2012) and Avila & Hof (2013). Not surprisingly, the correlation is not as strong as observed by Meseguer (2002), who only considered the data of Coles (1965). Moreover, Burin & Czarnocki (2012) have found their experimental results to depend significantly on the applied endcap configurations, where the sensitivity is stronger for wider gaps. Our results indeed range from G max ≈ 54 to G max ≈ 155. If we translate this to shear Reynolds numbers, the uncertainty roughly agrees with the observed endcap effects. Calculating the mean value of all computed threshold amplifications yields an a priori estimate for the threshold transient growth in an arbitrary Taylor-Couette flow setting of G max,T = 92 ± 26.
Applying the estimate formula for G max , we obtain a threshold Reynolds number of Re T = a(η) −3/2 (880 ± 370), giving, for instance, Re T = 67 000 ± 29 000 if η = 0.7. For quasi-Keplerian flows, recent experiments have proceeded up to Re = O(10 6 ), yielding contradictory results (see Ji et al. 2006;Paoletti & Lathrop 2011). However, Avila (2012) has shown the state-of-the-art Taylor-Couette apparatus to be possibly unsuited for such measurements because of axial endwall effects. On the other hand, our estimated Re T still lies within the range of direct numerical simulations. Hence, these may be able to resolve the controversy concerning the existence of hydrodynamic turbulence in the quasi-Keplerian regime. If turbulence were found, the value of the threshold Re T could be used to probe the significance of linear transient growth as a measure for subcritical instability.