Transient growth of perturbations on scales beyond the accretion disc thickness

Abstract

A turbulent state of spectrally stable shear flows can be developed and sustained according to the bypass scenario of transition. If it works in non-magnetized boundless and homogeneous quasi-Keplerian flow, then transiently growing shearing vortices should supply turbulence with energy. Employing the large shearing box approximation, as well as a set of global disc models, we study the optimal growth of the shearing vortices in such a flow in the whole range of azimuthal length-scales, λ_y, compared to the flow scaleheight, H. It is shown that with the account of the viscosity the highest possible amplification of shearing vortices, G_max, attains maximum at λ_y ≲ H and declines towards both the large scales λ_y ≫ H and the small scales λ_y ≪ H in good agreement with analytical estimations based on balanced solutions. Our main attention is directed to the large-scale vortices λ_y ≫ H, which produce G_max ∝ (Ω/κ)⁴, where Ω and κ denote local rotational and epicyclic frequencies, respectively. It is demonstrated that the large-scale vortices acquire high-density perturbation as they approach the instant of swing. At the same time, their growth is not affected by bulk viscosity. We check that G_max obtained globally is comparable to its local counterpart, and the shape and localization of global optimal vortices can be explained in terms of the local approach. The obtained results allow us to suggest that the critical Reynolds number of subcritical transition to turbulence in quasi-Keplerian flow, as well as the corresponding turbulent effective azimuthal stress should substantially depend on shear rate.

accretion, accretion discs, hydrodynamics, instabilities, turbulence, proto-planetary discs

1 INTRODUCTION

Homogeneous quasi-Keplerian shear flow, which is defined as a rotating flow with its specific angular momentum and angular velocity, respectively, increasing and decreasing outwards, has proved to be the most enigmatic case with regards to the transition to turbulence. To date, its hydrodynamic (HD) stability with respect to finite amplitude perturbations has been experimentally checked up to the Reynolds number equal to a few millions (see Edlund & Ji 2014). This is a high value, bearing in mind that the Reynolds number corresponding to the subcritical transition is R_T ∼ 350 for plane Couette flow and R_T ∼ 1000 for plane Poiseille flow – see Schmid & Henningson (2001), as well as numerous references therein – and it is R_T ∼ 7000 for the narrow gap rotating Taylor–Couette flow with angular velocity increasing towards the periphery (see Burin & Czarnocki 2012). However, there are a number of well-known cases of observational evidence that cold astrophysical discs, where magnetorotational instability (MRI) ceases to operate, are accreting. At the same time, the smallness of microscopic viscosity entails huge Reynolds numbers up to ∼10¹³ attained in astrophysical discs. Thereby, the possibility of a purely HD route to turbulence in quasi-Keplerian flow is not ruled out and still represents an important unresolved problem.

The transition to turbulence and the replenishment of turbulent energy budget cannot both take place without growth of perturbations (i.e. an increase of their amplitudes due to the energy transfer from the background shear). At the same time, it turns out that instantaneous growth rate is independent of perturbation amplitude (Henningson & Reddy 1994). This basic fact is obtained from the Reynolds–Orr equation for evolution of kinetic energy of perturbations obeying Navier–Stokes equations. Thus, the energy transfer from the background motion to turbulent fluctuations is given by mechanisms that are present in linear dynamical equations. For a subcritical transition to turbulence, which is the case also for a spectrally stable homogeneous quasi-Keplerian flow, this process has been attributed to a transient growth of perturbations, which remains the only way to draw energy from the background motion in the absence of exponentially growing modes of perturbations (Trefethen et al. 1993; Schmid 2007).

As has become clear in the past decades, there are two different kinds of transiently growing perturbations and two respective variants of growth mechanisms in homogeneous shear flows. The first are usually referred to as streamwise rolls (Butler & Farrell 1992), which generate growing streamwise streaks (i.e. the domains of the enhanced streamwise velocity) by the so-called lift-up effect (Ellingsen & Palm 1975). The lift-up effect plays a key role in the self-sustaining process (SSP), which provides a turbulent state of plane shear flows (see Waleffe 1997). The bypass transition via the transient growth of streamwise rolls has been investigated in depth in the literature (see, e.g. Reddy et al. 1998, and references therein), because of the possibility of simulating the breakdown to turbulence in plane shear flows at moderate Reynolds numbers. Recently, Mamatsashvili et al. (2016) has studied turbulence dynamics by means of the three-dimensional (3D) Fourier representation, enabling one to describe the SSP in terms of the non-linear transverse cascade, which is qualitatively different from the previously known direct (or inverse) turbulence cascade. Yet, the situation becomes much less clear as we proceed to spectrally stable rotating shear flows. Basically, the Coriolis force that comes into play stabilizes the rolls (i.e. their transformation to streaks turns into epicyclic motion of fluid particles), which prevents transient growth. In such a way, another type of transiently growing perturbations is recognized: shearing vortices.

Shearing vortices are streamwise elongated but streamwise-dependent structures enhanced via contraction by the shear, provided that vorticity is conserved. The latter is usually called the Orr mechanism or, alternatively, the swing amplification process. In order to be amplified in the rotating flow, shearing vortices must be in the form of leading (trailing) spirals with respect to the background flow of quasi-Keplerian (cyclonic) type; see, e.g. Lesur & Longaretti (2005) for a definition of rotating flow types. Thus, only certain spatial Fourier harmonics (SFHs) of shearing vortices can supply turbulence with energy. Note that shearing vortices are also incorporated in the bypass transition of two-dimensional (2D) models of plane shear flows, where the rolls are absent. Another picture of the non-linear feedback sustaining turbulence should be suggested in this case (see Chagelishvili et al. 2003). The novel ingredient here is a coupling of SFHs from certain quadrants of wave-vector space, which repopulates the quadrant with linearly growing shearing vortices. Lithwick (2007) has considered the particular case of such a coupling between axisymmetric and shearing vortices in a 2D model of Keplerian flow and has shown that it produces a new leading spiral with a larger amplitude than the existing one. In light of the present study, it is important to note that this kind of positive non-linear feedback is provided by the shearing vortices of a sufficiently large azimuthal length-scale, λ_y. Moreover, in later work, Lithwick (2009) has found that non-linear vortices emerging by this 2D mechanism survive in a 3D model, provided that their azimuthal length-scales are larger than the flow scaleheight, λ_y > H. Umurhan & Regev (2004) and Horton et al. (2010) also studied the non-linear 2D dynamics of rotating (and plane) shear flows, providing additional arguments in favour of the non-linear sustenance of linearly growing shearing vortices. However, presently, there is no theoretical understanding about whether the reservoir of linearly growing shearing vortices can be repopulated in the full 3D dynamics of perturbations with not necessarily λ_y > H. Nevertheless, Meseguer (2002) has revealed a strong correlation between a value of maximum transient growth, G_max, and a point of relaminarization of the Rayleigh-stable incompressible flow between counter-rotating cylinders. Apparently, this is an argument in favour of bypass transition via the swing amplification of perturbations, as G_max in this particular regime of Taylor–Couette flow is produced by shearing vortices rather than by streamwise rolls.

The quantitative characteristic of transient growth, G_max, is obtained employing the method of optimization for solving the problem of perturbation dynamics. Yecko (2004) and Mukhopadhyay, Afshordi & Narayan (2005) (hereafter MAN05) have also identified the shearing vortices, as well as the corresponding G_max, applying this method to the dynamics of 3D incompressible HD perturbations. They have considered the problem in the shearing box approximation invented earlier by Goldreich & Lynden-Bell (1965) in application to swing amplification in perturbed stellar discs. It has been confirmed that in the case of Keplerian shear, the optimal growth is attained by quasi-2D columnar leading spirals exposed to swing amplification. Additionally, MAN05 have considered super-Keplerian shear rates, which are intermediate between the Keplerian shear and the one with constant specific angular momentum distribution. They have found that exactly constant angular momentum shear provides the optimal vortices in the form of streamwise streaks, which are deformed by the background flow and give birth to growing streamwise rolls. This mechanism of non-modal dynamics is exactly the reverse to the lift-up effect in plane shear flows and, consequently, it was called the ‘anti-lift-up effect’ (Antkowiak & Brancher 2007). The corresponding maximum optimal growth factor is several orders of magnitude larger than that for the Keplerian shear. Another important result is that the maximum optimal growth suffers an exceedingly rapid drop as one proceeds to the shear flow with a non-zero (positive) angular momentum radial derivative; see fig. 4 of MAN05. Physically, this is caused by the fact that the growth of rolls degenerates into epicyclic motions, whereas the growth of small-scale shearing vortices is much lower and weakly depends on the shear rate. Further, the recent results of Maretzke, Hof & Avila (2014) confirm that in all spectrally stable regimes of Taylor–Couette flow, rather than only in the quasi-Keplerian regime, the incompressible 3D perturbations attain the highest growth factors by virtue of swinging shearing vortices. However, though in each regime they find growth factors of the same order of magnitude, there are qualitative features of cyclonic and counter-rotating flows that make them different from the quasi-Keplerian regime. Unlike in the quasi-Keplerian regime, in the cyclonic and counter-rotating regimes, those shearing vortices have (though weak) axial dependence and up to 80 per cent of kinetic energy transferred into the axial velocity perturbation near the transient growth maximum. This feature goes beyond the swing amplification mechanism and Maretzke et al. (2014) note that the reason for such a qualitative difference between quasi-Keplerian and other regimes remains unclear. Thereby, the two-dimensionality of shearing vortices attaining maximum transient growth is a unique peculiarity of the quasi-Keplerian flow, which might prove to be the reason for its extraordinary non-linear HD stability.

The studies mentioned above have dealt with incompressible perturbation dynamics, which in the context of accretion disc physics implies that perturbations have a characteristic length-scale much smaller than the disc thickness. Hereafter, we refer to such perturbations as ‘small-scale vortices’. In contrast, we consider how a finite disc thickness affects shearing vortices in the quasi-Keplerian flow. Zhuravlev & Razdoburdin (2014) (hereafter ZR14) showed analytically that the growth factors of the vortical leading spirals with azimuthal wavelength larger than the disc thickness, which we refer to as ‘large-scale vortices’ hereafter, should strongly depend on the shear rate as one approaches constant angular momentum rotation, and G_max ∝ κ⁻⁴, where κ is epicyclic frequency expressed in units of rotation frequency. In this paper, we check this analytical result numerically by employing a variational technique for the determination of linear optimal perturbations in the shearing box approximation. Our task is to take into account the dynamical influence of viscosity. Besides, we also consider the influence of bulk viscosity, which is especially important to do as the spiral swings from a leading to a trailing spiral, because as this happens, the radial wavenumber of SFH tends to zero and it is no longer possible to distinguish between vortices and density waves. Thereby, the field of the velocity perturbation acquires significant divergence, which might cause a suppression of transient growth of initially vortical SFH. Also, as far as the azimuthal wavenumber of SFH is of the order of the inversed disc thickness, the vortical SFH generates strong density wave at the instant of swing (Chagelishvili et al. 1997; Bodo et al. 2005; Heinemann & Papaloizou 2009), whilst the acoustic motion of the fluid is primarily affected by a finite time of the establishment of thermodynamic equilibrium, characterized by the bulk viscosity (see, e.g. Landau & Lifshitz 1987).

Whereas, in this work, such a problem is considered for the first time in application to rotating and, in particular, quasi-Keplerian flows, Hanifi, Schmid & Henningson (1996), Farrell & Ioannou (2000) and Malik, Alam & Dey (2006) have already studied it with an interest in compressible plane shear flows. They found that the lift-up effect remains valid with the account of a finite Mach number, and that the streamwise rolls remain the optimal perturbations in the flow. With regards to the shearing vortices in a compressible medium, Farrell & Ioannou (2000) have found their growth factors to be considerably larger than the corresponding incompressible values, as explained by the strong excitation of density waves by shearing vortices at high Mach numbers. Note, however, that none of these authors considered the influence of bulk viscosity, either setting it to zero or adopting the Stokes hypothesis.

This paper is organized as follows. We start describing the 2D model for local perturbations of quasi-Keplerian shear flow in Section 2. Next, in Section 3, we introduce the norm of perturbations and the definition of maximum optimal growth, G_max, to be determined further. We describe an iterative loop employed to obtain G_max. It consists of basic and adjoint dynamical equations – see equations (7)–(9) and (17)–(19), respectively. Next, we proceed with a numerical study of the optimal growth. First, we analyse the suppression of G_max by bulk viscosity and we find that shearing vortices are almost unaffected by the value of bulk viscosity (see Section 4.1). At the same time, the growth factors of large-scale vortices appear to be substantially larger than those of small-scale vortices. Secondly, in Section 4.2, we present the optimal growth in super-Keplerian shear demonstrating that G_max strongly depends on the shear rate, increasing as one proceeds to the constant specific angular momentum rotation. Moreover, this numerical result (see Fig. 3) is in good agreement with analytical estimations, which are relegated to Appendix B. In Section 4.3, we discuss possible consequences of our results, speculating that the transition Reynolds number, as well as the dimensionless azimuthal turbulent stress, may depend very much on the shear rate. Additionally, in Section 5, we consider the optimal growth of large-scale shearing vortices on a global spatial scale, when the azimuthal wavelength of perturbations becomes comparable to a radial disc scale. It is shown that a global initial vortex is a radially localized leading spiral, having a certain radial position and a degree of winding, both altering with the optimization time-span. Employing the general condition (60), we show that these properties of the leading spiral can be explained well in terms of the local approach, even for perturbations with the azimuthal wavenumbers of order of unity. However, the value of globally determined G_max substantially depends on particular radial profiles of surface density and viscosity in a flow with alternating shear rate (e.g. in a relativistic disc). We give a summary of our study and briefly discuss further developments in Section 6.

2 MODEL FOR DYNAMICS OF PERTURBATIONS

In order to study the transient dynamics of perturbations, taking into account a finite disc scaleheight, we employ the most basic model of a perturbed flow. As we noted in the introductory section, the most rapidly growing vortices in quasi-Keplerian flows have columnar structure, and the swing amplification of spirals occurs via the dynamics in the disc plane. This is why we consider the perturbations, which are independent of the disc vertical direction. Next, we leave out the details of energy balance between viscous heat and radiative (or/and thermal conductivity) cooling of fluid elements, adopting a polytropic equation of state with polytropic index n for both background and perturbed flows. Note that polytropic perturbations with a velocity field independent of vertical coordinate preserve hydrostatic equilibrium (i.e. they have a planar velocity field). Thus, the 3D problem with such a restriction is reduced to a 2D problem after one integrates equations along the vertical coordinate. However, one should bear in mind that vertical and planar perturbed motions can be fully separated from each other in the particular cases of vertically isothermal or vertically homogeneous flows only. Otherwise, even the initially columnar perturbation may acquire vertical dependence and, consequently, vertical velocity perturbation during its evolution. This issue must be addressed in our future study of optimal transient growth of shearing vortices in vertically inhomogeneous baratropic quasi-Keplerian flow, from which the restriction of perturbation columnar structure is going to be removed; see also the discussion in Section 6.1.

Substituting the polytropic relation between Eulerian perturbations of pressure and density, p₁ and ρ₁, as well as their background counterparts, p and ρ, into the continuity equation, we change to a new variable W = p₁/ρ in the dynamical equations for linear perturbations. These equations are written in the large shearing box approximation (see A.3 in Umurhan & Regev 2004) and integrated across the disc thickness. We finally obtain

\begin{equation} \left(\frac{\mathrm{\partial} }{\mathrm{\partial} t} - q\Omega _0 x\frac{\mathrm{\partial} }{\mathrm{\partial} y} \right) u_x - 2\Omega _0 u_y = -\frac{\mathrm{\partial} W}{\mathrm{\partial} x} + f_x + g_x, \end{equation}

(1)

\begin{equation} \left(\frac{\mathrm{\partial} }{\mathrm{\partial} t} - q\Omega _0 x\frac{\mathrm{\partial} }{\mathrm{\partial} y} \right) u_y + (2 - q)\Omega _0 u_x = -\frac{\mathrm{\partial} W}{\mathrm{\partial} y} + f_y + g_y, \end{equation}

(2)

\begin{equation} \left(\frac{\mathrm{\partial} }{\mathrm{\partial} t} - q\Omega _0 x\frac{\mathrm{\partial} }{\mathrm{\partial} y} \right) W + a_*^2 \left(\frac{\mathrm{\partial} u_x}{\mathrm{\partial} x} + \frac{\mathrm{\partial} u_y}{\mathrm{\partial} y} \right) = 0, \end{equation}

(3)

where u_x, u_y are the Eulerian perturbations of the velocity components excited in a small patch of the disc. Variables x and y are Cartesian coordinates, which locally correspond to radial and azimuthal directions in the disc, whereas Ω₀ is the angular velocity of the fluid rotation at x = 0 corresponding to some radial distance r₀ ≫ x. Also, q is a constant shear rate that defines the background azimuthal velocity as v_y = −qΩ₀x. All other background quantities, such as Σ ≡ ∫ρ dz and a_*, are assumed to be constant. Note that |$a_*^2=n a_{\rm eq}^2/(n+1/2)$|⁠, where a_eq is the sound speed in the equatorial plane of disc. The inviscid part of equations (1)–(3) coincides with the set of equations (28)–(30) given in ZR14. At the same time, the new ingredients are viscous terms, which we take as

\begin{equation} f_{x,y} = \nu \left(\frac{\mathrm{\partial} ^2}{\mathrm{\partial} x^2} + \frac{\mathrm{\partial} ^2}{\mathrm{\partial} y^2} \right) u_{x,y}, \end{equation}

(4)

\begin{equation} g_{x,y} = \left(\nu _{\rm b} + \frac{\nu }{3}\right) \frac{\mathrm{\partial} }{\mathrm{\partial} x,\mathrm{\partial} y} \left(\frac{\mathrm{\partial} u_x}{\mathrm{\partial} x} + \frac{\mathrm{\partial} u_y}{\mathrm{\partial} y} \right), \end{equation}

(5)

where ν and ν_b denote kinematic viscosity coefficients corresponding, respectively, to the shear and bulk viscosities of the background flow, η and ζ, integrated across the disc thickness and subsequently normalized by Σ.

In order to obtain f_{x, y} and g_{x, y} in the form of equations (4) and (5), we first collect perturbed parts of all dissipative terms entering the Navier–Stokes equations; see, for example, Appendix B in Kato, Fukue & Mineshige (2008). Note that those terms emerge not only due to viscous forces but also as a consequence of radial advection in the background flow. However, in the case of local perturbations with characteristic length λ ≪ r₀, only the dissipative terms containing the highest-order spatial derivatives of perturbations should be retained. This suggests that all contributions containing νρ₁/ρ or η₁/ρ, where η₁ is the Eulerian perturbation of dynamical viscosity, or proportional to the radial background velocity v_x ∼ ν/r₀ should be omitted. Indeed, for example, the advection term

\begin{equation*} v_x \frac{\mathrm{\partial} u_x}{\mathrm{\partial} x} \sim \nu \epsilon _u (r_0 \lambda )^{-1}, \end{equation*}

whereas one of the terms from equation (4)

\begin{equation*} \nu \frac{\mathrm{\partial} ^2 u_x}{\mathrm{\partial} x^2} \sim \nu \epsilon _u (\lambda )^{-2}, \end{equation*}

which is r₀/λ times larger (ε_u is an amplitude of u_x). In this way, we obtain expressions (4) and (5).

2.1 Dimensionless equations for single SFH

Shear flow with a_* < ∞ has a finite scaleheight H = a_*/Ω₀, which follows from the condition of vertical hydrostatic equilibrium in an unperturbed disc with a polytropic equation of state. Following ZR14, let us use the natural H in our choice of dimensionless comoving Cartesian coordinates:

\begin{equation} x^\prime = \frac{x}{H},\qquad y^\prime = \frac{y+q\Omega _0 xt}{H},\qquad t^\prime =\Omega _0 t. \end{equation}

(6)

The coordinates (6) lead us to a spatially homogeneous set of equations, which in the case of boundless shear have partial solutions in the form of SFH, |$f = \hat{f} (k_x,k_y,t^\prime ) \exp ({\rm i} k_x x^\prime + {\rm i} k_y y^\prime ),$| where f is any of the unknown quantities, |$\hat{f}$| is its Fourier amplitude and (k_x, k_y)¹ are the dimensionless wavenumbers along x and y axes, expressed in units of H⁻¹. Omitting the primes from now on and introducing the shearing radial wavenumber |$\tilde{k}_x \equiv k_x + q k_y t$| as well as the full wavenumber squared |$k^2 \equiv \tilde{k}_x^2 +k_y^2$|⁠, we arrive at the following dimensionless equations:

\begin{eqnarray} \frac{{\rm d} \hat{u}_x}{{\rm d} t} &=& 2\hat{u}_y - {\rm i}\tilde{k}_x \hat{W} - R^{-1} k^2 \hat{u}_x - \left(\frac{R^{-1}}{3} + R_{\rm b}^{-1}\right)\nonumber \\ &&\times\, \tilde{k}_x (\tilde{k}_x \hat{u}_x + k_y \hat{u}_y), \end{eqnarray}

(7)

\begin{eqnarray} \frac{{\rm d} \hat{u}_y}{{\rm d} t} &=& -(2-q)\hat{u}_x - {\rm i} k_y \hat{W} - R^{-1} k^2 \hat{u}_y - \left(\frac{R^{-1}}{3} + R_{\rm b}^{-1}\right)\nonumber \\ &&\times \,k_y (\tilde{k}_x \hat{u}_x + k_y \hat{u}_y), \end{eqnarray}

(8)

\begin{eqnarray} \frac{{\rm d} \hat{W}}{{\rm d} t} &=& - {\rm i} (\tilde{k}_x \hat{u}_x + k_y \hat{u}_y ). \end{eqnarray}

(9)

In equations (7)–(9), it is assumed that |$\hat{u}_x,\hat{u}_y$| are expressed in units of a_* and that |$\hat{W}$| is expressed in units of |$a_*^2$|⁠. Also, there are two different Reynolds numbers in equations (7)–(9), R and R_b, which parametrize dynamical and bulk viscosities, respectively. Explicitly,

\begin{equation} R\equiv a_*^2/(\nu \Omega _0), \qquad R_{\rm b} \equiv a_*^2 / (\nu _{\rm b} \Omega _0). \end{equation}

(10)

In Appendix A1, we describe the incompressible limit of the model presented above. Note that in this limit the wavenumbers are expressed in units of inversed auxiliary length, rather than in units of H⁻¹. Also, the dynamical problem becomes of the first order; that is, for particular values of k_x, k_y, q and R, one uniquely determines a normalized solution and the corresponding transient growth factor (see equations A7 and A9). On the contrary, the problem (7)–(9) has the third order, which is why one first has to find a particular solution of equations (7)–(9) corresponding to the highest transient growth factor for a specified time interval. This solution is referred to as the optimal solution. To find this optimal solution, in Section 3, we supply the set of equations (7)–(9) with the corresponding set of adjoint equations. Afterwards, we will be able to implement the iterative loop elucidated in Section 2 of ZR14.

3 TRANSIENT GROWTH FACTOR

In this section, we describe how the transient growth of perturbations is quantified in our study. As soon as finite compressibility is taken into account, the size of perturbations is commonly measured by surface density of their acoustic energy, E_ac, resulting in the following expression for dimensionless quantities constituting the SFH:

\begin{equation} E_{\rm ac} = \frac{1}{2} (|\hat{u}_x|^2 + |\hat{u}_y|^2 + |\hat{W}|^2). \end{equation}

(11)

The growth factor of an arbitrary SFH is then given by

\begin{equation} g(t) \equiv \frac{E_{\rm ac}(t)}{E_{\rm ac}(0)}. \end{equation}

(12)

Equation (11) implies that in the functional space of perturbations, with |$\boldsymbol {q}$| and |$\tilde{\boldsymbol {q}}$| being its members, the inner product |$(\boldsymbol {q}, \tilde{\boldsymbol {q}})$| is

\begin{equation} (\boldsymbol {q}, \tilde{\boldsymbol {q}}) = \frac{1}{2} \left(\hat{u}_x \hat{\tilde{u}}_x^* + \hat{u}_y \hat{\tilde{u}}_y^* + \hat{W} \hat{\tilde{W}}^*\right), \end{equation}

(13)

where an asterisk denotes complex conjugation and |$\boldsymbol {q} \equiv \lbrace \hat{u}_x, \hat{u}_y,\hat{W}\rbrace$|⁠, |$\tilde{\boldsymbol {q}} \equiv \lbrace \hat{\tilde{u}}_x,\hat{\tilde{u}}_y,\hat{\tilde{W}}\rbrace$|⁠.

In order to find perturbation yielding the highest growth factor, we use the iteration procedure (see ZR14 for details). Each step of the iteration procedure consists of the advance of the perturbation state vector forward in time, employing the dynamical equations

\begin{equation} \frac{\mathrm{\partial} \boldsymbol {q}}{\mathrm{\partial} t} = \boldsymbol {A}\, \boldsymbol {q}, \end{equation}

(14)

and the subsequent integration of adjoint equations

\begin{equation} \frac{\mathrm{\partial} \tilde{\boldsymbol {q}}}{\mathrm{\partial} t} = - \boldsymbol {A}^{\dagger} \, \tilde{\boldsymbol {q}} \end{equation}

(15)

backwards in time. Here, operators |$\boldsymbol {A}$| and |$\boldsymbol {A}^{\dagger}$| are related by the Lagrange rule:

\begin{equation} (\tilde{\boldsymbol {q}},\boldsymbol {A}\boldsymbol {q}) = (\boldsymbol {A}^{\dagger} \tilde{\boldsymbol {q}}, \boldsymbol {q}). \end{equation}

(16)

The explicit form of |$\boldsymbol {A}$| and |$\boldsymbol {A}^{\dagger}$| is defined by the particular model for perturbations and particular choice of their norm. When |$\boldsymbol {A}$| is represented by equations (7)–(9) and the inner product is given by equation (13), operator |$\boldsymbol {A}^{\dagger}$| corresponds to the following set of adjoint equations:

\begin{eqnarray} \frac{{\rm d} \hat{\tilde{u}}_x}{{\rm d} t} &=& (2-q)\hat{\tilde{u}}_y - {\rm i}\tilde{k}_x \hat{\tilde{W}} + R^{-1} k^2 \hat{\tilde{u}}_x + \left(\frac{R^{-1}}{3} + R_{\rm b}^{-1}\right)\nonumber \\ &&\times\, \tilde{k}_x (\tilde{k}_x \hat{\tilde{u}}_x + k_y \hat{\tilde{u}}_y), \end{eqnarray}

(17)

\begin{eqnarray} \frac{{\rm d} \hat{\tilde{u}}_y}{{\rm d} t} &=& -2\hat{\tilde{u}}_x - {\rm i} k_y \hat{\tilde{W}} + R^{-1} k^2 \hat{\tilde{u}}_y + \left(\frac{R^{-1}}{3} + R_{\rm b}^{-1}\right)\nonumber \\ &&\times\, k_y (\tilde{k}_x \hat{\tilde{u}}_x + k_y \hat{\tilde{u}}_y), \end{eqnarray}

(18)

\begin{eqnarray} \frac{{\rm d} \hat{\tilde{W}}}{{\rm d} t} &=& - {\rm i} (\tilde{k}_x \hat{\tilde{u}}_x + k_y \hat{\tilde{u}}_y). \end{eqnarray}

(19)

Note that using equations (48)–(50) and the rule (16) we derive another set of adjoint equations to consider the dynamics of global perturbations with λ ∼ r₀ (see Section 5.3).

The set of equations (17)–(19) looks similar to the set of equations (7)–(9), except that the coefficients in front of |$\hat{u}_{y,x}$| in the inviscid terms of equations (7) and (8), respectively, are flipped over, and the viscous terms acquire opposite signs in the adjoint equations. The latter implies that equations (17)–(18) describe viscous diffusion with a negative diffusion coefficient. However, because the adjoint equations are advanced backwards in time during the implementation of iterative loop, we do not face any sort of instability caused by a non-zero negative viscosity.

In this way, starting from an arbitrary perturbation state vector |$\boldsymbol {q}$| at t = 0, we use the following procedure:

1. integrate equations (7)–(9) forward in time up to the instant t > 0;

2. substitute the resulting |$\boldsymbol {q}(t)$| into equations (17)–(19) and integrate them backwards in time up to the instant t = 0;

3. renormalize the upgraded |$\boldsymbol {q}(0)$| to the unit norm;

4. repeat the previous three steps until |$\boldsymbol {q}(0)$| converges with a required accuracy.

As discussed in detail by ZR14, the final |$\boldsymbol {q}(0)$| represents the optimal SFH for particular values of k_x, k_y, R, R_b and q, which exhibits the highest possible transient growth factor

\begin{equation} g_{\rm opt}(t) \equiv \max _{\forall {\rm SFH}} g(t), \end{equation}

(20)

up to the time t among all other SFHs. We call g_opt the optimal growth hereafter.

Provided that k_y, R, R_b and q are fixed, g_opt has a global maximum over all k_x and t, which we denote as

\begin{equation} G_{\max } \equiv \max _{\forall k_x,\forall t}\, g_{\rm opt}, \end{equation}

(21)

as this quantity is analogous to the maximum optimal growth obtained for 3D incompressible perturbations by MAN05. Note that in Appendix A2 we define an incompressible counterpart of equation (21).

Additionally, we look at an intermediate quantity

\begin{equation} G(t)\equiv \max _{\forall k_x} g_{\rm opt} \end{equation}

(22)

in the case of compressible perturbations.

In what follows, G_max is determined using common numerical routines in the search for the extremum. In the rest of the paper, we present a parametric study of G_max depending on k_y, R and q.² This type of non-modal problem, with the full account of viscous terms in dynamical equations, is solved for the first time. Our goal is to determine accurate constraints on G_max, covering vortical perturbations with λ_y ≳ H along with the case λ_y ≪ H considered previously in the literature. It is also important to check the strong dependence of G_max on q estimated analytically by ZR14 – see also the review by Razdoburdin & Zhuravlev (2015), RZ15 hereafter – which is elucidated in Appendix B of this paper (see equation B3).

4 RESULTS

4.1 Bulk viscosity: suppression of shearing density waves

ZR14 showed that, as far as k_x < 0, the optimal SFH corresponds to a pure vortex, rather than to a shearing density wave or a mixture of both (see also Bodo et al. 2005). This vortex is a leading spiral contracted and amplified by the background shear. In other words, as far as k_x < 0, the optimization procedure for determining g_opt is equivalent to looking for a pure vortex. At the time t_s (which is the instant of swing), the leading spiral turns into a trailing one starting its decay. By definition, |$\tilde{k}_x(t_{\rm s})=0$|⁠. Additionally, at t = t_s the spiral excites a shearing density wave with an initial amplitude defined by k_y. This density wave starts growing because its frequency increases due to the trailing spiral being wound up by the shear. The strongest density waves are generated by vortices with k_y ∼ 1 (see Heinemann & Papaloizou 2009). On the one hand, as soon as k_y ∼ 1, they significantly enhance the optimal growth G(t) as defined by equation (22). On the other hand, it is not clear what role they play in the production and maintenance of turbulence. Additionally, as density waves produce a highly divergent velocity field, they must be strongly affected by the value of bulk viscosity, which is poorly known in astrophysical flows. For these reasons, we focus on the dynamics of vortices only. However, in order to be confident that a reliable quantitative criterion will be chosen to make a conclusion about the transient amplification of vortices irrespective of the excitation of density waves, we should first examine the optimal solutions to the general problem represented by equations (7)–(9).

In this respect, we present the result of the optimization procedure for the case of Keplerian shear in Fig. 1. For the sake of consistency with MAN05 and with analytical estimations shown in Appendix B, here we take R₀₅ = 2000. Note that R₀₅ is related to R by equation (A5). In the top panel of Fig. 1 we can see the highest possible enhancement of linear perturbations in the Keplerian shear, which is plotted versus k_y. First, we start with the case of zero bulk viscosity demonstrated by the solid curve and we find the following.

1. Despite the fact that we consider boundless shear, it does have a maximum at k_y ∼ 1 and wings declining for perturbations with both small and large azimuthal length-scale in comparison with the disc thickness.

2. These wings are matched well by the analytical expressions (B1) and (B3) derived in Appendix B employing the so-called balanced solutions; see, for example, equations (35)–(37) of RZ15.

3. Perturbations with k_y < 1 yield a considerably larger transient growth than previously analysed small-scale perturbations described in the incompressible limit (see Appendix A).

4. In the range 0.3 < k_y < 3, we find a bump that apparently corresponds to the effect of density wave generation by vortices and is referred to as the sonic bump hereafter.

Figure 1.

Top panel: curves of G_max as defined by equation (21) obtained via an iterative loop for R₀₅ = 2000, q = 3/2 and various R_b. The solid curve denotes R_b = ∞. The dashed curves denote (from the top down) R_b = R, 0.1R and 0.01R. The dotted curve shows g(t_s) versus k_y for optimal SFH corresponding to G_max represented by the solid curve. The dot-dashed and dot-dot-dashed curves are produced by equations (B1) and (B3), respectively. Bottom panel: curves of G(t) as defined by equation (22) obtained via an iterative loop for R₀₅ = 2000, q = 3/2 and k_y = 0.3. The solid, dashed, dotted and dot-dashed curves represent R_b = ∞, R, 0.1R and 0.01R, respectively.

Point (4) becomes more clear as we look at the bottom panel of Fig. 1. In order to plot it, we fix k_y = 0.3, which approximately corresponds to the left edge of the sonic bump for the case R_b = ∞ and we obtain the profile of G(t) (see the solid curve). It turns out that, while approaching the area of k_y ∼ 1, G(t) acquires an additional local maximum, which is located farther in time than t_max given by equation (B2). So, t_max ≈ 19 for the parameters taken to obtain the solid curve, which matches the location of its left maximum. If it were not for the prominent excitation of density waves, it would be a unique extremum of G(t) (see fig. 9 of MAN05). Let us denote it by |$G_{\max }^\prime$|⁠. However, as shown by ZR14 in the inviscid model (see their fig. 2), the vortices with k_y ∼ 1, which have already swung well in advance of some fixed time interval t₀, produce an additional bump in the profile of g_opt(k_x, t = t₀) due to their ability to generate sufficiently strong density waves. This additional bump is suppressed by non-zero viscosity at a considerably later time than the local peak located at k_x = −qk_yt₀ and produced by vortices swinging right at t = t₀. As a consequence, G(t) acquires the second local maximum, which we denote by |$G_{\max }^{\prime \prime }$|⁠. Further, as we approach k_y ≈ 1 from either small k_y or large k_y, we enter the domain where |$G_{\max }^{\prime \prime }> G_{\max }^\prime$|⁠, or even |$G_{\max }^{\prime \prime }$| ‘swallows up’ |$G_{\max }^\prime$|⁠. This domain corresponds to the sonic bump in the solid curve in the top panel of Fig. 1, where |$G_{\max } = G_{\max }^{\prime \prime }$|⁠.

Now, let us consider the non-zero bulk viscosity. In both panels of Fig. 1, another three curves of G_max and three curves of G(t) show the transient growth behaviour as we proceed from R_b = ∞ up to R_b = 0.01R. Apparently, such a huge change in the value of bulk viscosity wipes out the density waves and the corresponding sonic bump but does not affect the growth of vortices. In order to see this better, we add the dotted curve in the top panel of Fig. 1, which represents

\begin{equation*} g_{\rm s}\equiv g(t_{\rm s}) \end{equation*}

of the vortical SFH that yields G_max(k_y) for R_b = ∞. Clearly, in the limit R_b → 1 the curve of G_max just approaches the curve of g_s from above. The same is also confirmed by the profiles of |$\hat{u}_{x,y}(t)$| plotted in Fig. 2. These velocity perturbations belong to the optimal vortices producing G_max at k_y = 0.4 in the cases of R_b = ∞ and R_b = 0.001R. We find that the bulk viscosity enhancement leaves |$\hat{u}_{x,y}$| almost unchanged for t < t_s, whereas their prominent oscillations disappear right after the instant of swing at t_s ≈ 16. Also note that in the case of R_b = ∞ the oscillations of |$\hat{u}_x$| start damping approximately at t ≈ 30, which is the time corresponding to the value of optimal growth equal to |$G_{\max }^{\prime \prime }$|⁠. We do not present plots with t_max(k_y), but we do check that t_max, which is numerically determined for wings of G_max(k_y) corresponding to k_y ≪ 1 (k_y ≫ 1), fits the analytical estimation (B2) very well.

Figure 2.

Solid and dashed curves denote |$\hat{u}_x(t)$| for optimal SFH corresponding to G_max(k_y = 0.4) obtained for R₀₅ = 2000, q = 3/2 with R_b = ∞ and R_b = 0.001R, respectively. Dotted and dot-dashed curves are |$\hat{u}_y(t)$| for the same two SFHs.

The conclusion we have just made about the invariance of transient growth of vortical SFH with respect to the magnitude of bulk viscosity can be supported analytically in the case of tightly wound spirals |$|\tilde{k}_x|\gg 1$| (see Appendix C). However, the latter assumption breaks down as soon as the optimal spiral unwinds and enters the swing interval, where |$|\tilde{k}_x|\lesssim 1$|⁠, and one cannot even distinguish between vortices and density waves. Note that the swing interval may be many times longer than |${\sim }\Omega _0^{-1}$| for large-scale vortices. This is why the numerical analysis presented above appears to be necessary to check the robustness of analytical conclusions.

Thus, in the next section, we use g_s corresponding to perturbations exhibiting the maximum optimal growth, G_max, as a measure of the transient growth of shearing vortices.

4.2 Transient growth in super-Keplerian shear of a finite scaleheight

As shown by Butler & Farrell (1992), optimal vortices in an incompressible plane flow are streamwise rolls, which transiently grow into streamwise streaks due to the lift-up effect (see Schmid & Henningson 2001). In the rotating constant specific angular momentum flow, the optimal vortices are, in contrast, the streamwise streaks that generate streamwise rolls, which has been referred to as the anti-lift-up effect; see Antkowiak & Brancher (2007) and the discussion by Rincon et al. (2008). In both cases, the corresponding optimal growth ∝R² and attains several thousands for R₀₅ = 2000 (see MAN05). In other words, this type of transient growth is produced by a nearly axisymmetric (k_y ≈ 0) SFH of 3D perturbations and is limited only by viscous dissipation. However, the growth of rolls in the rotating flow is dramatically suppressed by the onset of epicyclic oscillations. The optimal growth drops steeply as one proceeds to the shear with a positive radial gradient of the specific angular momentum. Indeed, according to the analytical result of MAN05,G_max ≈ 2/(2 − q), which gives, for example, G_max ∼ 20 at q = 1.9. At the same time, the swing amplification of columnar (independent on z) small-scale vortices λ_y ≪ H do not produce a substantial growth before they are damped by the viscosity, so that the corresponding maximum growth has the value G_max ∼ 30 at q = 1.9 according to MAN05, which is only slightly larger than the value obtained for streamwise rolls. Nevertheless, as we have seen in the previous section, the swing amplification of 2D vortical perturbations with λ_y ≳ H proves to be significantly larger than that of vortices with λ ≪ H in the Keplerian case q = 1.5. Additionally, according to the analytical result (B3), the large-scale vortices should attain even larger growth factors as one proceeds to the super-Keplerian shear rate. In other words, the swing amplification of large-scale vortices enhances as the stabilization of fluid elements by Coriolis and tidal forces becomes weaker. For this reason, we present the results of numerical optimization for 1.5 < q < 2.0 involving the compressible dynamics in this section (see Fig. 3). There are solid curves that demonstrate the gradual increase of the growth factor as we approach q = 2. As we see, at q = 1.9, g_s ∼ 10³, which is two orders of magnitude larger than the corresponding incompressible result of MAN05. Also, these curves are matched well by their analytical counterparts, obtained using equation (B3); see the dotted lines in Fig. 3. For confidence, we numerically check that the dependence g_s ∝ R^2/3 is valid throughout the whole range of k_y. Thereby, in contrast to what we see in fig. 4 of MAN05, the transient growth of perturbations in shear flow of a finite scaleheight reduces much more gradually as one introduces a positive radial gradient of the specific angular momentum.

Figure 3.

Solid curves show g_s versus k_y for optimal SFH that attains G_max as defined by equation (21) for R₀₅ = 2000 and R_b = ∞; from bottom to top, q = 1.5, 1.6, 1.7, 1.8, 1.9 and 2.0. Dashed curves show g_kin(t_s) versus k_y for the same SFHs; from bottom to top, q = 1.5, 1.6, 1.7, 1.8, 1.9 and 2.0. Dotted curves represent the analytical estimation (B3) for the same parameters; from bottom to top, q = 1.5, 1.6, 1.7, 1.8 and 1.9.

Besides, g_s diverges as k_y → 0 in the constant specific angular momentum flow, q = 2. In other words, it obeys the same law as in the incompressible limit: |$g_{\rm s}\propto k_y^{-4/3}$|⁠. Using equation (57) from ZR14 at κ = 0, we find that inviscid balanced solutions (35)–(37) of RZ15 yield |$G_{\max } = t_{\max }^2$| in the limit k_y ≪ 1 and |k_x| ≫ 1, which also recovers the incompressible dependence on k_y after we substitute t_max (equation B2). However, the numerical curve for q = 2 displayed in Fig. 3 gives a larger value in the limit k_y ≪ 1, even though we have not included the factor of viscosity damping in the estimate of G_max. Presumably, this discrepancy is a consequence of prominent excitation of density wave up to k_y → 0, which enhances the growth factor of the optimal vortex compared to the value coming from a balanced solution at t = t_s. This claim is supported by the analytical result of Heinemann & Papaloizou (2009), which states that the amplitude of a generated density wave grows exponentially as the parameter ε ≪ 1 increases; see, for example, equation (34) of RZ15. As long as κ ∼ 1 and k_y is not too close to unity, ε is small, the theory constructed by Heinemann & Papaloizou (2009) is valid and the generated wave has a small amplitude at t = t_s. This is why the numerically obtained g_s is close to the estimate coming from the balanced solution for k_y > 1 as well as for k_y < 1. However, as κ → 0, the region of weak wave excitation shifts to smaller k_y (one can compare the pairs of solid and dotted curves for q = 1.5 and q = 1.9 in Fig. 3). For q = 2, the whole region 0 < k_y ≲ 1 suffers strong excitation of density waves, making the use of the balanced solution unjustified.³ However, it still reproduces the correct dependence on k_y.

4.2.1 Kinetic growth factor

Going back to the case q < 2, the question arises for both Keplerian and super-Keplerian shear rates, to what extent either the velocity perturbation or the enthalpy perturbation provide the growth of g_s as one proceeds to large azimuthal scales. To elucidate this, we additionally plot the kinetic part of the growth factor defined as

\begin{equation} g_{\rm kin} = \frac{E_{\rm kin}(t)}{E_{\rm kin}(0)}, \end{equation}

(23)

where |$E_{\rm kin} = 1/2 (|\hat{u}_x|^2 + |\hat{u}_y|^2)$| (see the dashed curves in Fig. 3). As we see, g_ks ≡ g_kin(t_s) of the optimal large-scale vortices ceases much faster than g_s as one shifts to k_y → 0 (i.e. a very large part of g_s is contributed by the amplitude of |$\hat{W}$|⁠). Nevertheless, it is interesting to note that the ratio g_ks/g_s substantially increases as we fix k_y and approach q = 2. This can also be seen in Fig. 4, where we compare the profiles of |$\hat{u}_x$|⁠, |$\hat{u}_y$| and |$\hat{W}$| in the Keplerian and super-Keplerian cases. The amplitude of |$\hat{u}_x$| becomes more than six times larger at the instant of swing, as we change from q = 1.5 to q = 1.8, whereas the amplitude of |$\hat{W}$| increases by a factor of 2 only. However, the kinetic energy is still only a few per cent of the acoustic energy of optimal SFH at the peak of its norm in the particular case introduced in Fig. 4.⁴

Figure 4.

Comparison of |$\hat{u}_{x,y}$| and |$\hat{W}$| for Keplerian and super-Keplerian shear. Dotted, short-dashed and short-dash-dotted curves show |$\hat{u}_x$|⁠, |$\hat{u}_y$| and |$\hat{W}$| versus t for optimal SFH producing G_max(k_y = 0.02) as defined by equation (21) for R₀₅ = 2000, R_b = 0.001R and q = 1.5. Solid, long-dashed and long-dash-dotted curves show |$\hat{u}_x$|⁠, |$\hat{u}_y$| and |$\hat{W}$| versus t for optimal SFH producing G(t) at k_y = 0.02, where q = 1.8, R₀₅ = 2000, R_b = 0.001R and t equals the swing time of the previous SFH.

We are not the first to find the optimal transient perturbations with a significant or even prevailing deposit of enthalpy (alternatively, density or temperature) perturbation amplitude. A number of studies devoted to inspection of compressible transient perturbations in plane shear flows have found this previously. For example, Hanifi et al. (1996), who showed that the lift-up mechanism survives in a moderately supersonic boundary layer, and that optimal perturbations are independent of the streamwise direction similarly to the incompressible case, also found that these optimals contain a large amplitude of temperature perturbation (see their fig. 12). Next, Farrell & Ioannou (2000) solved a 2D optimization problem in a hypersonic plane shear flow, which is more analogous to our model. Similarly to our results, they found that perturbation harmonics in the form of the leading spirals exhibit growth factors much larger than they would be in incompressible fluid because of the excitation of density waves at the instant of swing. For the particular situation, they considered that the initial optimal perturbation mostly consisted of a pressure perturbation (see their fig. 4). This general result is not surprising, at least in the case of 2D dynamics studied by Farrell & Ioannou (2000) and by us, as balanced solutions acquire the non-zero velocity divergence and non-zero density perturbation with an account of compressibility; see equations (35)–(37) of RZ15. It should be kept in mind, however, that the balanced solutions (as well as their accurate numerical counterparts) are indeed vortices because they are proportional to I_ν and are not affected by bulk viscosity (see Section 4.1). Also note that 2D optimal growth in a hypersonic plane flow obtained by Farrell & Ioannou (2000) would be effectively suppressed by high bulk viscosity down to the incompressible value, which is much smaller than the 3D optimal growth produced by the streamwise rolls. However, this is not the case in super-Keplerian rotating flow with q → 2, as we find the enhancement of transient growth of large-scale vortices irrespective of an excitation of density waves (see Fig. 3). Thus, for a super-Keplerian shear rate close enough to q = 2, 2D transient growth produced by the swing amplification of shearing spirals exceeds 3D transient growth of streamwise streaks produced by the anti-lift-up effect. Note that it is correct to compare the growth factors of both types of perturbations as for q = 2 the streamwise-independent solution in the form of the streaks (see equations 52–55 of MAN05) satisfies the set of equations (7)–(9) provided that the equation describing the evolution of the vertical velocity perturbation |$\hat{u}_z$| is added. Also, it can be checked that such a set of equations allows the existence of axisymmetric vortical solutions for q ≲ 2, which are akin to streamwise streaks with a non-modal growth strongly modified by epicyclic motions and slightly modified by divergent motions up to λ ∼ H.

Finally, we argue that it is questionable in what way the high-amplitude density perturbations, emerged from the transient growth of vortices with λ_y ≳ H, should affect the stage of non-linear interaction of different SFHs. Unlike density perturbations of the excited density waves, they can hardly lead to a generation of shocks, because at the peak of their amplitudes (i.e. near the vortex swing) they have the least spatial gradient and, in contrast, as soon as their spatial gradient starts to grow through the spiral winding, their amplitude drops back to smaller values. Contrary to this, the amplitude of W in the density wave grows simultaneously with its spatial gradient, which most likely ends with a shock. It is tempting to suppose that high amplitudes of W generated by large-scale vortices qualitatively change the picture of the non-linear dynamics through the term |$\nabla \rho _1 \boldsymbol {u}$| entering the continuity equation for finite amplitude perturbations. However, to the best of our knowledge, it is not clear whether this has any consequence with regards to the problem of positive non-linear feedback in the bypass transition to turbulence. To date, significant progress has been made on wind tunnel experiments (e.g. Schneider 2006) and on direct numerical simulations of the transition in hypersonic boundary layers; see the review by Zhong & Wang (2012) and one of the recent particular examples by Kudryavtsev & Khotyanovsky (2015). Furthermore, theoretical studies have been carried out with emphasis on the transition scenario based on spectral instabilities of the hypersonic flow; see the review by Fedorov (2011). At the same time, we are not aware of any research analysing possible qualitative changes in the transition scenarios involving the bypass mechanism.

Nevertheless, whatever part of the transient growth factor is contained in enthalpy perturbation, the kinetic part, g_ks, remains significant for optimal SFH with k_y ≲ 1, even in the Keplerian case. No less important, g_ks becomes larger as q → 2, even more steeply than g_s, attaining the value g_kin ∼ 600 for q = 1.9 because of the substantial increase of the |$\hat{u}_x$| amplitude.

4.3 Low limit of transition Reynolds number

Meseguer (2002) found a strong correlation between G_max for 3D incompressible perturbations and the subcritical stability boundary in a viscous counter-rotating spectrally stable Taylor– Couette flow. Recently, Maretzke et al. (2014) suggested that the threshold for a transition in any regime of Taylor–Couette rotation should be G_max ∼ 100, based on a variety of new experimental results and non-modal growth calculations. Of course, it is not clear yet whether this is true for the quasi-Keplerian regime as, to date, nobody has detected its turbulent state. Furthermore, it is not known whether the transition threshold for G_max exists in the case of vortices in the compressible medium considered in this work. Nevertheless, let us consider a conjecture about the value of the transition Reynolds number resulting from the assumption that G_max(R_T) ∼ const, where R_T is the Reynolds number of the subcritical transition in quasi-Keplerian case 1.5 ≤ q < 2. If so, equation (B3) yields

\begin{equation} R_{\rm T} \sim 10^8 \, \frac{(2-q)^3}{q}, \end{equation}

(24)

where we set R_T ∼ 10⁷ for the Keplerian rotation, which is stable at least up to the Reynolds number equal to a few million (see Edlund & Ji 2014).

According to equation (24), the bypass transition to HD turbulence can occur at significantly lower R, for example, in the inner region of relativistic accretion disc. Indeed, using the Paczynski–Wiita potential (Paczynsky & Wiita 1980)

\begin{equation} \Phi _{\rm PW} = - \frac{M}{r-2M}, \end{equation}

(25)

where M is the mass of central object, r is the radial coordinate in the disc plane and it is implied that G = c = 1, we obtain

\begin{equation} q \equiv \frac{r}{\Omega }\frac{{\rm d}\Omega }{{\rm d}r}= \frac{1}{2} + \left(1-\frac{2M}{r}\right)^{-1}. \end{equation}

(26)

Equation (26) is plotted in Fig. 5, which shows that for r < 20M, where the relativistic disc emits most of its energy, q > 1.6 and R_T ≲ 4 × 10⁶.

Figure 5.

Dependence of local shear rate q on radial distance in the Paszynski–Wiita potential (see equation 26).

Next, we note that kinematic viscosity is likely to increase with the radial distance in the standard thin disc model:

\begin{equation} \nu \sim v_{\rm p} l_{\rm p} \propto a_* / N \sim \Omega H^2 / \Sigma \propto r^{1/2+\beta }. \end{equation}

(27)

Here, v_p, l_p and N are, respectively, the mean velocity, the mean free path and the number density of gas particles, and we assume that the disc has a constant aspect ratio H/r and Σ ∝ r^−β with β > 0.⁵ Consequently, the situation is possible when the relativistic disc is turbulent inside some radius r = r_T, whereas turbulence ceases in the outer regions of the disc. Presumably, if this happens, the angular momentum flux is absent throughout the laminar flow, and the inner disc exhibits a non-stationary accretion. It is unclear whether this process might become recurrent, and this requires additional research.

Furthermore, it is known that in a disc with a super-Eddington accretion rate the radial pressure gradient becomes significant, flattening the specific angular momentum distribution of the rotating matter (see Abramowicz et al. 1988). Thereby, the actual q(r) should become even steeper compared to the curve shown in Fig. 5. Of course, the enhanced transient growth of large-scale vortices may drive turbulization and associated angular momentum outflow in the Newtonian discs, where the super-Keplerian rotation may arise solely due to the radial pressure gradient. More generally, the aforementioned mechanism may be at work whenever the super-Keplerian rotation takes place on streamwise scales comparable to or exceeding the vertical scaleheight in a flow. Thus, in principle, a kind of a secondary bypass transition can occur inside the elongated non-linear vortices, which produce a sufficient excess of vorticity on the smooth background shear they are embedded in. Those non-linear vortices have been found in a large number of simulations including those mentioned in the introductory section (Umurhan & Regev 2004; Lithwick 2009). If elongated in the streamwise direction, they are known to be almost linearly stable (see Lesur & Papaloizou 2009), and are efficient in trapping the dust particles in protoplanetary discs (e.g. Heng & Kenyon 2010).

Another consequence of the dependence of G_max on q takes place for the dimensionless azimuthal viscous stress α if it is caused by turbulent motions fed by transiently growing large-scale shearing vortices. Lesur & Longaretti (2005) suggested that |$\alpha \propto R_{\rm T}^{-1}$|⁠, and thus α should change with q as

\begin{equation} \alpha /\alpha _K \sim \frac{q}{12 (2-q)^3}, \end{equation}

(28)

where α = α_K is its Keplerian value. Local and global MHD simulations of relativistic accretion discs (e.g. Penna et al. 2013; Morales Teixeira et al. 2014) show a radially varying α decreasing to the periphery of disc, though this is what one finds in a supercritical turbulence fed by MRI. In Fig. 6, we compare equation (28) with the expressions given in equation (9) of Abramowicz (2009) and in the abstract of Penna et al. (2013). Because q increases as we approach the inner edge of the disc (see equation 26), we suppose that if the subcritical turbulence based on the bypass mechanism and fed by transiently growing quasi-2D optimal shearing vortices were resolved in simulations, one would obtain a steeper dependence of α on r compared to what has been obtained to date with regards to supercritical magnetic turbulence. In the case of a relativistic disc, such a substantial radial variation of α on r would result in a significant drop in temperature of the inner parts of disc for a fixed accretion rate |$\dot{M}$|⁠, thus increasing the limiting value of |$\dot{M}$| corresponding to the Eddington luminosity in the standard disc model. According to equation (28), α diverges in the vicinity of the inner edge of a standard relativistic disc. However, the finite aspect ratio of the disc, δ = H/r₀, puts a lower limit on the value of k_y ∼ δ, which constrains the value of optimal growth as q → 2 (see Fig. 3). Additionally, the non-zero pressure radial gradient close to the inner edge of a standard disc always results in the shear rate q < 2.

Figure 6.

Dependence of α on local shear rate q (see Section 4.3). The solid curve represents equation (28). The dashed and dotted curves correspond to the results of Abramowicz, Brandenburg & Lasota (1996) and Penna et al. (2013), respectively.

5 EXTENSION TO GLOBAL SPATIAL SCALE

So far, we have been considering the transient dynamics of large-scale shearing vortices in the local framework implying that r₀ ≫ λ_y ≳ h. However, as k_y ∼ δ, we have to take into account a global disc structure and the cylindrical geometry of the problem. At first, the disc generally is radially non-uniform as its surface density and its microscopic viscosity vary on the scale ∼r₀. Next, the transient amplification of large-scale vortices substantially increases along with the local shear rate, which in turn might be a function of the radial distance. For example, in a relativistic disc model with the Paczynsky–Wiita gravitational potential, q(r) is given by equation (26). Additionally, even if q is a constant throughout the disc, which is the case for the Newtonian potential, global perturbations may be affected by a non-zero radial gradient of the vorticity, which vanishes in the local spatial limit.

We have found in this work that the above factors can considerably alter a value of G_max for global shearing vortices or even make the flow unstable. In order to see this, we present a supplementary study of global optimal perturbations using a set of global disc models with two variants of surface density and microscopic viscosity radial distributions in both the Newtonian and the Paczynsky–Wiita gravitational potentials. At the same time, we finally come to the overall conclusion that, regardless of the particular radial structure of a disc, it is always able to generate prominent transient growth of global perturbations at huge Reynolds numbers realized in astrophysical conditions.

5.1 Background models

For the first variant of the background, we would like to take profiles of surface density, disc semi-thickness and disc radial velocity, which correspond to a stationary accretion with a radially constant |$\dot{M}$| occurring due to the fluid microscopic viscosity.

5.1.1 Stationary configuration

A geometrically thin disc in the Paczynsky–Wiita potential (25) rotates with the angular frequency

\begin{equation} \Omega = \frac{M^{1/2}}{r^{1/2}\left(r-2M\right)}, \end{equation}

(29)

and is subject to radial oscillations with the following epicyclic frequency

\begin{equation} \kappa = \frac{M^{1/2}}{\left(r-2M\right)^{3/2}}\left(1 - \frac{6M}{r}\right)^{1/2}. \end{equation}

(30)

Another quantity entering equations for perturbations is the sound speed at the equatorial plane of the disc, a_eq. Provided that ρ(z = H) = 0 and that the disc is polytropic along the vertical direction, there is the following exact relation between the disc thickness and the equatorial speed of sound:

\begin{equation} a_{\rm eq}=\frac{H\Omega }{(2n)^{1/2}}. \end{equation}

Because, in our 2D problem, we are dealing with vertically integrated equations for perturbations, which contain |$a_*^2=n a_{\rm eq}^2/(n+1/2)$| (see ZR14), we obtain

\begin{equation} a_* = \frac{H\Omega }{(2n+1)^{1/2}}. \end{equation}

(31)

Next, we have to specify the kinematic viscosity coefficients, ν and ν_b, across the flow. Similarly to what has been defined in Section 2, these are the shear and the bulk viscosities, η and ζ, integrated over the disc thickness and subsequently divided by surface mass density, Σ. As we have mentioned above (see equation 27), it is reasonable to assume that ν ∝ H²Ω/Σ. Provided that the value of ν is known at some distance r = r₁,

\begin{equation} \nu (r)=\nu _1 \frac{\Sigma _1}{H_1^2\Omega _1} \frac{H^2(r)\Omega (r)}{\Sigma (r)}, \end{equation}

(32)

where the variables marked by the subscript ‘1’ stand for the corresponding quantities taken at r₁. Let us choose r₁, which is the distance where ν(r) attains the minimum. A value of ν₁ can be specified via the Reynolds number R (cf. equation 10):

\begin{equation} R = \frac{a_{*1}^2}{\Omega _1 \nu _1} = \frac{H_1^2 \Omega _1}{\nu _1 (2n+1)}. \end{equation}

(33)

Combining equations (32) and (33), we obtain the final expression for the kinematic viscosity coefficient profile,

\begin{equation} \nu = \frac{q \Sigma _1}{(2n+1)R_{05}}\frac{H^2(r)\Omega (r)}{\Sigma (r)}, \end{equation}

(34)

where it is taken into account that R₀₅ = q(r_i)R (see Section A1) and r_i is the inner boundary of disc.

With regards to the second viscosity, we assume that the second (local) Reynolds number is constant throughout the disc, which implies that

\begin{equation} \nu _{\rm b}= \frac{a_{*}^2}{\Omega R_{\rm b}}. \end{equation}

(35)

In the next stage, we obtain an expression for the background radial velocity, v_r, and the disc half-thickness H. For this purpose we use the azimuthal projection of the Navier–Stokes equations integrated over z

\begin{equation} \Sigma v_r \frac{{\rm d}}{{\rm d}r} (\Omega r) + \Sigma v_r \Omega =\frac{1}{r^2}\frac{{\rm d}}{{\rm d} r}\left(\Sigma \nu r^3 \Omega ^{\prime }\right) \end{equation}

(36)

and mass conservation in the form

\begin{equation} v_r=-\frac{\dot{M}}{2\pi \Sigma r}, \end{equation}

(37)

with |$\dot{M}$| constant throughout the disc.

Equation (37) allows us to integrate equation (36) assuming that the rφ-component of the stress tensor vanishes at the inner boundary of the disc r = r_i:

\begin{equation} \dot{M} \Omega r^2 D=-2\pi \Sigma \nu r^3 \Omega ^{\prime }. \end{equation}

(38)

In equation (38), we introduce |$D=1-(\Omega _{\rm i} r_{\rm i}^2)/ (\Omega r^2)$|⁠. Note that for the Newtonian potential D = 1 − (r_i/r)^1/2, whereas for the Paczynsky–Wiita potential D = 1 − (r_i/r)^3/2(r − 2M)/(r_i − 2M).

Combining equations (37) and (38), we obtain the radial velocity

\begin{equation} v_r=\frac{\nu \Omega ^{\prime }}{\Omega D}. \end{equation}

(39)

Then, combining equations (32) and (38), we are able to obtain the radial profile of H:

\begin{equation*} H\propto \left(-\frac{D}{r\Omega ^{\prime }}\right)^{1/2}. \end{equation*}

In order to parametrize H, we use the disc aspect ratio |$\delta \equiv (H/r)|_{r=r_2}$|⁠, where r₂ is the distance at which H/r takes its maximum value. Finally,

\begin{equation} H=\delta r_2\left(\frac{r_2 \Omega _2^{\prime }}{r \Omega ^{\prime }}\right)^{1/2}\left(\frac{D}{D_2}\right)^{1/2}, \end{equation}

(40)

where the variables marked by the subscript ‘2’ denote the corresponding quantities taken at r₂.

At last, in order to specify the radial profile of surface density for such a disc we proceed in a similar way as Shakura & Sunyaev (1973) and assume that viscous dissipation is locally balanced by the radiation flux from the disc surface:

\begin{equation} Q=\frac{1}{2}\Sigma \nu r^2 \left(\Omega ^{\prime }\right)^2. \end{equation}

(41)

Next, the diffusive radiation transfer in the vertical direction yields an approximate relation between the surface and the equatorial temperatures, T_eff and T_eq:

\begin{equation*} T_{\rm eq}^4 / T_{\rm eff}^4 \propto \tau = H / l_{\rm ph}. \end{equation*}

Here, l_ph ∝ ρ⁻¹ is the mean free path of photons in the case when the electron scattering prevails the free–free absorption opacity. Thus, assuming that the non-relativistic gas pressure prevails, we conclude that

\begin{equation} T_{\rm eq}^4 \propto a_{\rm eq}^8 \propto Q \Sigma . \end{equation}

(42)

Combining equations (34), (40), (41) and (42), we finally obtain

\begin{equation} \Sigma =C_\Sigma \frac{D^3\Omega ^7}{r^5 (-\Omega ^{\prime })^5}, \end{equation}

(43)

where C_Σ is a dimensional constant. Because the problem for dynamics of perturbations is linear, an explicit form of |$C_\Sigma$| is unimportant for us.

Equations (29), (30), (31), (34), (35), (39), (40) and (43) specify our stationary background configuration. Hereafter, we denote this as ‘P1’ in the case of Ω and κ given by equations (29) and (30); otherwise (i.e. in the limit r ≫ 1 when Ω and κ acquire their Newtonian values) we denote the background as ‘N1’.

5.1.2 Homogeneous configuration

The models N1 and P1 describe the quasi-stationary stage of accretion that can be achieved by the inner region of an arbitrary rotating laminar flow with δ ≪ 1 evolving due to microscopic viscous relaxation only. However, some would argue that this stage is unrealistic, as microscopic viscosity is too small in typical astrophysical situations leading to a very long relaxation time, which most likely exceeds the turbulization time of any initial rotating configuration by orders of magnitude. For this reason, we regard the models N1 and P1 just as the limiting case, which we would like to counterpose to a simplified radially homogeneous configuration described below.

As before, the frequencies and the sound speed are given by equations (29)–(31). However, now the surface density is assumed to be a constant

\begin{equation} \Sigma = {\rm const}. \end{equation}

(44)

Accordingly, we assume that the kinematic viscosity coefficient is given by (cf. equation 34)

\begin{equation} \nu = \frac{a_{*}^2}{\Omega R}. \end{equation}

(45)

The bulk viscosity coefficient is given by equation (35) as before. The half-thickness profile is chosen to be

\begin{equation} H=\delta r. \end{equation}

(46)

Note that in a disc of such a simple shape the locally defined dimensionless azimuthal wavenumber at every distance corresponds to a global harmonic of perturbations with one and the same azimuthal wavenumber (see equation 58). At last, the radial velocity is set to zero:

\begin{equation} v_r=0. \end{equation}

(47)

Thereby, equations (29)–(31), (35) and (44)–(47) specify our homogeneous configuration. Similarly, hereafter we denote this configuration as ‘P2’ in the case of Ω and κ corresponding to the Paczynsky–Wiita potential, whereas for the Newtonian case the model is denoted as ‘N2’.

5.2 Equations for perturbations

In order to derive the set of equations for global perturbations, we linearize the Navier–Stokes equations written in the cylindrical coordinates {r, φ, z}; see, for example, Appendix B of Kato et al. (2008). Following Section 2, we do not consider the energy balance in viscous fluid adopting the equation of state in the perturbed flow to have the same polytropic index as in the background equation of state. Similarly to what has been done in the local case, we restrict ourselves by considering the planar velocity perturbations independent of vertical coordinate. Thus, perturbations are described by three variables: the Eulerian perturbations of velocity components, δv_r and δv_φ, and the Eulerian perturbation of enthalpy δh. Because of the rotational symmetry of background, it is sufficient to consider complex amplitudes of the Fourier harmonic ∝e^imφ, which depend on r and t only. Integrating equations over z, we obtain the non-viscous terms to be exactly the same as in section 3.1 of ZR14. The viscous part of the equations contains terms emerging due to the non-zero perturbation of shear (which we denote by N_{r, φ}) and terms emerging due to the non-zero perturbation of the velocity divergence (which we denote by B_{r, φ}); additionally, the advective terms ∝v_r must be retained on a global scale. Note that the advective terms do not affect the dynamics of local perturbations (see comments after equations 4 and 5). Explicitly,

\begin{eqnarray} \frac{\mathrm{\partial} \delta v_r}{\mathrm{\partial} t}+{\rm i}m\Omega \delta v_r&=&2\Omega \delta v_{\varphi } -\frac{\mathrm{\partial} \delta h}{\mathrm{\partial} r}-\frac{\mathrm{\partial} }{\mathrm{\partial} r}\left(v_r \delta v_r\right)\nonumber \\ &&+\,N_r+B_r, \end{eqnarray}

(48)

\begin{eqnarray} \frac{\mathrm{\partial} \delta v_{\varphi }}{\mathrm{\partial} t}+{\rm i}m\Omega \delta v_{\varphi }&=&-\frac{\kappa ^2}{2\Omega }\delta v_r - \frac{im\delta h}{r}-\frac{v_r}{r}\frac{\mathrm{\partial} }{\mathrm{\partial} r}\left(r \delta v_{\varphi }\right)\nonumber \\ &&+\,N_{\varphi }+B_{\varphi }, \end{eqnarray}

(49)

\begin{eqnarray} \frac{\mathrm{\partial} \delta h}{\mathrm{\partial} t}+{\rm i}m\Omega \delta h&=&-\frac{a_*^2}{\Sigma r}\frac{\mathrm{\partial} }{\mathrm{\partial} r}\left[r\Sigma \left(\delta v_r+\frac{v_r \delta h}{a_*^2}\right)\right]\nonumber \\ &&-\,\frac{{\rm i}m a_*^2}{r}\delta v_{\varphi }, \end{eqnarray}

(50)

where the viscous terms N_i and B_i are shown in Appendix D1. In equations (48)–(50) we have also assumed that the Eulerian perturbations of dynamical and bulk viscosities are equal to zero δη = δζ = 0. This is done to exclude viscous overstability of perturbations leading to their modal growth irrelevant to (both modal and non-modal) amplification of non-axisymmetric perturbations (see, e.g. Kato 2001).

5.3 Adjoint equations

In order to obtain optimal perturbations via the iterative loop, we have to derive the set of adjoint equations (i.e. an explicit form of the adjoint operator |$\boldsymbol {A}^{\dagger}$| is needed). In the global problem, |$\boldsymbol {A}^{\dagger}$| is defined by the Lagrange rule (see equation 16) through the inner product

\begin{equation} (\boldsymbol {q},\tilde{\boldsymbol {q}})=\pi \int \Sigma \left(\delta v_{r} \delta \tilde{v}_{r}^{*} + \delta v_{\varphi } \delta \tilde{v}_{\varphi }^{*} + \frac{\delta h \delta \tilde{h}^{*}}{a_*^2}\right)r\,{\rm d}r, \end{equation}

(51)

corresponding to a norm equal to total acoustic energy of perturbations

\begin{equation} ||\boldsymbol {q}||^2=E_a=\pi \int \Sigma \left(|\delta v_r|^2+|\delta v_{\varphi }|^2+\frac{|\delta h|^2|}{a_*^2}\right)r\,{\rm d}r. \end{equation}

(52)

Derivation of the adjoint equations is described in detail in section 3.4 of RZ15 for the inviscid case. However, it can be straightforwardly generalized on to viscous problem. We obtain

\begin{eqnarray} \frac{\mathrm{\partial} \delta \tilde{v}_r}{\mathrm{\partial} t}+{\rm i}m\Omega \delta \tilde{v}_r&=&\frac{\kappa ^2}{2\Omega }\delta \tilde{v}_{\varphi }-\frac{\mathrm{\partial} \delta \tilde{h}}{\mathrm{\partial} r}-\frac{v_r}{\Sigma r}\frac{\mathrm{\partial} }{\mathrm{\partial} r}\left(\Sigma r \delta \tilde{v}_r\right)\nonumber \\ &&+\,\tilde{N}_r+\tilde{B}_r, \end{eqnarray}

(53)

\begin{eqnarray} \frac{\mathrm{\partial} \delta \tilde{v}_{\varphi }}{\mathrm{\partial} t}+{\rm i}m\Omega \delta \tilde{v}_{\varphi }&=&-2\Omega \delta \tilde{v}_r-\frac{{\rm i}m\delta \tilde{h}}{r}-\frac{1}{\Sigma }\frac{\mathrm{\partial} }{\mathrm{\partial} r}\left(v_r \Sigma \delta \tilde{v}_{\varphi }\right)\nonumber \\ &&+\,\tilde{N}_{\varphi }+\tilde{B}_{\varphi }, \end{eqnarray}

(54)

\begin{eqnarray} \frac{\mathrm{\partial} \delta \tilde{h}}{\mathrm{\partial} t}+{\rm i}m\Omega \delta \tilde{h}&=&-\frac{a_*^2}{r\Sigma }\frac{\mathrm{\partial} }{\mathrm{\partial} r}\left(r\Sigma \delta \tilde{v}_r\right)-\frac{{\rm i}m a_*^2}{r}\delta \tilde{v}_{\varphi }\nonumber \\ &&-\,v_r\frac{\mathrm{\partial} \delta \tilde{h}}{\mathrm{\partial} r}+\tilde{N}_h, \end{eqnarray}

(55)

where the viscous terms |$\tilde{N}_{\rm i}$| and |$\tilde{B}_{\rm i}$| are shown in Appendix D1.

In order to solve the dynamical and the adjoint equations, we formulate the boundary conditions for viscous perturbations (see Appendix D2).

5.4 Numerical method

For the numerical advance of perturbations obeying equations (48)–(50) and (53)–(55), we use the so-called leap-frog method (see, e.g. Frank & Robertson 1988). This second-order method employs four uniform grids in the space (r, t) shifted to each other, which allows us to use central differences for the approximation of spatial derivatives of the unknown variables. However, the approximation of the advection terms (that ∝v_r) by central differences leads to numerical instability. In order to check this, one should apply Von Neumann stability analysis method (see Charney, Fjörtoft & von Neumann 1950). This is why we use forward differences to approximate those terms in equations (48)–(50), whereas in equations (53)–(55) advanced back in time the advection terms are approximated by backward differences (a description of finite differences of various types can be found in Olver 2014).

Also note that in order to evaluate |$\mathrm{\partial} \delta \tilde{h} / \mathrm{\partial} t$| in the grid node neighbouring the boundary, the value of |$\delta \tilde{h}$| at the boundary is required because we use backward differences in the scheme for the adjoint equations. On the contrary, in the basic equations, the forward differences are used, so the boundary condition for δh is not required.

The allocation of the variables in the grids is the same as in the inviscid problem (see fig. 10 of RZ15).

5.5 Results

As a result of the global optimization problem, we obtain curves of optimal growth

\begin{equation} G(t) \equiv \max \frac{E_a(t)}{E_a(0)}, \end{equation}

(56)

where E_a is given by equation (52). Also, it is common practice to consider another quantity

\begin{equation} G_{\max } \equiv \max _{\forall t} G(t), \end{equation}

(57)

which is similar to G_max introduced by MAN05 as well as Yecko (2004) and Maretzke et al. (2014) (also see equation 21 as its local analogue). Note that G_max depends on the parameters of the disc such as δ, R₀₅ and R_b and on the azimuthal wavenumber m of perturbations. In this section, we present the results related to the first azimuthal harmonic of global perturbations, m = 1, with polytropic index n = 3/2. Also, to exclude the global analogue of density wave excitation during the swing of a global spiral of perturbations (see Section 4.1 for details about the corresponding effect in the local framework), we use a high value of bulk viscosity R_b = 0.001R₀₅ in all calculations. It is implied that time is expressed in units of Ω⁻¹(r_i), the frequencies are expressed in units of Ω(r_i) and radial distance is expressed in units of r_i.

To be able to make comparisons with the corresponding local results on transient growth of large-scale vortices, we additionally introduce a global analogue of k_y

\begin{equation} \bar{k}_y=\frac{m H}{(2n+1)^{1/2} \, r} \end{equation}

(58)

and a global analogue of R

\begin{equation} \bar{R}=\frac{a_*^2}{\Omega \nu } \end{equation}

(59)

(cf. equation 10). In general, |$\bar{k}_y$| and |$\bar{R}$| are functions of r. However, in the models N2 and P2 they are constant throughout the disc. Obviously, the transient growth of the respective local perturbations is determined by values of those quantities at a certain radial distance. Also, it makes sense to give an interpretation of some features of global transient growth with the help of equations (58) and (59).

5.5.1 Models N1 and N2

Let us consider the optimal growth of global perturbations in the case of a Newtonian background. We start by showing the behaviour of G as a function of time in the model N2 (see Fig. 7). At first glance, the optimal growth for m = 1 still attains maximum, which is expected given the results of local analysis (see the bottom panel of Fig. 1), though for some reason G tends to a non-zero value as t → ∞. However, we additionally show the curve of G(t) describing the behaviour of perturbations with m = 10, which has to reproduce better the known bell-shape of the local optimal growth curve obtained previously, for example, by MAN05. By contrast, it can be seen that as m → ∞ the hump in G(t) tends to disappear and the optimal growth almost monotonically approaches its plateau value across the considered time. A clue to this discrepancy comes from checking profiles of optimal initial perturbations (see the other two panels in Fig. 7). Clearly, the phase of growth of G corresponds to the upwind of global leading spiral of initial perturbations. Using terms of the local approach, we say that the upwind continues until the time of order of t_max (see equation B2), evaluated near the inner boundary of the disc. Equation (B2) gives t_max ≈ 25 for the parameters used in Fig. 7, which is, indeed, close to what we find for perturbations with m = 10. As the characteristic time of viscous dissipation for spiral swinging at t becomes shorter than t (i.e. as t > t_max) its transient growth decreases, which is a single outcome for local shearing harmonics as well as in the case of perturbations confined between the walls considered in previous studies. However, in the radially infinite flow similar to N2 the optimal spiral starts shifting outwards instead of further upwind. In terms of local approach, this means that the instant of swing stays approximately equal to t_max taken in units of Ω⁻¹ at the new r = r_sp where the spiral has shifted. The time interval until the instant of swing, t_maxΩ⁻¹(r_sp), becomes longer and longer and it is always approximately equal to the optimization time-span, t. In fact, the latter is a general condition for the approximate localization of global spiral:

\begin{equation} t_{\max }[\bar{R}(r_{\rm sp}), \bar{k}_y(r_{\rm sp}),q(r_{\rm sp})] \Omega ^{-1}(r_{\rm sp}) \approx t. \end{equation}

(60)

At the same time, the global optimal spiral with m ≫ 1 corresponding to |$\bar{k}_y$| should recover the local value of G in whatever part of the disc it was localized. However, from the above, it follows that this local value of G must be equal to nothing but G_max. Then, as |$\bar{k}_y$| and |$\bar{R}$| (as well as q) are constant in N2, the locally determined t_max and G_max are also constant and, for m = 10, we do see the plateau of G(t) at t ≳ t_maxΩ⁻¹(r_i), which is close to the analytical value of G_max given by equation (B3) with the values of |$\bar{k}_y$| and |$\bar{R}$| substituted there (see the dot-dashed curve in the top-left panel of Fig. 7).

Figure 7.

Optimal transient growth of global perturbations with m = 1 in the model N2. Curves are presented for parameters δ = 0.4, n = 1.5, R₀₅ = 2000 and R_b = 0.001R. The top-left panel (solid curve) shows the optimal growth as a function of time, whereas the top-right and bottom panels show the profiles of ℜ [δh] expressed in arbitrary units versus the radial coordinate. These profiles correspond to the optimal initial perturbations. The solid and dashed curves in the top-right panel correspond to the optimization time-spans denoted, respectively, by the cross and the square in the top-left panel. The solid, dashed and dot-dashed curves in the bottom panel correspond to the optimization time-spans denoted, respectively, by the circle, triangle and diamond in the top-left panel. Additionally, the dotted curve in the top-left panel shows the optimal growth for m = 10 and δ = 0.04, whereas the dot-dashed curve is given by analytical estimation (B3) with |$k_y=\bar{k}_y$| and |$R=\bar{R}$| corresponding to the set of parameters used on this plot.

The picture described above indicates that the hump in the curve of G(t) and the related suppress of plateau for perturbations m ∼ 1 in comparison with local value of G are purely global features. Actually, this extra suppression of optimal growth takes place because of an additional spiral upwind after G has already reached its maximum value corresponding to the equality of viscous dissipation time and optimization time (compare the initial optimal profiles corresponding to the square and circle in the top-left panel of Fig. 7). As has been discussed, such an additional upwind vanishes when m → ∞. So, this purely global effect can probably be ascribed to a non-zero vorticity gradient as the most simple variant of background configuration corresponding to |$\bar{k}_y=const$| and |$\bar{R}\,=\,const$| has been considered herein. However, a clear explanation remains to be given. As Fig. 7 illustrates, the optimal spirals with m = 1 also upwind/shift in good agreement with the general picture that has been suggested and expressed through equation (60): for example, equation (60) yields r_sp ≈ 3 and r_sp ≈ 4, respectively, for the dashed and dot-dashed curves in the bottom panel of Fig. 7.

Another purely global effect that is absent in the local formulation of the problem is modal growth of perturbations with m ∼ 1 in the range of sufficiently small δ ≲ 0.05. In Fig. 8, the dependence G(T) and corresponding radial profiles of optimal initial perturbations are shown for δ = 0.02. In contrast to the previous case, a sufficiently thin disc represented by model N2 is unstable with respect to linear perturbations, which is indicated by G exponentially growing at t → ∞. Obviously, at large t the optimal growth recovers the growth factor of the most unstable global mode. The amplitude of this mode increases ∝exp (γt), with γ ≈ 0.001 in Fig. 8.

Figure 8.

Same as in Fig. 7 but for δ = 0.02. The curves in the top-right and bottom panels correspond to the optimization time-spans marked by symbols in the top-left panel in the same way as in Fig. 7.

Along with this fact, at small t the optimal growth G represents the value of transient growth (i.e. the transient growth dominates the modal growth). This conclusion is drawn because the corresponding profiles of optimal initial perturbations are nothing but leading spirals (see the top-right panel in Fig. 8). Also, it can be checked that the growth factor of those spirals damps at large t. Further, at longer t (see the bottom panel of Fig. 8), the optimal initial profile has nothing in common with leading spirals and converges to the unique profile, which is the shape of the most unstable mode. We have checked that it corresponds to a pattern rigidly rotating with angular velocity ω_m ≈ 0.98 (i.e. with corotation radius located inside the flow). It is beyond the scope of this work to study the physical nature of the revealed instability, though presumably it is akin to the well-known Papaloizou–Pringle instability (Papaloizou & Pringle 1987), which has also been studied in detail in 2D semi-infinite rotating flows by Glatzel (1987a) and Glatzel (1987b). Taking into account the results obtained in those papers, we suppose that our unstable mode is a global non-axisymmetric sonic mode arising due to the existence of the boundary and supplied with energy by the background flow at the critical layer localized at the corotation radius. The latter mechanism of instability is usually called the Landau mechanism; see the book by Fabrikant, Stepanyants & Stepaniants (1998) for details. Our suggestion is confirmed by the fact that instability rapidly disappears as we go to perturbations with m > 1, whereas the energy exchange between the mode and the flow is possible only in presence of the non-zero vorticity gradient, which primarily affects perturbations on the largest spatial scale.

Next, we proceed with model N1 keeping the same values of all parameters of the problem (see Fig. 9). Clearly, the shape of G(t) has substantially changed compared to model N2. The plateau has disappeared and instead we obtain a curve with a global maximum and G → 0 as t → ∞ for perturbations with both m = 1 and m = 10. Yet, at the time after G attains maximum, this variant of G(t) does not resemble the bell-shaped curve known previously (see, e.g. MAN05). Again, we see a complex behaviour of the initial optimal profile of perturbations combined of the spiral upwind and the spiral shift along the radial direction. The curves in the bottom panel of Fig. 9 demonstrate that as the optimization time-span increases, the initial optimal spiral shifts outwards and unwinds. This occurs because for r > r₁ (see equation 32) the viscosity grows to the periphery of disc, so that |$\bar{R}$| decreases to the periphery of disc for r > r₁. Consequently, the dimensionless t_max expressed in units of Ω⁻¹(r) decreases outwards as well. Furthermore, as we have suggested above, the optimization time is equal to the instant of swing of the optimal spiral, which, in turn, is approximately equal to t_max. Because the degree of spiral winding (i.e. the initial ratio of radial and azimuthal wavenumbers in terms of the local approach) is proportional to the dimensionless t_max expressed in units of Ω⁻¹(r), we see that the optimal spiral must unwind as it shifts to the periphery of disc specified by the model N1. Note that t_max expressed in units of Ω⁻¹(r_i), which is given by equation (60), remains a quantity that grows outwards; otherwise, the optimal spiral would not be shifted at larger r_sp for larger optimization time-spans.

Figure 9.

Same as in Fig. 7 but for the model N1. The curves in the top-right and bottom panels correspond to the optimization time-spans marked by symbols in the top-left panel in the same way as in Fig. 7. The dot-dashed curve in the top-left panel has been obtained using the local analytical estimations (B2) and (B3) (see details in text).

In order to produce the analytical G(t) corresponding to model N1, we fix an arbitrary t and first estimate a location of global optimal spiral, r_sp, using the condition (60). After that, we find the respective value of |$G_{\max }[\bar{k}_y(r_{\rm sp}), \bar{R}(r_{\rm sp}), q(r_{\rm sp})]$|⁠. This procedure allows us to obtain the dot-dashed curve in the top-left panel of Fig. 9 using the definitions (58) and (59) together with specifications of the model given in Section 5.1.1. As we see, the analytical curve obtained by means of simple estimations of local optimal growth is in reasonable agreement with the numerical curve obtained by employing the iterative loop for global perturbations with m = 10. Once again, even the global spirals with m = 1 upwind/unwind and shift along the radial direction in correspondence with interpretation produced in terms of local approach: for example, equation (60) yields r_sp ≈ 7 and r_sp ≈ 12, respectively, for the dashed and dot-dashed curves in the bottom panel of Fig. 9.

We have numerically checked that the modal growth of perturbations with m = 1 does not exist for any value of δ in the case of model N1. Presumably, the unstable mode, which is located near the inner boundary of the disc (see Fig. 8), is suppressed by a sufficiently large viscosity therein, because in the model N1 |$\bar{R} \rightarrow 0$| as r → r_i.

To complete results related to the flow in the Newtonian potential, we compare the local dependence G_max(k_y) with its global analogue (see Fig. 10). In the case of global perturbations with m = 1, G_max can be defined as the height of the hump in the curve of G(t). As can be seen in Fig. 10, both global and local large-scale vortices with |$\bar{k}_y <1$| produce the transient growth of approximately the same value in discs with different aspect ratios and also in discs specified by different models. Moreover, the disc of the uniform surface density specified by model N2 gives rise to transiently growing spirals with m = 1, which are able to exceed the local value of optimal growth as soon as |$\bar{k}_y \lesssim 0.1$|⁠. At the same time, for larger values of |$\bar{k}_y$| corresponding to geometrically thicker discs, as one considers the case of m = 1 only, the model N2 demonstrates that the global spirals exhibit a substantially weaker transient growth in comparison with the local value of G. The latter is in agreement with the results obtained in ZR14, where in the inviscid problem it is revealed that the optimal growth of large-scale vortices (i.e. λ_y > h) stays almost at the same level as one passes from local perturbations to global perturbations with m ∼ 1, whereas small-scale vortices (i.e. λ_y < h) with m ∼ 1 reach substantially smaller growth factors than their local counterparts.

Figure 10.

Dependence of G_max on |$\bar{k}_y\equiv m\delta /(2n+1)^{1/2}$| for n = 1.5, m = 1, R₀₅ = 2000 and R_b = 0.001R. The solid and dashed curves represent models N1 and N2, respectively, whereas the dot-dashed curve represents G_max(k_y) for local perturbations with q = 3/2, R₀₅ = 2000 and R_b = 0.001R. The termination of the dashed curve corresponds to the occurrence of modal growth.

At last, let us make an important note about the coexistence of modal and transient growth of perturbations with m ∼ 1 in the sufficiently thin discs specified by the model N2. In Figs 7–10, we fix R₀₅ = 2000 as is done in the local problem for easier comparison with the results of MAN05. However, one should keep in mind that the actual Reynolds numbers in astrophysical discs are orders of magnitude larger and approach values such as ∼10¹⁰–10¹³. This implies that, in practice, the transient growth dominates the modal growth by many orders of magnitude for any value of δ, as the optimal growth |$\propto R_{05}^{2/3}$|⁠, whereas the increment increases by a few per cent only as R₀₅ → ∞.

5.5.2 Models P1 and P2

Let us see how the optimal growth changes as we consider the Paczynski–Wiita gravitational potential. In this case, the epicyclic frequency (30) vanishes at the inner boundary, κ(r_i) = 0 (i.e. r_i is the last stable circular orbit in the vicinity of gravitating centre). According to the results obtained in the local approach, the optimal growth of large-scale perturbations significantly enhances as q → 2. As represented in Fig. 3, for the case k_y = 0.2 the local transient growth attains the value g(t_s) ∼ 10³, whereas at the same time G_max ≈ 2250 for R₀₅ = 2000 and R_b ≪ R₀₅ (see Fig. 14). At first, we check the transient growth of global spirals with m = 1 in the model of a disc with uniform surface density, which has been denoted as P2 (see Fig. 11). Following the strategy adopted for a disc in the Newtonian potential, let us give an interpretation in terms of the local approach. Similar to the case of N2, the optimal initial spiral upwinds at short time-spans while t < t_max, where t_max ≈ 20 is estimated at r_i. After that, the spiral starts shifting outwards, which causes a decline of optimal growth. Now this decline of G is much stronger than in the Newtonian case (cf. Fig. 7) because the value of optimal growth of large-scale vortices falls as one approaches the Keplerian shear rate. In the top-left panel of Fig. 11, it is shown that for t ≳ 20 the analytical curve is in moderate agreement with the numerical curve of G(t) generated for global perturbations with m = 10, and the agreement becomes excellent when t → ∞. The analytical curve corresponds to |$G_{\max }[ \bar{k}_y=const, \bar{R}=const, q(r_{\rm sp}) ]$| given by equation (B3), where we substitute r_sp evaluated from equation (60) into equation (26). G_max is not constant, due to the dependence q on r in P2. Also, because the analytical profile defined in this way diverges as r_sp → r_i, we truncate it for t ≲ 20 and replace by the numerical value of local G_max = 2250 obtained for k_y = 0.2 and parameters fixed in Fig. 11. Note that the maximum of the hump in G(t) corresponding to m = 10 is in good agreement with the latest value. Additionally, as t → ∞, the optimal growth approaches its Newtonian values for both m = 1 and m = 10 (cf. Figs 7 and 11).

Figure 11.

Same as Fig. 7 but for the model P2.

So, we see that the global perturbations with m = 1 generate a much larger hump in P2 compared with N2. Obviously, this happens because the optimal large-scale vortex corresponding to short time-spans is localized in the inner region of disc where q is close to 2. Also, equation (60) yields r_sp ≈ 3 and r_sp ≈ 4, respectively, for the dashed and dot-dashed curves in the bottom panel of Fig. 11, which is again in good agreement with the global numerical results obtained even for m = 1.

In the case of model P2, the modal growth emerges for δ < 0.3. Similar to the model N2, the optimal perturbation at short time is a leading spiral exhibiting a transient growth of energy (see the top-right panel of Fig. 12). Then, for a sufficienly long t, we find that the optimal initial profile of perturbations converges to a rigidly rotating pattern (see the bottom panel of Fig. 12). We check that its increment has the maximum γ ≈ 0.03 located at δ ≈ 0.05 for the particular values of free parameters chosen to produce Fig. 12. The angular velocity of the unstable mode pattern represented in Fig. 12 can be obtained via the Fourier decomposition of, say, ℜ[δh(t)] at some radius. We find that it corresponds to the corotation radius at r_c ≈ 1.04, which is close to the inner boundary of the disc. Yet, the critical layer of the mode is located inside the flow, so the necessary condition for modal energy transfer from background to perturbations is fulfilled. In spite of the fact that for R₀₅ = 2000 the instability exists in the substantial range of δ and the unstable mode has a noticeable increment, as one proceeds to larger Reynolds numbers typical for astrophysical discs, the transient growth starts to dominate the modal growth at any δ for the same reasons as in the case of model N2.

Figure 12.

Same as Fig. 8 but for the model P2.

Finally, the numerical consideration of global perturbations reveals that at R₀₅ = 2000 the disc specified by the model P1 is linearly stable for all δ, whereas the optimal growth is substantially suppressed in comparison with the model P2, and even more in comparison with the local G_max for q = 2 (see Fig. 14). The curve of G(t) and the behaviour of optimal initial spirals with m = 1 can be seen in Fig. 13. Clearly, the optimal growth is suppressed in P1 because the optimal spiral corresponding to its maximum value is located farther from the inner boundary of the disc than a similar spiral in P2 (see the dashed curve in the top-right panel of Fig. 11 and the solid curve in the bottom panel of Fig. 13). In the former situation, the spiral is centred at r_sp ≈ 1.5, where q ≈ 1.8, while in the latter situation the spiral is centred at r_sp ≈ 3, where q ≈ 1.6. For the local large-scale vortices, such a difference in shear rate corresponds to a drop in transient growth by a factor of ∼3. We see that for perturbations with m = 1 it provides a larger difference in the sizes of humps in G(t) in the top-left panels of Figs 11 and 13. The reason for a spiral to shift out of the inner region of the disc is that in P1 viscosity enhances as r → r_i, which excludes the contribution of the region with super-Keplerian shear rate close to q = 2 to the transient growth of global perturbations.

Figure 13.

Same as Fig. 9 but for the model P1.

As in the other models considered previously, we see a moderate agreement of G(t) obtained for global perturbations m = 10 with its analytical counterpart defined by equation (B3) together with definitions (58) and (59), the condition (60) and specifications of the model given in Section 5.1.1 (see the top-left panel of Fig. 13). Additionally, equation (60) yields r_sp ≈ 4.5 and r_sp ≈ 6.5, respectively, for the dashed and dot-dashed curves in the bottom panel of Fig. 13, which is also in agreement with global calculations for m = 1.

We finish the description of global optimal growth in model P1 by comparing its value for perturbations m = 1 with the local analogue being G_max(k_y) for q = 2 (see Fig. 14). In the case of global vortices, we define G_max as the height of the hump in the curve of G(t) (see the top-left panels in Figs 11 and 13). For P2, G_max is a few times smaller than the local value, yet it rises up as |$\bar{k}_y$| decreases similar to G_max(k_y) until the modal growth replaces the transient growth. On the contrary, with regards to vortices with m = 1, geometrically thin discs specified by the model P1 are subject neither to modal growth nor to transient growth, which would exceed its Newtonian value. Thus, in the model P1 the optimal growth of global perturbations, |$G_{\max }(\bar{k}_y)$|⁠, is similar to local optimal growth for the Newtonian shear rate rather than for q = 2.

Figure 14.

Same as Fig. 10 but the solid and dashed curves represent models P1 and P2, respectively, and the dot-dashed curve represents G_max(k_y) for local perturbations with q = 2.

6 CONCLUSIONS

We study the transient growth of shearing vortices at Keplerian and super-Keplerian shear rates, taking into account a finite accretion disc thickness H. This allows us to consider perturbations of an arbitrary azimuthal length-scale, λ_y, with respect to H and, particularly, vortices with λ_y ≳ H, which we refer to as large-scale vortices in contrast to small-scale vortices λ_y ≪ H that have received most of the attention previously. In order to obtain reliable conclusions on their ability to exhibit a transient growth, we include viscous forces acting on fluid elements as a result of both their shear motions and their divergent motions. Thus, our dynamical model is parametrized not only by a Reynolds number, but additionally by a second Reynolds number corresponding to a non-zero bulk viscosity. At first, it is necessary to take into account the influence of bulk viscosity on the dynamics of shearing vortices, because close to the instant of swing a characteristic radial scale of perturbations is always larger than H and velocity divergence becomes significant. This is of special importance for large-scale vortices as their swing interval may be much longer than the rotational period of the flow. Secondly, vortices with λ_y ∼ H generate shearing density waves that are primarily influenced by the bulk viscosity.

Mathematically, we identify the optimal perturbations exhibiting the highest possible transient growth among all perturbations of the same azimuthal wavelength. The value of maximum growth factor is denoted by G_max. This is an important quantity usually determined in the literature in the context of subcritical transition to turbulence in shear flows. Particularly, it has been demonstrated previously that in spectrally stable rotating shear flows G_max correlates with the point of flow relaminarization. In order to obtain optimal solutions, we employ an iterative loop constructed of the basic and adjoint equations advanced, respectively, forwards and backwards in time. Note that the optimization performed via the iterative loop allows us to consider an unbounded shear flow. In this work, the adjoint equations have been derived for both local and global problems; see the sets of equations (17)–(19) and (53)–(55), respectively.

We confirm that initial perturbations producing G_max are vortical leading spirals, which is strictly shown by ZR14 for inviscid perturbations. In the large shearing box, these spirals are represented by SFH with k_x < 0. Being tightly wound, |$|\tilde{k}_x|\gg 1$|⁠, they can be approximated by balanced solutions first used in Heinemann & Papaloizou (2009). With the help of balanced solutions, we derive an analytical estimate of G_max for the large-scale vortices in viscous flow (see equation B3). It complements the known expression for G_max (see equation B1), which is valid for the small-scale vortices only, in the range of λ_y ≫ H. Because equation (B3) yields |$G_{\max }\propto k_y^{2/3}$|⁠, there is no problem of divergence of optimal growth for long wavelength perturbations in the unbounded shear. It is also important that the account of finite H (and, thus, finite sound speed a_*) in the perturbed flow enables one to choose natural dimensionless variables in the problem, when k_y is expressed in units of H⁻¹ rather than in units of an artificial L⁻¹ (cf. Appendix A1 and Section 2.1).

We start the numerical analysis by inspecting the influence of bulk viscosity on the growth of optimal SFH. The results are shown in Fig. 1, where we find a sonic bump on the curves G_max(k_y). As soon as the bulk viscosity becomes larger and larger, the sonic bump caused by the excitation of shearing density waves shrinks and the profile of G_max converges to the profile of g_s(k_y), which is the growth factor of the same optimal SFH evaluated at the instant of swing. The suppression of excited density waves is illustrated in Fig. 2. Importantly, even a huge value of bulk viscosity does not affect the stage of transient growth as long as the vortical spiral is being contracted by the shear (i.e. it does not cause damping of the shearing vortex itself). Moreover, this holds throughout the whole range of azimuthal wavenumbers. Note that in Appendix C we additionally give an analytical interpretation of this property of shearing vortices. Having obtained this result, we choose g_s to be a measure of transient growth of vortices in a compressible rotating flow. Obviously, g_s is a quantity independent of a certain value of bulk viscosity.

It is shown that for R₀₅ = 2000 the growth of shearing vortices in the boundless Keplerian shear attains a maximum value g_s ∼ 100 at k_y ∼ 0.3 expressed in units of H⁻¹, which is substantially larger in comparison with the corresponding G_max ∼ 20 obtained by MAN05 in an incompressible flow between the walls. The optimal growth declines slower as k_y → 0, rather than as k_y → ∞, which is in agreement with analytical estimates (B1) and (B3) (see the top panel of Fig. 1). In the case of a much larger Reynolds number typical for astrophysical discs, the turbulent energy reservoir can be, in principle, replenished by perturbations with λ_y ≳ H, provided that there is a subsequent delivery of their energy to SFH with negative k_x < 0 on smaller azimuthal scales λ_y < H by means of a non-linear interaction between SFHs. Of course, there must be a positive non-linear feedback for regeneration of large-scale leading spirals. Shen, Stone & Gardiner (2006) performed high-resolution 3D numerical simulations in the large shearing box approximation with horizontal size L > H; however, it was still not enough to resolve the region k_y < 1 well, where the transient growth appears to be the most prominent. Yet, in the model of incompressible fluid, Lithwick (2009) had come to the conclusion that vortices of Shen et al. (2006) were not elongated enough in the azimuthal direction to survive for longer times. In light of our results, the high-resolution simulations in the box extended in the azimuthal direction could resolve the tightly wound large-scale leading spirals as well as their possible repopulation due to the non-linear feedback on the scale beyond the disc thickness.

At the same time, Fig. 3 demonstrates that approximately in the range k_y < 0.1 in the Keplerian case G_max is produced mostly at the expense of enthalpy perturbation rather than velocity perturbation. It should be noted, however, that this perturbation of enthalpy has nothing in common with its counterpart existing in density waves: its amplitude increases as |$\tilde{k}_x \rightarrow 0$| and, vice versa, decreases as |$\tilde{k}_x$| of SFH starts growing back. Also, as suggested by balanced solutions, the perturbation of enthalpy entering the large-scale vortex is proportional to the potential vorticity perturbation. For these reasons, it can be argued that even the large-scale SFH consisting mostly of |$\hat{W}$| does not lead to generation of shocks but more likely reorganizes a picture of non-linear interaction between SFHs, presumably, via the term |$\nabla (\rho _1 \boldsymbol {u})$| in the continuity equation. For example, the last term must be significant if it is evaluated for density perturbation generated by the swinging leading spiral with k_y ≪ 1 and velocity perturbation generated by the swinging leading spiral with k_y ≳ 1.

Considering super-Keplerian shear rates 1.5 < q < 2.0, we find that the numerically obtained g_s is also matched well by our analytical estimate of the optimal growth (B3) (see Fig. 3). Thus, with the account of all perturbation scales as compared to H, it is confirmed that there is no such abrupt depression of G_max that takes place for perturbations λ_y ≪ H with the onset of stable epicyclic motions in the flow. For example, the shear rate q = 1.9 provides a maximum growth of shearing vortices two orders of magnitude larger than found by MAN05 in incompressible flow. Notably, the fraction of kinetic energy, g_ks/g_s increases as we approach q = 2 for fixed k_y. This is shown also in Fig. 4, where the amplitude of |$\hat{u}_x$| taken at the instant of swing of optimal SFH rises faster than the amplitude of |$\hat{W}$| as q → 2. Additionally, as soon as q → 2, maxima of g_s(k_y) and g_ks(k_y) shift to smaller k_y and strictly for q = 2 the growth factors diverge almost as |$\propto k_y^{-4/3}$| exceeding the value of G_max produced by streamwise streaks (see fig. 2 d of MAN05), at k_y ≲ 0.1. At this point, let us remark that Rincon, Ogilvie & Cossu (2007) employed a sophisticated non-linear continuation method to show rigorously that already in the case of uniform angular momentum rotation, q = 2, the non-linear feedback sustaining the linear transient growth of perturbations cannot be similar to the one described for plane shear flows and observed in weakly rotating cyclonic flow as well as in Rayleigh-unstable anticyclonic flow. Partly, this is related to the fact that the finite-amplitude streamwise rolls generated by the anti-lift-up effect distort the rotational symmetry of the flow in sharp contrast to what takes place in the non-rotating lift-up case. At the same time, the anti-lift-up itself ceases to operate shortly after we shift to q < 2. This suggests that, at least with the account of disc finite scaleheight, the bypass mechanism of transition may be based on the swing amplification of shearing spirals not only at super-Keplerian shear rate, but also in the uniform angular momentum rotation.

The above results allow us to put a lower limit on the Reynolds number of subcritical transition to HD turbulence (see equation 24). We suggest that the profile R_T(q) may be followed by α(q) (see equation 28), where α is the dimensionless azimuthal stress emerged from a hypothetical turbulence driven by growing shearing vortices. The dependence (28) is steeper than obtained previously in simulations of supercritical MHD turbulence and, for example, implies a substantially larger accretion rate corresponding to the Eddington luminosity in relativistic thin disc models.

An azimuthal wavenumber of local perturbations cannot be as small as the disc aspect ratio, because k_y ∼ δ would correspond to SFH with λ_y ∼ r₀, which is no longer local. Therefore, a study of transient growth of large-scale shearing vortices would not be complete without an investigation of global azimuthal harmonics of perturbations ∝exp (imφ) with m ∼ 1. In Section 5, we additionally consider a global 2D dynamics of shearing vortices described by the set of equations (48)–(50) supplemented by global viscous terms (D1)–(D4). We provide a detailed construction of several global background models, we introduce boundary conditions for both the basic and the adjoint variables and we give brief notes on the numerical method. Note that here we leave out the problem of global density wave excitation, setting R_b ≪ R in all global calculations. In contrast to the local vortex, the global one must be determined not only by the shear rate q and the disc scaleheight H, but additionally by radial profiles of (i) a rotational frequency, (ii) a background vorticity, (iii) a surface density and a viscosity and (iv) a shear rate itself, for example, in relativistic discs. The global optimal transient growth turns out to be more or less sensitive to all these factors, so an exact shape of global |$G_{\max }(\bar{k}_y)$| depends on a particular disc model. Nevertheless, we manage to understand the behaviour of global transient spirals with m ≫ 1, employing interpretation in terms of the local approach and we obtain a quantitative agreement between the shape of G(t) for a variety of disc models and corresponding analytical estimations. The optimal spirals with m ∼ 1 produce G(t), which deviates from the analytical curve. However, they basically preserve the qualitative properties of spirals with m ≫ 1. It turns out that global optimal vortex is always radially localized and, as one goes to longer time-spans, the behaviour of its initial profile is a certain combination of the spiral upwind/unwind and spiral shift along the radial direction. In particular, we manage to predict analytically the locations of initial optimal vortices m = 1 obtained for different optimization time-spans in any of the disc models we used (see Figs 7, 9, 11 and 13). In order to construct an analytical G(t) we assume that (i) the optimization time t is equal to the time elapsed before the instant of swing t_s, and that (ii) as soon as t becomes not less than t_max evaluated at the inner boundary of disc r_i, it remains equal to t_max evaluated at r_sp > r_i, where r_sp is an approximate localization of global spiral corresponding to this optimization time-span. Note that assumption (i) is well justified for shearing vortices. The assumptions (i) and (ii) allow us to formulate the condition (60), which provides us with an estimation of r_sp(t) and, subsequently, with an estimation of G(t), which we suppose to be equal to local value of optimal growth at r_sp. We find that optimal spirals corresponding to global G_max(m = 1) obtained in models with non-uniform viscosity and surface density profiles are shifted outwards in comparison with similar spirals obtained in models with a constant viscosity and surface density. As soon as one includes relativistic effects (cf. the models P1 and P2), this causes the suppression of G_max down to a Keplerian value (see Fig. 14). At the same time, in the case of Newtonian dynamics we find the global G_max to attain approximately the same value as its local analogue (cf. the models N1 and N2 in Fig. 10). Consequently, the Keplerian value of optimal growth produced by the large-scale shearing vortices is robust with the account of global effects, whereas the enhancement of locally determined G_max obtained at the super-Keplerian shear rate does not manifest itself on a global scale in a disc with viscosity attaining relatively high values in regions with super-Keplerian rotation. Additionally, we consider the modal growth of global perturbations in sufficiently thin discs specified by the models N2 and P2. The most unstable mode is obtained by the optimization method applied at long time-spans. Its profile is shown is Figs 8 and 12 for models N2 and P2, respectively. As can be seen, the unstable mode is concentrated close to the inner boundary of disc. It is beyond the scope of this work to investigate the physical nature of modal growth, but we check that corotation of the mode is located inside the flow and the increment attains a few per cent larger value in the inviscid case. Also, the modal growth ceases rapidly as we increase m. Thus, we suspect that the revealed instability is explained by the Landau mechanism, implying that the energy flux from the background flow to the acoustic mode is caused by a non-zero gradient of background vorticity at the corotation point. Let us point out that this sort of modal growth competes with transient growth at Reynolds numbers equal to a few thousand only. As G_max ∝ R^2/3, we expect that in astrophysical situations R ∼ 10¹⁰–10¹³ the transient growth of global perturbations dominates the modal growth in a disc of an arbitrary thickness.

6.1 Future prospects

This work has demonstrated the importance of considering the transient growth at all scales with respect to the disc scaleheight. However, it has been carried out in a simplified 2D approximation of perturbation dynamics. Accordingly, the optimal vortices should be subsequently studied with the account of vertical motions. It is worth checking whether 3D optimal perturbations with either λ_y ≪ H or λ_y ≳ H in radially unbounded flow have a non-trivial axial structure and kinetic energy transferring into vertical motions in connection with the results obtained by Maretzke et al. (2014). Indeed, within the model of fluid between the rotating cylinders, Maretzke et al. (2014) have compared the optimal growth in various regimes of a spectrally stable incompressible Taylor–Couette flow. They found that in contrast to the quasi-Keplerian regime, the cyclonic and the counter-rotating regimes allow for the axial dependence of optimal vortices as well as for a significant part of the kinetic energy to be contained in vertical motion. Though the physical reasons for such a difference remain unknown, this can probably be related to a subcritical HD transition in cyclonic and counter-rotating regimes occurring at ordinary values of Reynolds number, contrary to the case of a quasi-Keplerian regime. Thus, an interesting issue would be to find out whether the optimal – or nearly optimal – perturbations in an unbounded quasi-Keplerian flow have the non-trivial axial dependence, and to explain its physical nature. Note that with regards to the large-scale vortices, λ_y ≳ H, the corresponding result may be different as one formally considers the vertically homogeneous or more realistic vertically baratropic flow.

While there is expectation to reveal the bypass transition within the most basic model of quasi-Keplerian flow, astrophysical discs can be stratified and there are a number of refined models shown to exhibit linear and non-linear HD activity due to various thermodynamical inhomogeneities. As a rule, this HD activity results in a generation of vortices of considerable amplitude, which are usually invoked to resolve the problems of angular momentum transport and planetesimal formation inside non-ionized parts of protoplanetary discs. Some of these are the product of modal instabilities, such as Rossby wave instability (Lovelace et al. 1999; Li et al. 2000, 2001) or vertical shear instability (Urpin 2003; Nelson, Gressel & Umurhan 2013; Richard, Nelson & Umurhan 2016), whereas others are attributed to subcritical instabilities emerging either due to the radial entropy gradient (Lesur & Papaloizou 2010), or due to the vertical entropy gradient (Marcus et al. 2015). In view of the subcritical behaviour of stratified Keplerian flow, it is important to note that a novel type of linearly growing vortical 3D perturbations appears to be the result of the so-called strato-rotational balance (see Tevzadze et al. 2003; Tevzadze, Chagelishvili & Zahn 2008). The existence of these 3D vortices is closely related to a special invariant of stratified fluid dynamics equal to the scalar product of the entropy gradient and the vorticity. These vortices have been suggested as transient perturbations for the bypass transition scenario in stratified discs; however, note that they degenerate into a trivial solution within the unstratified model considered in this paper. Also, interestingly, the presence of radial stratification allows for the production of vortices by means of purely linear dynamics (see Tevzadze et al. 2010). The 2D shearing density waves acquire the ability to generate 2D shearing vortices through indirect coupling with the assistance of shearing entropy waves. It is not yet clear what role this novel and complicated non-modal linear dynamics can play in the subcritical instability of stratified Keplerian flows, motivating us to consider the corresponding 3D optimization problem for perturbations with the account of the finite disc thickness.

Furthermore, it seems relevant to study the excitation of density waves within the 3D optimization problem in baratropic and vertically stratified backgrounds. Particularly, this would be necessary to do in the large-scale limit λ ≫ H in order to establish the link with the work done by Umurhan et al. (2006) and, later, by Rebusco et al. (2009). Employing the asymptotic time-dependent model of geometrically thin viscous Keplerian disc, these authors found that the vertical sound waves exhibit non-modal growth that is coupled to planar epicyclic oscillations.

At last, it is desirable to perform high-resolution HD simulations in the large shearing box approximation in order to check whether the critical value of the shear rate q < 2 corresponding to turbulence decay found by Hawley, Balbus & Winters (1999) becomes smaller as one considers the non-linear dynamics of shearing vortices with k_y ≲ 1. It would be particularly interesting to address the issue of significant density perturbations produced by large-scale vortices in the context of non-linear interaction between SFH with different k_x and k_y. We point out that these density perturbations are proportional to a perturbation of the potential vorticity being time-independent in 2D inviscid compressible dynamics, so they are qualitatively different from the density perturbations generated by shearing density waves. Also, it is of course an open question in which way these density perturbations would affect a hypothetical non-linear transverse cascade, which probably provides a positive feedback for transiently growing leading spirals at very high Reynolds numbers in the quasi-Keplerian flows.

6.2 Final remarks

Let us note that the issues we have mentioned in this paper are closely related to the applied problem of stability and transition to turbulence in hypersonic shear flows. The stability and transition in hypersonic shear flows is a developing area of experimental, theoretical and numerical fluid dynamics with a particular emphasis on the design of flying vehicles; see the book by Gatski & Bonnet (2009) and the reviews by Fedorov (2011) and Zhong & Wang (2012). These flows are usual phenomena in space and should be understood better in astrophysical science and, mainly, in the theory of astrophysical discs.

We would like to finish this paper by making link with the neighbouring field of magnetohydrodynamic (MHD) turbulence in accretion discs. Actually, it has become obvious in the past few years that the bypass concept of the transition to turbulence and mechanisms of transient growth are relevant not only to the HD problem. The work by Squire & Bhattacharjee (2014a) and Squire & Bhattacharjee (2014b), and see also Nath & Mukhopadhyay (2015), has shown that the established practice to oppose the subcritical HD transition to supercritical MHD transition to turbulence in astrophysical discs was misleading in a particular sense. It turned out that even in a spectrally unstable magnetized Keplerian shear the transient growth prevails the modal growth by orders of magnitude at certain short time-spans. Thus, in fact, it is not yet clear what linear mechanism gives rise to a sustained MHD turbulence even in hot ionized accretion discs: either the MRI or the transient growth of magnetized shearing vortices. Moreover, a non-modal amplification of the magnetic field and subsequent positive non-linear feedback analogous to SSP in the context of HD transition in plane shear flows has been suggested by Rincon et al. (2008) and Riols et al. (2013) to give a basic explanation of the subcritical magnetorotational dynamo action. The latter is crucial for turbulization of a spectrally stable Keplerian background with zero net magnetic flux. The above extension of non-modal concepts on to MHD turbulence in astrophysical flows should be completed by the subcritical transition to turbulence via the non-linear transverse cascade, providing positive feedback to transiently growing SFH in 2D magnetized plane flow, recently demonstrated by Mamatsashvili et al. (2014). Bearing in mind the current situation in the problem of the transition to turbulence in (quasi-)Keplerian flow, it seems to be fruitful to address the remaining issues in both HD and MHD problems using a general non-modal approach.

Acknowledgments

We thank P. B. Ivanov for careful reading the manuscript and for important remarks. DNR was supported by grant RSF 14-12-00146 when writing Section 5 of this paper. VVZ was supported by RFBR grant 15-02-08476. Also, the study was in part supported by the M. V. Lomonosov Moscow State University Programme of Development and in part by the Programme 7 of the Presidium of the Russian Academy of Science.

1

It is assumed hereafter that k_y > 0.

2

We also show the growth factors of corresponding optimals at the instant of swing in order to extinguish the excited density waves (see §4.2).

3

This explanation should be taken with a caution that one cannot distinguish between vortices and waves in the swing interval and, in particular, at the instant of swing. Strictly, we can only say that g_s becomes considerably larger than the estimate extracted from the balanced solution as soon as ε is not small. This remark is not unnecessary because, as we have already checked, g_s is not affected by high bulk viscosity, R_b ≪ R, even for q = 2 in the region ε ∼ 1.

4

Additionally, we have checked that the change to norm of perturbations with term |$|\hat{W}|^2$| in equation (11) multiplied by a constant p ∈ (0, 1] and the corresponding change of the adjoint equations leads to the same optimal solutions, which are the leading spirals of the vortical SFH with negative k_x.

5

In equation (27), it is assumed that the cross-section of collisions between the particles is independent of a_* and N.

APPENDIX A: SUBCASE OF INCOMPRESSIBLE DYNAMICS

A1 Shearing vortex solution

Let us consider the incompressible limit of transient dynamics also studied by Yecko (2004) and MAN05 in the case of a flow confined between two walls separated from one another by distance 2L across the shear. Here, we deal with an unbounded shear. As a_* → ∞, the continuity equation (3) yields the following restriction on the velocity components

\begin{equation*} \frac{\mathrm{\partial} u_x}{\mathrm{\partial} x} = - \frac{\mathrm{\partial} u_y}{\mathrm{\partial} y}, \end{equation*}

and we are left with two projections of Euler equation (1) and (2).

We change to dimensionless comoving Cartesian coordinates in equations (1)–(3),

\begin{equation} x^\prime = x/L,\qquad y^\prime = (y+q\Omega _0 xt)/L,\qquad t^\prime =\Omega _0 t, \end{equation}

(A1)

where the quantity L is an auxiliary radial scale, which does not have any physical meaning. The coordinates (A1) lead us to a spatially homogeneous set of equations, which have partial solutions in the form of SFH corresponding to the dimensionless wavenumbers k_x and k_y expressed in units of arbitrary L⁻¹. Similarly to Section 2.1, we omit the primes and introduce the shearing radial wavenumber |$\tilde{k}_x \equiv k_x + q k_y t$| and the full wavenumber squared |$k^2 \equiv \tilde{k}_x^2 +k_y^2$|⁠, arriving at the following dimensionless equations for a single SFH:

\begin{eqnarray} \frac{{\rm d} \hat{u}_x}{{\rm d} t} &=& 2\hat{u}_y - {\rm i}\tilde{k}_x \hat{W} - R^{-1} k^2 \hat{u}_x, \end{eqnarray}

(A2)

\begin{eqnarray} \frac{{\rm d} \hat{u}_y}{{\rm d} t} &=& -(2-q)\hat{u}_x - {\rm i}\tilde{k}_y \hat{W} - R^{-1} k^2 \hat{u}_y, \end{eqnarray}

(A3)

provided that

\begin{equation} \tilde{k}_x \hat{u}_x = -k_y \hat{u}_y. \end{equation}

(A4)

It is assumed that |$\hat{u}_x,\hat{u}_y$| are expressed in units of U ≡ Ω₀L, and |$\hat{W}$| is expressed in units of U². Correspondingly, R ≡ UL/ν is the Reynolds number. Note that MAN05 used Reynolds number

\begin{equation} R_{05} = q R. \end{equation}

(A5)

Using equations (A2)–(A4), we obtain the equation for |$\hat{u}_x$|⁠, which is

\begin{equation} \frac{{\rm d}\hat{u}_x}{{\rm d}t} = -2q \frac{k_y \tilde{k}_x}{k^2} \hat{u}_x -R^{-1} k^2 \hat{u}_x. \end{equation}

(A6)

It can be checked that equation (A6) is reproduced from equation (16) of MAN05 and has the following analytical solution:

\begin{equation} \hat{u}_x (t) = \hat{u}_x(0) \frac{k_x^2+k_y^2}{k^2} \exp \left(-\frac{\tilde{k}^3-k_x^3}{3qk_y R} - \frac{k_y^2 t}{R} \right). \end{equation}

(A7)

The solution (A7) is used to obtain the transient growth factor of incompressible perturbations for particular values of k_x, k_y, R and q; see Section A2 for further consideration.

A2 Transient growth factor

As usual, the size of incompressible perturbations in an unbounded shear flow is measured by a surface density of their kinetic energy, E_k, which leads to a simple expression for a single SFH:

\begin{equation} E_k = \frac{1}{2}(|\hat{u}_x|^2 + |\hat{u}_y|^2) = \frac{1}{2} \frac{k^2}{k_y^2} |\hat{u}_x|^2. \end{equation}

(A8)

Substituting equation (A7) into equation (A8) yields the following growth factor of SFH:

\begin{equation} g \equiv \frac{E_k(t)}{E_k(0)} = \frac{k_x^2+k_y^2}{k^2} \exp \left(-2\frac{\tilde{k}^3-k_x^3}{3qk_y R} - 2\frac{k_y^2 t}{R} \right). \end{equation}

(A9)

Given the particular values of k_y, R and q, equation (A9) has a global maximum over all k_x and t, which we denote here as

\begin{equation} G_{\max } \equiv \max _{\forall k_x, \forall t}\, g, \end{equation}

(A10)

because this quantity is analogous to the maximum optimal growth defined in Section 3 (see equation 21).

APPENDIX B: ANALYTICAL ESTIMATIONS OF OPTIMAL GROWTH

By means of the analytical supplement to their numerical optimization results, MAN05 demonstrate that in the case of incompressible dynamics G_max(k_y) acquires a maximum owing to the presence of rigid boundaries (see also the results of Yecko 2004 in their fig. 3b). In Fig. B1, we take the Keplerian shear, q = 3/2, R₀₅ = 2000 and plot equation (77) of MAN05 (dot-dashed curve) together with our quantity (A10) obtained numerically (solid curve). It can also be checked that equation (77) of MAN05 taken for a boundless shear,

\begin{equation} G_{\max } = \left(\frac{qR}{k_y^2} \right)^{2/3} \, {\rm e}^{-2/3}, \end{equation}

(B1)

virtually coincides with our numerical curve in Fig. B1. Equation (B1) can be easily obtained from equation (A9) setting |$\tilde{k}_x=0$| (the swing instant of SFH) and taking the limit |k_x| ≫ k_y. Thus, in general, G_max → ∞ as k_y → 0 in the absence of boundaries. Physically, this takes place due to the simple fact that for a tightly wound SFH with |k_x| ≫ k_y both viscous dissipation time-scale, t_ν, and transient growth duration, t_tg, become longer as we proceed to small azimuthal wavelengths, whereas the (inviscid) growth factor increases monotonically with time (see equation A9). As discussed, for example, by RZ15, equation (B1) can be approximately derived suggesting that G_max corresponds to equality t_ν = t_tg, which gives the time when the maximum growth occurs as

\begin{equation} t_{\max } \approx R^{1/3} (qk_y)^{-2/3}. \end{equation}

(B2)

Note that estimation (B2) provides exactly the same value that is obtained when looking for a maximum of equation (A9) in order to derive an approximate G_max given by equation (B1). Then, using the inviscid form of g ≈ (qt)², we recover equation (B1) up to the damping factor exp (−2/3).

Figure B1.

Curves of G_max in the subcase of incompressible dynamics of 2D perturbations for R₀₅ = 2000 and q = 3/2. The solid curve is obtained numerically using the definition (A10). The dashed curve represents equation (B1) and the dot-dashed curve corresponds to equation (77) of MAN05 taken with k_{X, min}L = 1.7 in their notations. The dotted curve represents the approximate G_max given by equation (B3) for the case of finite disc thickness and k_y normalized by H rather than by L.

Of course, the growth factor (A9) and so G_max do not depend on the scaling factor L. As we see from equation (B1), the transient growth value is controlled by the ratio |$R/k_y^2 \sim (\lambda _y/l_{\rm p})^2$|⁠, where λ_y is the characteristic azimuthal length of perturbations and l_p is the mean free path of gas particles that determine microscopic viscosity. In fact, |$\lambda _y^{-1} \sim k_y L^{-1}$| of SFH exhibiting G_max cannot be too low as |k_x| ∼ (qk_yR)^1/3 also decreases; however, |k_x|L⁻¹ cannot be less than H⁻¹, as in the last case we start dealing with a patch of disc comparable to its thickness in the direction of shear, which forces us to take into account the effects of compressibility (see, e.g. the discussion about the SFH swing interval in section 3.2.1 of ZR14). Thus, of course, we cannot reach an arbitrary high G_max at the expense of the λ_y value. It turns out (see section 2.3 of RZ15) that changing to a model with the account of finite disc thickness, we can derive an approximate expression for G_max in the opposite case of λ_y ≫ H (i.e. k_y ≪ 1) employing the balanced solutions to equations (7)–(9) (see Heinemann & Papaloizou 2009). The inviscid result is g_opt ≈ 4κ⁻⁴(k_yqt)², where κ² = 2(2 − q) is the dimensionless epicyclic frequency squared expressed in units of |$\Omega _0^2$|⁠. Substituting t_max given by equation (B2) into this result, we obtain

\begin{equation} G_{\max } \approx \frac{1}{(2-q)^2} (k_y q R)^{2/3} \, {\rm e}^{-2/3}, \end{equation}

(B3)

where the damping factor exp (−2/3) has been added because we expect a similar action of viscous forces on shearing vortices having both k_y ≫ 1 and k_y ≪ 1; see a detailed justification of this assumption in Section 4.1. As we see, one obtains an opposite dependence of the maximum optimal growth on k_y (i.e. G_max decreases as we approach larger azimuthal scales). Thus, there is no problem of divergence of transient growth factor towards k_y → 0 (see Fig. B1). At the same time, as seen in Fig. B1, one expects much larger transient growth of columnar perturbations (i.e. perturbations independent of vertical direction) in the domain of λ_y ≳ H, than it might seem within the model of Yecko (2004) and MAN05.

APPENDIX C: SUPPRESSION OF SHEARING DENSITY WAVES FOR TIGHTLY WOUND SPIRALS

In the limit |$|\tilde{k}_x|\gg 1$| – or, equivalently, outside the swing interval of SFH (see ZR14) – the inviscid shearing vortices are approximately given by balanced solutions obtained from wave equations for |$\hat{u}_{x,y}$| and |$\hat{W}$|⁠; see, for example, equations (27)–(29) of RZ15. Their explicit form is given by equations (35)–(37) of RZ15. The non-zero viscosity leads to the following wave equations for |$\hat{u}_y$| and |$\hat{u}_x$|⁠:

\begin{eqnarray} \ddot{\hat{u}}_x &=& - k_y I_\nu - K \hat{u}_x - 2{\rm i}qk_y \hat{W} \nonumber \\ &&-\, R^{-1} k^2 \left[ \dot{\hat{u}}_x + \hat{u}_y + 2q \frac{\tilde{k}_x k_y}{k^2} \hat{u}_x \right] \nonumber \\ &&-\, \left(R^{-1}/3 + R_{\rm b}^{-1}\right) \left[ \frac{{\rm d}}{{\rm d}t}\left(\tilde{k}_x \, {\rm i} \dot{\hat{W}} \right) + 2k_y \, {\rm i} \dot{\hat{W}} \right], \end{eqnarray}

(C1)

\begin{eqnarray} \ddot{\hat{u}}_y &=& \tilde{k}_x I_\nu - K \hat{u}_y \nonumber \\ &&-\, R^{-1} k^2 \left[ \dot{\hat{u}}_y - (2-q)\hat{u}_x + 2q \frac{\tilde{k}_x k_y}{k^2} \hat{u}_y \right] \nonumber \\ &&-\, \left(R^{-1}/3 + R_{\rm b}^{-1}\right) [k_y \, {\rm i} \ddot{\hat{W}} - (2-q)\tilde{k}_x \, {\rm i} \dot{\hat{W}}]. \end{eqnarray}

(C2)

Here, K ≡ k² + κ² and |$I_\nu \equiv \tilde{k}_x \hat{u}_y - k_y \hat{u}_x + {\rm i}(2-q)\hat{W}$| is SFH of the potential vorticity perturbation; see the corresponding expressions, for example, of Bodo et al. (2005) or Tevzadze et al. (2003) for the case of barotropic perturbations. I_ν does not depend on time in the inviscid limit; however, in the presence of non-zero viscosity this is no longer true. To emphasize this fact, we add a subscript ν to the notation of the potential vorticity perturbation. It is not difficult to construct the equation for I_ν using the set of equations (7)–(9):

\begin{equation} \dot{I}_\nu = - R^{-1} k^2 I_\nu . \end{equation}

(C3)

Here, the second term on the right-hand side of equation (C3), explicitly |${\rm i} (2-q) R^{-1} k^2 \hat{W}$|⁠, has been omitted for the following reason. In the problem we consider, the viscosity is small, R ≫ 1, R_b ≫ 1. Thus, in the leading orders in R⁻¹ and |$R_{\rm b}^{-1}$| one can use an inviscid balanced solution for |$\hat{W}$| (see equation 37 of RZ15), which gives |$|\hat{W}/I_\nu |\ll 1$| provided that |$|\tilde{k}_x| \gg 1$|⁠. This means that the omitted term is much less than the main term on the right-hand side of equation (C3). Similarly, in the leading orders in R⁻¹ and |$R_{\rm b}^{-1}$|⁠, it is sufficient to take the slow varying solution of the inhomogeneous equation (C2) in its inviscid form,

\begin{equation} \hat{u}_y = \frac{\tilde{k}_x}{K} I_\nu , \end{equation}

(C4)

but with I_ν decaying as described by equation (C3). This is the decay of I_ν, which yields an additional factor exp (−2/3) in equation (B3) that is evaluated at t_s = −k_x/(qk_y) with t_s additionally equated to equation (B2). Correspondingly, |$\hat{u}_x$| is given by

\begin{equation} \hat{u}_x = -\frac{K+4q}{K^2+4q^2k_y^2} k_y I_\nu ; \end{equation}

(C5)

see equation (35) of RZ15. We see that as far as |$|\tilde{k}_x|$| is large, the radial velocity perturbation is much less than the azimuthal one in the case of tightly wound vortices (i.e. when |$|\tilde{k}_x| \gg k_y$|⁠) and the velocity perturbation divergence is small:

\begin{equation} {\rm i}\dot{\hat{W}} \sim \frac{k_y}{\tilde{k}_x^2} \hat{u}_y. \end{equation}

(C6)

The tightly wound vortical spirals are of primary interest here, as we are dealing with high R and we focus on k_y ≲ 1, which leads to t_max ≫ 1 (see equation B2). Further, because t_max approximately equals t_s of an optimal SFH exhibiting G_max, at the initial moment |k_x| ≫ k_y and the same is about |$|\tilde{k}_x|$| until t > 0 is not too close to the instant of swing.

Equation (C6) implies that terms in the second square brackets on the right-hand sides of equations (C1) and (C2) are smaller at least by a factor |$|\tilde{k}_x|^{-3}$| than the first and the second terms therein. Taking into account that terms proportional to |$R_{\rm b}^{-1}$| are absent in equation (C3), we come to the conclusion that bulk viscosity cannot affect the dynamics of vortical SFH even if R_b ∼ 1.

The opposite situation takes place with density waves. These are rapidly oscillating solutions to the homogeneous equations (27)–(29) of RZ15. This is why, as long as |$|\tilde{k}_x|\gg 1$|⁠, |$\dot{\hat{u}}_x \sim |\tilde{k}_x| \hat{u}_x$| and so is about |$\dot{\hat{u}}_y \sim |\tilde{k}_x| \hat{u}_y$| in the inviscid case. We see that in this situation at least the term |$R_{\rm b}^{-1} \tilde{k}_x^2 \dot{\hat{u}}_x$| becomes the leading term among all viscous terms in equation (C1) and it exceeds the inviscid terms as soon as |$R_{\rm b} \lesssim |\tilde{k}_x|^2$|⁠.

APPENDIX D: SUPPLEMENTARY MATERIAL TO THE GLOBAL PROBLEM

D1 Viscous terms

Viscous terms entering the dynamical equations (48)–(50) explicitly are

\begin{eqnarray} N_r&=&\frac{\nu }{r}\frac{\mathrm{\partial} }{\mathrm{\partial} r}\left(r\frac{\mathrm{\partial} \delta v_r}{\mathrm{\partial} r}\right)-\frac{\nu \left(m^2-1\right)}{r^2}\delta v_r+\frac{2}{\Sigma r}\frac{\mathrm{\partial} \delta v_r}{\mathrm{\partial} r}\frac{{\rm d}}{{\rm d}r}\left(\nu \Sigma r\right), \nonumber\\ \end{eqnarray}

(D1)

\begin{eqnarray} N_{\varphi }&=&-\frac{\nu r\Omega ^{\prime }}{\Sigma a_*^2}r^{-({2n}/{n+1})}\frac{{\rm d}}{{\rm d}r}\left(\Sigma r^{({2n}/{n+1})}\right)\delta h+\nonumber\\ &&\!\!+\frac{1}{\Sigma r^2}\frac{\mathrm{\partial} }{\mathrm{\partial} r}\left(\!\nu \Sigma r^3 \frac{\mathrm{\partial} }{\mathrm{\partial} r}\frac{\delta v_{\varphi }}{r}\!\right)\!+\!\frac{{\rm i}m}{\Sigma r^3}\frac{{\rm d}}{{\rm d} r}\left(\nu \Sigma r^2\right)\delta v_r\!-\!\frac{\nu m^2}{r^2}\delta v_{\varphi },\nonumber\\ \end{eqnarray}

(D2)

\begin{eqnarray} B_r=\frac{1}{\Sigma }\frac{\mathrm{\partial} }{\mathrm{\partial} r}\left[\Sigma \left(\nu _{\rm b}+\frac{\nu }{3}\right){\rm div}\ \delta \boldsymbol {v}\right]-\frac{1}{\Sigma r^2}\frac{{\rm d}}{{\rm d} r}\left(\nu \Sigma r^2\right) {\rm div}\, \delta \boldsymbol {v}, \nonumber\\ \end{eqnarray}

(D3)

\begin{equation} B_{\varphi }=\frac{{\rm i}m}{r}\left(\nu _{\rm b}+\frac{\nu }{3}\right){\rm div}\, \delta \boldsymbol {v}. \end{equation}

(D4)

The norm (52) yields the following viscous terms in the adjoint equations (53)–(55):

\begin{eqnarray} \tilde{N_r}= -\frac{\nu }{r}\frac{\mathrm{\partial} }{\mathrm{\partial} r}\left(r\frac{\mathrm{\partial} \delta \tilde{v}_r}{\mathrm{\partial} r}\right) +\frac{\nu \left(m^2-1\right)}{r^2}\delta \tilde{v}_r -\frac{2}{\Sigma r}\frac{\mathrm{\partial} \delta \tilde{v}_r}{\mathrm{\partial} r}\frac{{\rm d}}{{\rm d}r}\left(\nu \Sigma r\right); \!\!\!\!\nonumber\\ \end{eqnarray}

(D5)

\begin{eqnarray} \tilde{N}_{\varphi }&=&-\frac{1}{\Sigma r^2}\frac{\mathrm{\partial} }{\mathrm{\partial} r}\left(\nu \Sigma r^3 \frac{\mathrm{\partial} }{\mathrm{\partial} r}\frac{\delta \tilde{v}_{\varphi }}{r}\right)-\frac{{\rm i}m}{\Sigma r^3}\frac{{\rm d}}{{\rm d}r}\left(\nu \Sigma r^2\right)\delta \tilde{v}_r\nonumber \\ &&+\,\frac{\nu m^2}{r^2}\delta \tilde{v}_{\varphi }, \end{eqnarray}

(D6)

\begin{equation} \tilde{N}_h=\frac{\nu r \Omega ^{\prime }}{\Sigma }r^{-({2n}/{n+1})}\frac{{\rm d}}{{\rm d}r}\left(\Sigma r^{({2n}/{n+1})}\right) \delta \tilde{v}_{\varphi }; \end{equation}

(D7)

\begin{eqnarray} \tilde{B}_r=-\frac{1}{\Sigma }\frac{\mathrm{\partial} }{\mathrm{\partial} r}\left[\Sigma \left(\nu _{\rm b}+\frac{\nu }{3}\right) {\rm div}\, \delta \tilde{\boldsymbol {v}}\right] +\frac{1}{\Sigma r^2}\frac{{\rm d}}{{\rm d}r}\left(\nu \Sigma r^2\right) {\rm div}\, \delta \tilde{\boldsymbol {v}};\nonumber\\ \end{eqnarray}

(D8)

\begin{equation} \tilde{B}_{\varphi }=-\frac{{\rm i}m}{r}\left(\nu _{\rm b}+\frac{\nu }{3}\right) {\rm div}\, \delta \tilde{\boldsymbol {v}}. \end{equation}

(D9)

Note that boundary terms that emerge in the course of derivation of equations (53)–(55) are used to formulate boundary conditions for the adjoint variables (see Appendix D2). Also, similar to what we have in the local problem, the terms in equations (D5), (D6), (D8) and (D9) acquire opposite signs compared to the corresponding terms in equations (D1)–(D4). Additionally, the term reciprocal to term ∝δh in equation (D2) appears in the adjoint equation for |$\delta \tilde{h}$| (see equation 55). Note that the last term vanishes as one turns to local perturbations.

D2 Boundary conditions

Let us impose the boundary conditions for perturbations necessary to advance equations (48)–(50) forward in time and equations (53)–(55) backward in time. As the disc is considered to be radially infinite, only a condition at the inner boundary r = r_i is relevant for the dynamics at a finite time-span.

D2.1 Models N1 and P1

First, let us consider the set of equations for projections of the displacement vector |$\boldsymbol {\xi }$|⁠. This can be done using a general relationship between |$\boldsymbol {\xi }$| and the velocity Lagrangian perturbation (Lynden-Bell & Ostriker 1967, equation 16):

\begin{equation} \frac{\mathrm{\partial} \xi _r}{\mathrm{\partial} t}+{\rm i}m\Omega \xi _r=\delta v_r+\frac{\mathrm{\partial} }{\mathrm{\partial} r}\left(v_r \xi _r\right); \end{equation}

(D10)

\begin{equation} \frac{\mathrm{\partial} \xi _{\varphi }}{\mathrm{\partial} t}+{\rm i}m\Omega \xi _{\varphi }=\delta v_{\varphi }+r\Omega ^{\prime }\xi _r+r v_r \frac{\mathrm{\partial} }{\mathrm{\partial} r}\left(\frac{\xi _{\varphi }}{r}\right). \end{equation}

(D11)

Because v_r tends to infinity at r_i (see equation 37), ξ_r and ξ_φ must vanish at r_i. This implies that the boundary remains unperturbed, which simplifies the boundary conditions. Now, as v_r ∝ D⁻³ in the vicinity of the boundary, the components of the displacement vector must behave at least as

\begin{equation} \xi _r \propto D^3, \end{equation}

(D12)

\begin{equation} \xi _{\varphi } \propto D^4, \end{equation}

(D13)

close to r = r_i.

Boundary conditions to the set (48)–(50) must be that the Lagrangian perturbation of surface force acting on the inner boundary of the flow is equal to zero. A relation between this force and the full stress tensor, σ_ik, can be found, for example, in Landau & Lifshitz (1987, section 15). Because ξ_r = 0, the Lagrangian perturbation of the normal to the boundary is equal to zero. Thus, the boundary condition turns into two independent conditions on the Lagrangian perturbation of σ_ik:

\begin{equation} \Delta \sigma _{rr}=0, \end{equation}

(D14)

\begin{equation} \Delta \sigma _{r\varphi }=0. \end{equation}

(D15)

Integration of equation (D15) over z yields

\begin{eqnarray} {\Sigma \nu \left[r\frac{\mathrm{\partial} }{\mathrm{\partial} r}\left(\frac{\delta v_{\varphi }}{r}\right)+\frac{{\rm i}m}{r}\delta v_r\right]+ \xi _r \frac{{\rm d}}{{\rm d}r}\left(\Sigma \nu r \Omega ^{\prime }\right)}\nonumber \\ {\quad-2\xi _r\frac{{\rm d}H}{{\rm d}r}\nu r \Omega ^{\prime } \rho {|}_{z=H}=0.} \end{eqnarray}

(D16)

In equation (D16), all variables are taken at the inner boundary of the disc. However, because of the relation (38) and limitations on the displacement vector (D12) and (D13), equation (D16) is satisfied for arbitrary δv_r and δv_φ. Thus, equation (D16) does not provide us with any restriction on perturbations at r_i.

At the same time, equation (D14) results in the following equality:

\begin{eqnarray} &&{2\rho \nu \frac{\mathrm{\partial} \delta v_r}{\mathrm{\partial} r}+\rho \left(\nu _{b}-\frac{2}{3}\nu \right)\left[\frac{1}{r}\frac{\mathrm{\partial} }{\mathrm{\partial} r}\left(r \delta v_r\right)+\frac{{\rm i}m}{r}\delta v_{\varphi }\right]}\nonumber \\ &&{\qquad\quad-2\rho \nu \Omega ^{\prime } \xi _{\varphi }-\delta p -\frac{\mathrm{\partial} p}{\mathrm{\partial} r}\xi _r=0.} \end{eqnarray}

(D17)

In equation (D17), all variables are taken at the inner boundary of disc. Equation (D17) allows us to formulate the necessary restrictions on δv_r and δv_φ at r_i, provided that it is supplied by the regularity condition for δh = δp/ρ at the boundary. Indeed, equation (D17) yields

\begin{eqnarray} \delta h&=& \frac{\delta p}{\rho }= \left(\nu _{\rm b}+\frac{4}{3}\nu \right)\left[\frac{\nu _{\rm b}-({2}/{3})\nu }{\nu _{\rm b}+({4}/{3})\nu }\frac{\delta v_r}{r}+\frac{\mathrm{\partial} \delta v_r}{\mathrm{\partial} r}\right]\nonumber \\ &&+\left(\nu _{\rm b}-\frac{2}{3}\nu \right)\frac{{\rm i}m}{r}\delta v_{\varphi }- 2\nu \Omega ^{\prime } \xi _{\varphi }- \frac{\mathrm{\partial} p}{\mathrm{\partial} r}\frac{\xi _r}{\rho }. \end{eqnarray}

(D18)

The last two terms vanish due to limitations (D12) and (D13). As ν diverges as r → r_i, δh is regular at r_i as soon as

\begin{equation} -\!\frac{1}{2}\frac{\delta v_r}{r}+\frac{\mathrm{\partial} \delta v_r}{\mathrm{\partial} r}=0, \end{equation}

(D19)

\begin{equation} \delta v_{\varphi }=0. \end{equation}

(D20)

As mentioned in Section 5.3, in order to obtain the boundary conditions for the adjoint variables, it is necessary to put the boundary terms emerging in derivation of |$\boldsymbol {A}^{\dagger}$| equal to zero for each adjoint equation independently. This results in the following set of the adjoint boundary conditions:

\begin{equation} \delta \tilde{v}_r=0; \end{equation}

(D21)

\begin{equation} \delta \tilde{v}_{\varphi }=0; \end{equation}

(D22)

\begin{equation} \delta \tilde{h}=0. \end{equation}

(D23)

D2.2 Models N2 and P2

In this case, we impose the no-slip boundary conditions at r = r_i:

\begin{equation} \delta v_r=0; \end{equation}

(D24)

\begin{equation} \delta v_{\varphi }=0. \end{equation}

(D25)

For a disc of uniform surface density, the substitutions in derivation of the adjoint equations vanish under the conditions

\begin{equation} \delta \tilde{v}_r=0, \end{equation}

(D26)

\begin{equation} \delta \tilde{v}_{\varphi }=0, \end{equation}

(D27)

at r = r_i. Note that for models N2 and P2 the boundary condition on to |$\delta \tilde{h}$| is not required, as the advective terms are absent in this case; see the explanation in Section 5.4.