Transient growth of perturbations on scales beyond the accretion disc thickness

TL;DR

This paper investigates the transient growth of shearing vortices in quasi-Keplerian flows, revealing how vortex scale and viscosity influence turbulence development beyond the accretion disc thickness.

Contribution

It provides a comprehensive analysis of optimal vortex growth across scales, highlighting the role of large-scale vortices and viscosity effects in shear flow turbulence.

Findings

Maximum vortex amplification occurs at scales near the disc height.

Large-scale vortices can reach high density perturbations near swing.

Growth of vortices is unaffected by bulk viscosity at large scales.

Abstract

Turbulent state of spectrally stable shear flows may be developed and sustained according to the bypass scenario of transition. If it works in non-magnetised boundless and homogeneous quasi-Keplerian flow, 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, $\lambda_y$, as compared to the flow scale-height, $H$. It is shown that with the account of the viscosity the highest possible amplification of shearing vortices, $G_{max}$, attains maximum at $\lambda_y\lesssim H$ and declines towards both the large scales $\lambda_y\gg H$ and the small scales $\lambda_y\ll H$ in a good agreement with analytical estimations based on balanced solutions. We pay main attention to the large-scale vortices $\lambda_y\gg H$, which produce $G_{max}\propto (\Omega/\kappa)^4$, where $\Omega$ and $\kappa$ 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 localisation 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.

Equations

$$\left(\frac{\partial }{\partial t}-q{\Omega }_{0}x\frac{\partial }{\partial y}\right)u_x-2{\Omega }_0{u}_y=-\frac{\partial W}{\partial x}+f_x+g_x,$$

$$\left(\frac{\partial }{\partial t}-q{\Omega }_{0}x\frac{\partial }{\partial y}\right)u_y+(2-q){\Omega }_0{u}_x=-\frac{\partial W}{\partial y}+f_y+g_y,$$

$$\left(\frac{\partial }{\partial t}-q{\Omega }_{0}x\frac{\partial }{\partial y}\right)W+a_{\ast }^2\left(\frac{\partial u_x}{\partial x}+\frac{\partial u_y}{\partial y}\right)=0,$$

$$f_{x,y}=\nu \left(\frac{{\partial }^2}{\partial x^2}+\frac{{\partial }^2}{\partial y^2}\right)u_{x,y},$$

$$g_{x,y}=({\nu }_b+\nu /3)\frac{\partial }{\partial x,\partial y}\left(\frac{\partial u_x}{\partial x}+\frac{\partial u_y}{\partial y}\right),$$

$$v_x\frac{\partial u_x}{\partial x}\sim \nu {ϵ}_u(r_0\lambda {)}^{-1},$$

$$\nu \frac{{\partial }^2{u}_x}{\partial x^2}\sim \nu {ϵ}_u(\lambda {)}^{-2},$$

$$x^{\prime }=x/H,y^{\prime }=(y+q{\Omega }_{0}xt)/H,t^{\prime }={\Omega }_{0}t.$$

$$\begin{gathered} \frac{d{\hat{u}}_x}{dt}=2{\hat{u}}_y-\mathrm{i}{\tilde{k}}_x\hat{W}-R^{-1}{k}^2{\hat{u}}_x- \\ (R^{-1}/3+R_b^{-1}){\tilde{k}}_x({\tilde{k}}_x{\hat{u}}_x+k_y{\hat{u}}_y), \end{gathered}$$

$$\begin{gathered} \frac{d{\hat{u}}_y}{dt}=-(2-q){\hat{u}}_x-\mathrm{i}{k}_y\hat{W}-R^{-1}{k}^2{\hat{u}}_y- \\ (R^{-1}/3+R_b^{-1})k_y({\tilde{k}}_x{\hat{u}}_x+k_y{\hat{u}}_y), \end{gathered}$$

$$\frac{d\hat{W}}{dt}=-\mathrm{i}({\tilde{k}}_x{\hat{u}}_x+k_y{\hat{u}}_y).$$

$$R\equiv a_{\ast }^2/(\nu {\Omega }_0),R_b\equiv a_{\ast }^2/({\nu }_b{\Omega }_0).$$

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

$$g(t)\equiv \frac{E_{ac}(t)}{E_{ac}(0)}.$$

$$(\mathbf{q},\tilde{\mathbf{q}})=\frac{1}{2}({\hat{u}}_x{\hat{\tilde{u}}}_x^{\ast }+{\hat{u}}_y{\hat{\tilde{u}}}_y^{\ast }+\hat{W}{\hat{\tilde{W}}}^{\ast }),$$

$$\frac{\partial \mathbf{q}}{\partial t}=\mathbf{A}\mathbf{q},$$

$$\frac{\partial \tilde{\mathbf{q}}}{\partial t}=-{\mathbf{A}}^{†}\tilde{\mathbf{q}}$$

$$(\tilde{\mathbf{q}},\mathbf{A}\mathbf{q})=({\mathbf{A}}^{†}\tilde{\mathbf{q}},\mathbf{q}).$$

$$\begin{gathered} \frac{d{\hat{\tilde{u}}}_x}{dt}=(2-q){\hat{\tilde{u}}}_y-\mathrm{i}{\tilde{k}}_x\hat{\tilde{W}}+R^{-1}{k}^2{\hat{\tilde{u}}}_x+ \\ (R^{-1}/3+R_b^{-1}){\tilde{k}}_x({\tilde{k}}_x{\hat{\tilde{u}}}_x+k_y{\hat{\tilde{u}}}_y), \end{gathered}$$

$$\begin{gathered} \frac{d{\hat{\tilde{u}}}_y}{dt}=-2{\hat{\tilde{u}}}_x-\mathrm{i}{k}_y\hat{\tilde{W}}+R^{-1}{k}^2{\hat{\tilde{u}}}_y+ \\ (R^{-1}/3+R_b^{-1})k_y({\tilde{k}}_x{\hat{\tilde{u}}}_x+k_y{\hat{\tilde{u}}}_y), \end{gathered}$$

$$\frac{d\hat{\tilde{W}}}{dt}=-\mathrm{i}({\tilde{k}}_x{\hat{\tilde{u}}}_x+k_y{\hat{\tilde{u}}}_y).$$

$$g_{opt}(t)\equiv \underset{\forall \mathrm{SFH}}{\max }g(t)$$

$$G_{max}\equiv \underset{\forall k_x,\forall t}{\max }{g}_{opt},$$

$$G(t)\equiv \underset{\forall k_x}{\max }{g}_{opt}$$

$$g_s\equiv g(t_s)$$

$$g_{kin}=\frac{E_{kin}(t)}{E_{kin}(0)},$$

$$R_T\sim 10^8\frac{(2-q{)}^3}{q},$$

$${\Phi }_{PW}=-\frac{M}{r-2M},$$

$$q\equiv \frac{r}{\Omega }\frac{d\Omega }{dr}=\frac{1}{2}+{\left(1-\frac{2M}{r}\right)}^{-1}.$$

$$\nu \sim v_p{l}_p\propto a_{\ast }/N\sim \Omega H^2/\Sigma \propto r^{1/2+\beta },$$

Full text

Key words— hydrodynamics — accretion, accretion discs — instabilities — turbulence — protoplanetary 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 pretty high value bearing in mind that the Reynolds number corresponding to the subcritical transition $R_{T}\sim 350$ for plane Couette flow and $R_{T}\sim 1000$ for plane Poiseille flow, see Schmid & Henningson (2001) as well as numerous references therein, and $R_{T}\sim 7000$ for the narrow gap rotating Taylor-Couette flow with angular velocity increasing towards the periphery, see Burin & Czarnocki (2012). However, there is a number of well-known observational evidences that cold astrophysical discs, where magnetorotational instability ceases to operate, are accreting. At the same time, the smallness of microscopic viscosity entails huge Reynolds numbers up to $\sim 10^{13}$ 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 both cannot take place without growth of perturbations, i.e. 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 perturbations amplitude (Henningson & Reddy, 1994). This basic fact is obtained from the Reynolds-Orr equation for evolution of kinetic energy of perturbations obeying Navier-Stockes equations. Thus, the energy transfer from the background motion to turbulent fluctuations is given by mechanisms which 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 to be 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 it 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 one is 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 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 well in literature (see, e.g., Reddy et al. (1998) and its citations): thanks to a possibility to simulate 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 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 coming into play stabilises 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 comes out of the shadow: 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 (look, e.g., Lesur & Longaretti (2005) for 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 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 wavevector space, which repopulates the quadrant with linearly growing shearing vortices. Lithwick (2007) has considered a particular case of such coupling between axisymmetric and shearing vortices in 2D model of Keplerian flow and showed 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, $\lambda_{y}$. Moreover, in later work Lithwick (2009) has found that non-linear vortices emerging by this 2D mechanism survive in 3D model provided that their azimuthal length-scales are larger than the flow scale-height, $\lambda_{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 non-linear sustenance of linearly growing shearing vortices. However, presently there is no theoretical understanding, whether the reservoir of linearly growing shearing vortices may be repopulated in full 3D dynamics of perturbations with not necessarily $\lambda_{y}>H$. Nevertheless, Meseguer (2002) has revealed a strong correlation between a value of maximum transient growth, $G_{max}$, and a point of relaminarisation 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, since $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 optimisation for solving the problem of perturbation dynamics. Yecko (2004) and Mukhopadhyay et al. (2005) (hereafter MAN05) have also identified the shearing vortices, as well as 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 case of the 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 one 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 reverse to the lift-up effect in plane shear flows, consequently, it was termed as the ’anti lift-up effect’ (Antkowiak & Brancher, 2007). The corresponding maximum optimal growth factor is several orders of magnitude larger than those for the Keplerian shear. Another important result is that the maximum optimal growth suffers 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 et al. (2014) confirm that in all spectrally stable regimes of Taylor-Couette flow, rather than only in quasi-Keplerian one, 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 making them different from quasi-Keplerian one. Unlike in quasi-Keplerian regime, in the cyclonic and counter-rotating regimes those shearing vortices have (though weak) axial dependence and up to 80% of kinetic energy transferred into the axial velocity perturbation near the transient growth maximum. This features go beyond swing amplification mechanism and Maretzke et al. (2014) note that the reason for such qualitative difference between quasi-Keplerian and other regimes remains unclear. Thereby, two-dimensionality of shearing vortices attaining maximum transient growth is a unique peculiarity of the quasi-Keplerian flow, which may prove to be the reason for its extraordinary non-linear HD stability.

The studies mentioned above have been dealing with incompressible perturbation dynamics, which in the context of accretion disc physics implies that perturbations have characteristic length-scale much smaller than the disc thickness. We will refer to such perturbations as ’small-scale vortices’ below in the text. In contrast, we are going to 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 will refer to as ’large-scale vortices’ below in the text, should strongly depend on the shear rate as one approaches constant angular momentum rotation, and $G_{max}\propto\kappa^{-4}$, where $\kappa$ is epicyclic frequency expressed in units of rotation frequency. In this work we are going to check numerically this analytical result employing a variational technique for 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, what is especially important to do as the spiral swings from leading to trailing one, since as this happens, the radial wavenumber of SFH tends to zero and one cannot distinguish between vortices and density waves anymore. Thereby, the field of the velocity perturbation acquires significant divergence, which may cause a suppression of transient growth of initially vortical SFH. Also, as far as the azimuthal wavenumber of SFH is of 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, characterised by the bulk viscosity, see e.g. Landau & Lifshitz (1987).

Whereas in this work such kind of problem is considered for the first time in application to rotating and, particularly, quasi-Keplerian flows, Hanifi et al. (1996), Farrell & Ioannou (2000) and Malik et al. (2006) have already studied it with an interest to compressible plane shear flows. They found that the lift-up effect remains to be valid with the account of finite Mach number, and that the streamwise rolls remain to be the optimal perturbations in the flow. With regards to the shearing vortices in compressible medium, Farrell & Ioannou (2000) have found their growth factors being considerably larger than corresponding incompressible values, which is explained by 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 organised as follows. We start describing the 2D model for local perturbations of quasi-Keplerian shear flow in §2. Next, in §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 find that shearing vortices are almost unaffected by the value of bulk viscosity, see §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 §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 a good agreement with analytical estimations, relegated to Appendix B. In §4.3 we discuss possible consequences of our results speculating that the transition Reynolds number, as well as the dimensionless azimuthal turbulent stress, may strongly depend on the shear rate. Additionally, in §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 localised leading spiral having a certain radial position and a degree of winding both altering with the optimisation timespan. 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 Conclusions.

2 Model for dynamics of perturbations

In order to study the transient dynamics of perturbations taking into account a finite disc scale-height, 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. That 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 velocity field independent of vertical coordinate preserve hydrostatic equilibrium, i.e. they have planar velocity field. Thus, the 3D problem with such a restriction is reduced to 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 §6.1.

Substituting the polytropic relation between Eulerian perturbations of pressure and density, $p_{1}$ and $\rho_{1}$, as well as their background counterparts, $p$ and $\rho$, into the continuity equation, we change to a new variable $W=p_{1}/\rho$ 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 get:

[TABLE]

[TABLE]

[TABLE]

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 $\Omega_{0}$ is angular velocity of the fluid rotation at $x=0$ corresponding to some radial distance $r_{0}\gg x$. Also, $q$ is a constant shear rate that defines the background azimuthal velocity as $v_{y}=-q\Omega_{0}x$. All other background quantities, such as $\Sigma\equiv\int\rho dz$ and $a_{*}$, are assumed to be constant. Note that $a_{*}^{2}=na_{eq}^{2}/(n+1/2)$ with $a_{eq}$ being sound speed in the equatorial plane of disc. The inviscid part of (1-3) coincides with the set (28-30) given in ZR14. At the same time, the new ingredients are viscous terms, which we take as follows:

[TABLE]

[TABLE]

where $\nu$ and $\nu_{b}$ denote kinematic viscosity coefficients corresponding, respectively, to the shear and bulk viscosities of the background flow, $\eta$ and $\zeta$, integrated across the disc thickness and subsequently normalised by $\Sigma$.

In order to obtain $f_{x,y}$ and $g_{x,y}$ in the form (4-5), we first collect perturbed parts of all dissipative terms entering the Navier-Stockes equations, see for example Appendix B in Kato et al. (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 $\lambda\ll r_{0}$, only the dissipative terms containing the highest order spatial derivatives of perturbations should be retained. This suggests that all contributions containing $\nu\rho_{1}/\rho$ or $\eta_{1}/\rho$, where $\eta_{1}$ is the Eulerian perturbation of dynamical viscosity, or proportional to the radial background velocity $v_{x}\sim\nu/r_{0}$ should be omitted. Indeed, for example the advection term

[TABLE]

whereas one of the terms from (4)

[TABLE]

which is $r_{0}/\lambda$ times larger ($\epsilon_{u}$ is an amplitude of $u_{x}$). In this way, we obtain expressions (4-5).

2.1 Dimensionless equations for single SFH

Shear flow with $a_{*}<\infty$ has a finite scale-height $H=a_{*}/\Omega_{0}$ what follows from the condition of vertical hydrostatic equilibrium in unperturbed disc with a polytropic equation of state. Following ZR14, let us use the natural $H$ in our choice of dimensionless comoving Cartesian coordinates:

[TABLE]

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})$111it is assumed everywhere below that $k_{y}>0$. are the dimensionless wavenumbers along $x$ and $y$ axes expressed in units of $H^{-1}$. Omitting the primes from now on and introducing the shearing radial wavenumber $\tilde{k}_{x}\equiv k_{x}+qk_{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:

[TABLE]

[TABLE]

[TABLE]

In equations (7-9) it is assumed that $\hat{u}_{x},\hat{u}_{y}$ are expressed in units of $a_{*}$ and $\hat{W}$ is expressed in units of $a_{*}^{2}$. Also, there are two different Reynolds numbers in (7-9), $R$ and $R_{b}$, which parametrise dynamical and bulk viscosities, respectively. Explicitly,

[TABLE]

In Appendix A.1 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^{-1}$. Also, the dynamical problem becomes of the first order, i.e. for particular values of $k_{x},k_{y}$, $q$ and $R$, one uniquely determines a normalised solution and the corresponding transient growth factor, see equations (67) and (69). On the contrary, the problem (7-9) has the third order, that is why, one first has to find a particular solution of (7-9) corresponding to the highest transient growth factor for a specified time interval. This solution is referred to as the optimal one. To find this optimal solution, in §3 we are going to supply the set of equations (7-9) with the corresponding set of adjoint equations. After that we will be able to implement the iterative loop elucidated in §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:

[TABLE]

The growth factor of an arbitrary SFH is then given by

[TABLE]

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

[TABLE]

where asterisk stands for complex conjugation and ${\bf q}\equiv\{\hat{u}_{x},\hat{u}_{y},\hat{W}\}$, $\tilde{\bf q}\equiv\{\hat{\tilde{u}}_{x},\hat{\tilde{u}}_{y},\hat{\tilde{W}}\}$.

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

[TABLE]

and the subsequent integration of adjoint equations

[TABLE]

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

[TABLE]

The explicit form of ${\bf A}$ and ${\bf A}^{\dagger}$ is defined by the particular model for perturbations and particular choice of their norm. In the case when ${\bf A}$ is represented by equations (7-9) and the inner product is given by (13), operator $\mathbf{A}^{{\dagger}}$ corresponds to the following set of adjoint equations:

[TABLE]

[TABLE]

[TABLE]

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 $\lambda\sim r_{0}$, see §5.3.

The set (17-19) looks rather similar to the set (7-9) except that the coefficients standing 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 an opposite signs in the adjoint equations. The latter implies that (17-18) describe viscous diffusion with a negative diffusion coefficient. However, since the adjoint equations are advanced backwards in time during the implementation of iterative loop, one does not face any sort of instability caused by a non-zero negative viscosity.

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

integrate equations (7-9) forward in time up to the instant $t>0$;

• 2)

substitute the resulting ${\bf q}(t)$ into equations (17-19) and integrate them backwards in time up to the instant $t=0$;

• 3)

renormalise the upgraded ${\bf q}(0)$ to the unit norm;

• 4)

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

As it is discussed in detail by ZR14, the final ${\bf 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

[TABLE]

up to the time $t$ among all other SFH’s. 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

[TABLE]

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

Additionally, we are going to look at an intermediate quantity

[TABLE]

in the case of compressible perturbations.

In what follows, $G_{max}$ is determined using common numerical routines for search of extremum. In the rest of the paper we are going to present a parametric study of $G_{max}$ depending on $k_{y}$, $R$ and $q$222We also show the growth factors of corresponding optimals at the instant of swing in order to extinguish the exited density waves, see §4.2.. Such a kind of nonmodal 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 $\lambda_{y}\gtrsim H$ along with the case $\lambda_{y}\ll 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 work, see equation (73).

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 then 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 optimisation procedure for determining $g_{opt}$ is equivalent to looking for a pure vortex. At the time $t_{s}$ called the instant of swing the leading spiral turns into a trailing one starting its decay. By definition, $\tilde{k}_{x}(t_{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}\sim 1$, see Heinemann & Papaloizou (2009). On the one hand, as soon as $k_{y}\sim 1$, they significantly enhance the optimal growth $G(t)$ as defined by equation (22), but on the other hand, it is not clear what role they play in production and maintenance of turbulence. Additionally, as density waves produce highly divergent velocity field, they must be strongly affected by value of bulk viscosity which is poorly known in astrophysical flows. For these reasons, we are going to focus on the dynamics of vortices only. However, in order to be confident that a reliable quantitative criterion will be chosen to conclude about the transient amplification of vortices irrespective to excitation of density waves, one should first examine the optimal solutions to general problem represented by equations (7-9).

In this regard, we present the result of optimisation 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_{05}=2000$. Note that $R_{05}$ is related to $R$ by equation (65). On the top panel in Fig. 1 one can see the highest possible enhancement of linear perturbations in the Keplerian shear which is plotted vs. $k_{y}$. First, we start with the case of zero bulk viscosity demonstrated by the solid curve and find the following:

despite the fact we consider boundless shear, it does have a maximum at $k_{y}\sim 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 (71) and (73) derived in Appendix B employing the so called balanced solutions, see e.g. (35-37) by 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$ one finds a bump which apparently corresponds to the effect of density waves generation by vortices and is referred to as the ’sonic bump’ hereafter.

The point (4) becomes more clear as we look at the bottom panel in 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}=\infty$ and obtain the profile of $G(t)$, see solid curve. In turns out that while approaching the area of $k_{y}\sim 1$, $G(t)$ acquires an additional local maximum which is located farther in time than $t_{max}$ given by equation (72). So, $t_{max}\approx 19$ for the parameters taken to obtain the solid curve, which matches the location of its left maximum. If it were not for prominent excitation of density waves, it would be a unique extremum of $G(t)$, see Fig. 9 by MAN05. Let us denote it by $G_{max}^{\prime}$. However, as it is shown by ZR14 in the inviscid model, see their Fig. 2, the vortices with $k_{y}\sim 1$ that have already swung well in advance of some fixed time interval $t_{0}$ produce an additional bump in the profile of $g_{opt}(k_{x},t=t_{0})$ 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_{y}t_{0}$ and produced by vortices swinging right at $t=t_{0}$. As a consequence, $G(t)$ acquires the second local maximum, which we denote by $G_{max}^{\prime\prime}$. Further, as we approach $k_{y}\approx 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 on the top panel in Fig. 1, where $G_{max}=G_{max}^{\prime\prime}$.

Now, let us consider the non-zero bulk viscosity. On both panels in 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}=\infty$ 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 on the top panel in Fig. 1, which represents

[TABLE]

of vortical SFH that yields $G_{max}(k_{y})$ for $R_{b}=\infty$. Clearly, in the limit $R_{b}\to 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}=\infty$ 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}\approx 16$. Also note that in the case of $R_{b}=\infty$ the oscillations of $\hat{u}_{x}$ start damping approximately at $t\approx 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 it is checked that $t_{max}$ numerically determined for wings of $G_{max}(k_{y})$ corresponding to $k_{y}\ll 1$ ($k_{y}\gg 1$) fits the analytical estimation (72) very well.

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. That is why the numerical analysis presented above appears to be necessary to check the robustness of analytical conclusions.

Thus, in the next section we are going to 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 scale-height

As it is 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, oppositely, the streamwise streaks which generate streamwise rolls what 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 $\propto R^{2}$ and attains several thousands for $R_{05}=2000$, see MAN05. In other words, this type of transient growth is produced by a nearly axisymmetric ($k_{y}\approx 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 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}\approx 2/(2-q)$, which gives, for example, $G_{max}\sim 20$ at $q=1.9$. At the same time, the swing amplification of columnar (independent on $z$) small-scale vortices $\lambda_{y}\ll H$ do not produce a substantial growth before they are damped by the viscosity, so that the corresponding maximum growth has value $G_{max}\sim 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 $\lambda_{y}\gtrsim H$ proves to be significantly larger than that of vortices with $\lambda\ll H$ in the Keplerian case $q=1.5$. Additionally, according to the analytical result (73), 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 stabilisation of fluid elements by Coriolis and tidal forces becomes weaker. For this reason, we present results of numerical optimisation for $1.5<q<2.0$ involving the compressible dynamics in this Section, see Fig. 3. There are solid curves that demonstrate gradual increase of growth factor as we approach $q=2$. As we see, at $q=1.9$ $g_{s}\sim 10^{3}$ 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 (73), see the dotted lines in Fig. 3. For a confidence, we numerically check that the dependence $g_{s}\propto 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 scale-height reduces much more gradually as one introduces a positive radial gradient of the specific angular momentum.

Besides, $g_{s}$ diverges as $k_{y}\to 0$ in the constant specific angular momentum flow, $q=2$. In other words, it obeys the same law as in the incompressible limit: $g_{s}\propto k_{y}^{-4/3}$. Using equation (57) from ZR14 at $\kappa=0$ we find that inviscid balanced solutions (35-37) of RZ15 yield $G_{max}=t_{max}^{2}$ in the limit $k_{y}\ll 1$ and $|k_{x}|\gg 1$, which also recovers the incompressible dependence on $k_{y}$ after we substitute $t_{max}$ (72). However, the numerical curve for $q=2$ displayed in Fig. 3 gives a somewhat larger value in the limit $k_{y}\ll 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}\to 0$ which enhances the growth factor of optimal vortex comparing to the value coming from balanced solution at $t=t_{s}$. This claim is supported by the analytical result of Heinemann & Papaloizou (2009) which states that the amplitude of generated density wave grows exponentially as the parameter $\epsilon\ll 1$ increases, see e.g. equation (34) of RZ15. As long as $\kappa\sim 1$ and $k_{y}$ is not too close to unity, $\epsilon$ is small, the theory constructed by Heinemann & Papaloizou (2009) is valid and the generated wave has a small amplitude at $t=t_{s}$. That 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 $\kappa\to 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}\lesssim 1$ suffers strong excitation of density waves, making the use of the balanced solution unjustified333This 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 $\epsilon$ 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}\ll R$, even for $q=2$ in the region $\epsilon\sim 1$. Though, still it 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

[TABLE]

where $E_{kin}=1/2(|\hat{u}_{x}|^{2}+|\hat{u}_{y}|^{2})$, see the dashed curves in Fig. 3. As we see, $g_{ks}\equiv g_{kin}(t_{s})$ of the optimal large-scale vortices ceases much faster than $g_{s}$ as one shifts to $k_{y}\to 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 be seen also 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 larger by more than 6 times 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 two only. But still the kinetic energy is only a few percents of the acoustic energy of optimal SFH at the peak of its norm in the particular case introduced in Fig. 4444Additionally, we have checked that change to norm of perturbations with term $|\hat{W}|^{2}$ in equation (11) multiplied by a constant $p\in(0,1]$ and corresponding change of the adjoint equations leads to the same optimal solutions being the leading spirals of the vortical SFH with negative $k_{x}$..

We are not the first who 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 optimal perturbations are independent of the streamwise direction similarly to the incompressible case, found also that these optimals contain large amplitude of temperature perturbation, see their Fig. 12. Next, Farrell & Ioannou (2000) solved a 2D optimisation problem in a hypersonic plane shear flow, which is more analogous to our model. Similarly to our results, they obtained 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 density waves excitation at the instant of swing. For the particular situation they considered 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, since balanced solutions acquire the non-zero velocity divergence and non-zero density perturbation with an account of compressibility, see equations (35-37) by RZ15. It should be kept in mind, however, that the balanced solutions (as well as their accurate numerical counterparts) are indeed vortices since they are proportional to $I_{\nu}$ and are not affected by bulk viscosity, see §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. On the contrary, this is not the case in super-Keplerian rotating flow with $q\to 2$ since we find the enhancement of transient growth of large-scale vortices irrespectively 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 since 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 set of equations allows the existence of an axisymmetric vortical solutions for $q\lesssim 2$, which are akin of streamwise streaks with a non-modal growth strongly modified by epicyclic motions and slightly modified by divergent motions up to $\lambda\sim H$.

Finally, we argue that it is questionable which way the high amplitude density perturbations emerged from the transient growth of vortices with $\lambda_{y}\gtrsim H$ should affect the stage of non-linear interaction of different SFHs. Unlike density perturbations of the excited density waves, they hardly can lead to generation of shocks, because at the peak of their amplitudes, i.e. near the vortex swing, they have the least spatial gradient and, oppositely, 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 density wave grows simultaneously with its spatial gradient, what 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}{\bf u}$ entering the continuity equation for finite amplitude perturbations. However, to the best of our knowledge, it is not clear, whether this has any consequences with regards to the problem of positive non-linear feedback in the bypass transition to turbulence. To date, a 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) . Further, theoretical studies has been carried out with emphasise 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 transient growth factor is contained in enthalpy perturbation, the kinetic part, $g_{ks}$, remains to be significant for optimal SFH with $k_{y}\lesssim 1$ even in the Keplerian case. Not less importantly, $g_{ks}$ gets larger as $q\to 2$ even more steeply than $g_{s}$ attaining the value $g_{kin}\sim 600$ for $q=1.9$ due to the substantial increase of $\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 subcritical stability boundary in a viscous counter-rotating spectrally stable Taylor-Couette flow. Recently, Maretzke et al. (2014) suggested that the threshold for transition in any regime of Taylor-Couette rotation should be $G_{max}\sim 100$ basing on a variety of new experimental results and nonmodal growth calculations. Of course, it is not clear yet, whether this is true for the quasi-Keplerian regime, since to date, nobody has detected its turbulent state. Even more, it is not known, whether the transition threshold for $G_{max}$ exists in the case of vortices in 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})\sim const$, where $R_{T}$ is the Reynolds number of the subcritical transition in quasi-Keplerian case $1.5\leq q<2$. If so, equation (73) yields

[TABLE]

where we set $R_{T}\sim 10^{7}$ for the Keplerian rotation which is stable at least up to the Reynolds number equal to a few millions, see Edlund & Ji (2014).

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

[TABLE]

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 that

[TABLE]

Equation (26) is plotted in Fig. 5 which shows that for $r<20M$, where relativistic disc emits most its energy, $q>1.6$ and $R_{T}\lesssim 4\cdot 10^{6}$.

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

[TABLE]

where $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 disc has a constant aspect ratio $H/r$ and $\Sigma\propto r^{-\beta}$ with $\beta>0$555In equation (27) it is assumed that cross section of collisions between the particles is independent of $a_{*}$ and $N$.. Consequently, the situation is possible when relativistic disc is turbulent inside some radius $r=r_{T}$, whereas turbulence ceases in the outer regions of disc. Presumably, if this happens, the angular momentum flux is absent throughout the laminar flow, and the inner disc exhibits a non-stationary accretion. Whether this process may become recurrent, is unclear and requires an additional research.

Furthermore, it is known that in a disc with 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, actual $q(r)$ should become even steeper comparing to the curve shown in Fig. 5. Of course, the enhanced transient growth of large-scale vortices may drive turbulisation 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 scale-height 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 a smooth background shear they are embedded in. Those non-linear vortices have been found in a large number of simulations including mentioned in 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 dics, see e.g. Heng & Kenyon (2010).

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

[TABLE]

where $\alpha=\alpha_{K}$ is its Keplerian value. Local and global MHD simulations of relativistic accretion discs, see e.g. Penna et al. (2013) and Morales Teixeira et al. (2014), show a radially varying $\alpha$ 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) by Abramowicz (2009) and in the abstract of Penna et al. (2013). Since $q$ increases as we approach the inner edge of 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 $\alpha$ on $r$ comparing to what has been obtained to date with regards to supercritical magnetic turbulence. In the case of relativistic disc, such a substantial radial variation of $\alpha$ 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), $\alpha$ diverges in the vicinity of the inner edge of standard relativistic disc. However, the finite aspect ratio of disc, $\delta=H/r_{0}$, puts a lower limit on the value of $k_{y}\sim\delta$, what constrains the value of optimal growth as $q\to 2$, see Fig. 3. Additionally, the non-zero pressure radial gradient close to the inner edge of standard disc always results in the shear rate $q<2$.

5 Extension to global spatial scale

So far, we have been considering transient dynamics of large-scale shearing vortices in the local framework implying that $r_{0}\gg\lambda_{y}\gtrsim h$. However, as $k_{y}\sim\delta$, one has to take into account a global disc structure and the cylindrical geometry of the problem. At first, disc generally is radially non-uniform as its surface density and its microscopic viscosity vary on the scale $\sim r_{0}$. Next, the transient amplification of large-scale vortices substantially increases along with the local shear rate, which in turn may be a function of the radial distance. For example, in 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 may 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 overall conclusion that regardless of a particular radial structure of disc it is always able to generate prominent transient growth of global perturbations at huge Reynolds numbers realised 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

Geometrically thin disc in the Paczynsky-Wiita potential (25) rotates with the angular frequency

[TABLE]

and is subject to radial oscillations with the following epicyclic frequency

[TABLE]

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

[TABLE]

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

[TABLE]

Next, one has to specify the kinematic viscosity coefficients, $\nu$ and $\nu_{b}$, across the flow. Similarly to what has been defined in §2, those are the shear and the bulk viscosities, $\eta$ and $\zeta$, integrated over the disc thickness and subsequently divided by surface mass density, $\Sigma$. As we have mentioned above, see equation (27), it is reasonable to assume that $\nu\propto H^{2}\Omega/\Sigma$. Provided that the value of $\nu$ is known at some distance $r=r_{1}$,

[TABLE]

where the variables marked by the subscript ’1’ stand for the corresponding quantities taken at $r_{1}$. Let us choose $r_{1}$ being the distance, where $\nu(r)$ attains the minimum. A value of $\nu_{1}$ can be specified via the Reynolds number $R$ (cf. equation (10)):

[TABLE]

Combining (32) with (33) we obtain the final expression for the kinematic viscosity coefficient profile:

[TABLE]

where it is taken into account that $R_{05}=q(r_{i})R$, see §A.1, 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

[TABLE]

At the next stage we are going to get an expressions 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$

[TABLE]

and mass conservation in the form

[TABLE]

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

Equation (37) allows us to integrate (36) assuming that $r\varphi$-component of the stress tensor vanishes at the inner boundary of disc $r=r_{i}$:

[TABLE]

In equation (38) we introduce $D=1-\left(\Omega_{i}r_{i}^{2}\right)/\left(\Omega r^{2}\right)$. Note that for the Newtonian potential $D=1-(r_{i}/r)^{1/2}$, whereas for the Paczynsky-Wiita potential $D=1-\left(r_{i}/r\right)^{3/2}\left(r-2M\right)/\left(r_{i}-2M\right)$.

Combining (37) with (38) we obtain the radial velocity

[TABLE]

Then, combining (32) with (38) allows us to obtain the radial profile of $H$:

[TABLE]

In order to parametrise $H$, we use the disc aspect ratio $\delta\equiv(H/r)|_{r=r_{2}}$ with $r_{2}$ being the distance, where $H/r$ takes its maximum value. Finally,

[TABLE]

where the variables marked by the subscript ’2’ denote the corresponding quantities taken at $r_{2}$.

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

[TABLE]

From the other hand, the diffusive radiation transfer in vertical direction yields an approximate relation between the surface and the equatorial temperatures, $T_{eff}$ and $T_{eq}$:

[TABLE]

where $l_{ph}\propto\rho^{-1}$ 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, one concludes that

[TABLE]

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

[TABLE]

where $C_{\Sigma}$ is a dimensional constant. Since 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. Everywhere below we denote it as ’P1’ in the case of $\Omega$ and $\kappa$ given by equations (29) and (30), otherwise, i.e. in the limit $r\gg 1$ when $\Omega$ and $\kappa$ 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 which can be achieved by the inner region of an arbitrary rotating laminar flow with $\delta\ll 1$ evolving due to microscopic viscous relaxation only. However, someone would argue that this stage is quite unrealistic, since microscopic viscosity is too small in typical astrophysical situation leading to a very long relaxation time which most likely exceeds the turbulisation 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), (30), (31). However, now the surface density is assumed to be a constant

[TABLE]

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

[TABLE]

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

[TABLE]

Note that in a disc of such 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) below. At last, the radial velocity is set to zero

[TABLE]

Thereby, equations (29), (30), (31), (45), (35), (44), (46) and (47) specify our homogeneous configuration. Similarly, everywhere below we denote this configuration as ’P2’ in the case of $\Omega$ and $\kappa$ 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 linearise the Navier – Stokes equations written in the cylindrical coordinates $\{r,\varphi,z\}$ (see e.g. Appendix B of Kato et al. (2008)). Following §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 considering the planar velocity perturbations independent on vertical coordinate. Thus, perturbations are described by three variables: the Eulerian perturbations of velocity components, $\delta v_{r}$ and $\delta v_{\varphi}$, and the Eulerian perturbation of enthalpy $\delta h$. Owing to the rotational symmetry of background, it is sufficient to consider complex amplitudes of the Fourier harmonic $\propto e^{{\rm i}m\varphi}$ 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, whereas the viscous part of equations contains terms emerging due to the non-zero perturbation of shear, which we denote by $N_{r,\varphi}$, terms emerging due to the non-zero perturbation of the velocity divergence, which we denote by $B_{r,\varphi}$, and, additionally, the advective terms $\propto 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,

[TABLE]

[TABLE]

[TABLE]

where the viscous terms $N_{i}$ and $B_{i}$ are shown in Appendix D.1. In equations (48-50) we have also assumed that the Eulerian perturbations of dynamical and bulk viscosities are equal to zero $\delta\eta=\delta\zeta=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 $\mathbf{A}^{{\dagger}}$ is needed. In the global problem $\mathbf{A}^{\dagger}$ is defined by the Lagrange rule (see equation (16) ) through the inner product:

[TABLE]

corresponding to a norm equal to total acoustic energy of perturbations

[TABLE]

Derivation of the adjoint equations is described in detail in §3.4 RZ15 for the inviscid case. However, it can be straightforwardly generalised onto viscous problem. We obtain:

[TABLE]

[TABLE]

[TABLE]

where the viscous terms $\tilde{N}_{i}$ and $\tilde{B}_{i}$ are shown in Appendix D.1.

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

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 2nd order method employs 4 uniform grids in the space $(r,t)$ shifted to each other, which allows to use central differences for the approximation of spatial derivatives of the unknown variables. However, the approximation of the advection terms (that $\propto v_{r}$) by central differences leads to numerical instability. In order to check this, one should apply Von Neumann stability analysis method, see Charney et al. (1950). That 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 $\partial\delta\tilde{h}/\partial t$ in the grid node neighbouring to the boundary, the value of $\delta\tilde{h}$ at the boundary is required since 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 $\delta h$ is not required.

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

5.5 Results

As a result of the global optimisation problem we are going to obtain curves of optimal growth

[TABLE]

where $E_{a}$ is given by equation (52). Also, it is a common practice to consider another quantity

[TABLE]

which is similar to $G_{max}$ introduced by MAN05 as well as Yecko (2004) and Maretzke et al. (2014), also cf. equation (21) as its local analogue. Note that $G_{max}$ depends on the parameters of disc such as $\delta$, $R_{05}$ and $R_{b}$ and on 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 waves excitation during the swing of a global spiral of perturbations (see §4.1 for details about the corresponding effect in the local framework), we use high value of bulk viscosity $R_{b}=0.001R_{05}$ in all calculations. It is implied that time is expressed in units of $\Omega^{-1}(r_{i})$, the frequencies are expressed in units of $\Omega(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}$

[TABLE]

and a global analogue of $R$

[TABLE]

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 (58) and (59).

5.5.1 Models N1 and N2

Let us consider the optimal growth of global perturbations in the case of Newtonian background. We start 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, cf. bottom panel in Fig. 1, though for some reason $G$ tends to a non-zero value as $t\to\infty$. 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 local optimal growth curve obtained previously, e.g., by MAN05. By contrast, it can be seen that as $m\to\infty$ 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 (72), evaluated near the inner boundary of disc. Equation (72) gives $t_{max}\approx 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 outskirts 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 $\Omega^{-1}$ at the new $r=r_{sp}$ where the spiral has shifted. The time interval until the instant of swing, $t_{max}\Omega^{-1}(r_{sp})$, becomes longer and longer and it is always approximately equal to the optimisation timespan, $t$. In fact, the latter is a general condition for the approximate localisation of global spiral:

[TABLE]

At the same time, the global optimal spiral with $m\gg 1$ corresponding to $\bar{k}_{y}$ should recover the local value of $G$ in whatever part of disc it was localised. But from the above it follows that this local value of $G$ must be equal to nothing but $G_{max}$. Then, since $\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\gtrsim t_{max}\Omega^{-1}(r_{i})$ which is close to the analytical value of $G_{max}$ given by equation (73) with the values of $\bar{k}_{y}$ and $\bar{R}$ substituted in there, see the dot-dashed curve on top left panel in Fig. 7.

The picture described above indicates that the hump in the curve of $G(t)$ and the related suppress of plateau for perturbations $m\sim 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 optimisation time: compare the initial optimal profiles corresponding to square and circle on the top left panel in Fig. 7. As has been discussed, such an additional upwind vanishes when $m\to\infty$. So, this purely global effect may be probably ascribed to a non-zero vorticity gradient since the most simple variant of background configuration corresponding to $\bar{k}_{y}=const$ and $\bar{R}=const$ has been considered herein. However, the clear explanation remains to be done. Anyhow, as Fig. 7 illustrates, the optimal spirals with $m=1$ also upwind/shift in a good agreement with the general picture that has been suggested and expressed through equation (60) just above: for example, equation (60) yields $r_{sp}\approx 3$ and $r_{sp}\approx 4$, respectively, for dashed and dot-dashed curves on bottom panel in Fig. 7.

Another purely global effect that is absent in the local formulation of the problem is modal growth of perturbations with $m\sim 1$ in the range of sufficiently small $\delta\lesssim 0.05$. In Fig. 8 the dependence $G(T)$ and corresponding radial profiles of optimal initial perturbations are shown for $\delta=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\to\infty$. Obviously, at large $t$ the optimal growth recovers the growth factor of the most unstable global mode. The amplitude of this mode increases $\propto\exp(\gamma t)$, with $\gamma\approx 0.001$ in Fig. 8.

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 in Fig. 8, the optimal initial profile has nothing in common with leading spirals and converges to the unique profile being the shape of the most unstable mode. We have checked that it corresponds to a pattern rigidly rotating with angular velocity $\omega_{m}\approx 0.98$, i.e. with corotation radius located inside of the flow. It is beyond the scope of this work to study the physical nature of the revealed instability, though presumably it is akin of the well-known Papaloizou-Pringle instability, see Papaloizou & Pringle (1987), which has been also studied in details 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 localised at corotation radius. The latter mechanism of instability is usually called the Landau mechanism, see the book by Fabrikant et al. (1998) for the 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 comparing to model N2. The plateau has disappeared and instead we obtain a curve with a global maximum and $G\to 0$ as $t\to\infty$ 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-shape 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. Curves on the bottom panel in Fig. 9 demonstrate that as the optimisation timespan increases, the initial optimal spiral shifts outskirts and unwinds. This occurs because for $r>r_{1}$, see equation (32), the viscosity grows to the periphery of disc, so that $\bar{R}$ decreases to the periphery of disc for $r>r_{1}$. Consequently, the dimensionless $t_{max}$ expressed in units of $\Omega^{-1}(r)$ decreases outskirts as well. Further, as we have suggested above, the optimisation time is equal to the instant of swing of the optimal spiral which, in turn, is approximately equal to $t_{max}$. Since 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 $\Omega^{-1}(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 $\Omega^{-1}(r_{i})$, which is given by equation (60), remains to be a quantity which grows outwards, otherwise, the optimal spiral would not shifted at larger $r_{sp}$ for larger optimisation timespans.

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_{sp}),\bar{R}(r_{sp}),q(r_{sp})\,)$. This procedure allows us to obtain the dot-dashed curve on top left panel in Fig. 9 using the definitions (58) and (59) together with specifications of the model given in §5.1.1. As we see, the analytical curve obtained by means of simple estimations of local optimal growth is in a reasonable agreement with the numerical one obtained 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}\approx 7$ and $r_{sp}\approx 12$, respectively, for dashed and dot-dashed curves on the bottom panel in Fig. 9.

We have numerically checked, that the modal growth of perturbations with $m=1$ does not exist for any value of $\delta$ in the case of model N1. Presumably, unstable mode which is located near the inner boundary of disc, see Fig. 8, is suppressed by a sufficiently large viscosity therein, since in the model N1 $\bar{R}\to 0$ as $r\to 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 hump in 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. $\lambda_{y}>h$) stays almost at the same level as one passes from local perturbations to global perturbations with $m\sim 1$, whereas small-scale vortices (i.e. $\lambda_{y}<h$) with $m\sim 1$ reach substantially smaller growth factors than their local counterparts.

At last, let us make an important note about coexistence of modal and transient growth of perturbations with $m\sim 1$ in the sufficiently thin discs specified by the model N2. In the Fig. 7-10 we fix $R_{05}=2000$ as it is done in the local problem for easier comparison with the results of MAN05. However, one should keep in mind that actual Reynolds numbers in astrophysical discs are orders of magnitude larger and approach values like $\sim 10^{10}\div 10^{13}$. This implies that in practice the transient growth dominates the modal growth by many orders of magnitude for any value of $\delta$, since the optimal growth $\propto R_{05}^{2/3}$ whereas the increment increases by a few percent only as $R_{05}\to\infty$.

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, $\kappa(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\to 2$. As represented in Fig. 3, for the case $k_{y}=0.2$ the local transient growth attains the value $g(t_{s})\sim 10^{3}$, whereas at the same time $G_{max}\approx 2250$ for $R_{05}=2000$ and $R_{b}\ll R_{05}$, see Fig. 14 below. At first, we check the transient growth of global spirals with $m=1$ in the model of disc with uniform surface density which has been denoted as P2, see Fig. 11. Following the strategy adopted to 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 timespans while $t<t_{max}$, where $t_{max}\approx 20$ is estimated at $r_{i}$. After that, the spiral starts shifting outskirts which causes a decline of optimal growth. Now this decline of $G$ is much stronger than in the Newtonian case, cf. Fig. 7, since the value of optimal growth of large-scale vortices falls down as one approaches the Keplerian shear rate. On the top left panel in Fig. 11 it is shown that for $t\gtrsim 20$ the analytical curve is in a moderate agreement with the numerical curve of $G(t)$ generated for global perturbations with $m=10$, and the agreement becomes excellent when $t\to\infty$. The analytical curve corresponds to $G_{max}(\,\bar{k}_{y}=const,\bar{R}=const,q(r_{sp})\,)$ given by equation (73), where we substitute $r_{sp}$ evaluated from equation (60) into equation (26). $G_{max}$ is not constant, indeed, due to the dependence $q$ on $r$ in P2. Also, since the analytical profile defined in this way diverges as $r_{sp}\to r_{i}$ we truncate it for $t\lesssim 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 hump in $G(t)$ corresponding to $m=10$ is in a good agreement with the latest value. Additionally, as $t\to\infty$ the optimal growth approaches its Newtonian values for both $m=1$ and $m=10$, cf. figs. 7 and 11.

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

In the case of model P2, the modal growth emerges for $\delta<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 in Fig. 12. Then, for a sufficienly long $t$, we find that the optimal initial profile of perturbations converges to a rigidly rotating pattern, see bottom panel in Fig. 12. We check that its increment has the maximum $\gamma\approx 0.03$ located at $\delta\approx 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, $\Re[\delta h(t)]$ at some radius. We find that it corresponds to the corotation radius at $r_{c}\approx 1.04$, which is close to the inner boundary of disc. Yet, the critical layer of 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_{05}=2000$ the instability exists in the substantial range of $\delta$ and the unstable mode has a noticeable increment, as one proceeds to larger Reynolds numbers typical for astrophysical discs, the transient growth starts dominating the modal growth at any $\delta$ for the same reasons as in the case of model N2.

Finally, the numerical consideration of global perturbations reveals that at $R_{05}=2000$ the disc specified by the model P1 is linearly stable for all $\delta$, 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 below the 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 somewhat farther from the inner boundary of disc than a similar spiral in P2, see the dashed curve on the top right panel in Fig. 11 and the solid curve on the bottom panel in Fig. 13. In the former situation the spiral is centered at $r_{sp}\approx 1.5$, where $q\approx 1.8$, while in the latter situation spiral is centered at $r_{sp}\approx 3$, where $q\approx 1.6$. For the local large-scale vortices such a difference in shear rate corresponds to drop in transient growth by a factor of $\sim 3$. We see that for perturbations with $m=1$ it provides a somewhat larger difference in the sizes of humps in $G(t)$ on top left panels in Fig. 11 and 13. The reason for a spiral to shift out of the inner region of disc is that in P1 viscosity enhances as $r\to r_{i}$ which excludes the contribution of region with super-Keplerian shear rate close to $q=2$ to the transient growth of global perturbations.

Like in the other models considered before, we see a moderate agreement of $G(t)$ obtained for global perturbations $m=10$ with its analytical counterpart defined by equation (73) together with definitions (58), (59), the condition (60) and specifications of the model given in §5.1.1, see the top left panel in Fig. 13. Additionally, equation (60) yields $r_{sp}\approx 4.5$ and $r_{sp}\approx 6.5$, respectively, for dashed and dot-dashed curves on the bottom panel in Fig. 13, which is also in agreement with global calculations for $m=1$.

We finish the description of global optimal growth in the model P1 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 hump in curve of $G(t)$, see the top left panel in Fig. 11 and 13. For P2 $G_{max}$ is a few times smaller that local value, yet it raises 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$.

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, $\lambda_{y}$, with respect to $H$ and, particularly, vortices with $\lambda_{y}\gtrsim H$ which we refer to as ’large-scale vortices’ in contrast to small-scale vortices $\lambda_{y}\ll H$ received most of 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 parametrised 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, since 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 a special importance for large-scale vortices as their swing interval may be much longer than the rotational period of the flow. Secondly, vortices with $\lambda_{y}\sim 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 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 relaminarisation. In order to obtain optimal solutions we employ an iterative loop constructed of the basic and adjoint equations advanced, respectively, forward and backward in time. Note that the optimisation 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 (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 (73). It complements the known expression for $G_{max}$, see equation (71), which is valid for the small-scale vortices only, in the range of $\lambda_{y}\gg H$. Since equation (73) 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^{-1}$ rather than in units of an artificial $L^{-1}$, cf. Appendix A.1 and §2.1.

We start numerical analysis inspecting the influence of bulk viscosity on the growth of optimal SFH. The results are shown in Fig. 1, where one finds 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})$ being 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 vortical spiral is being contracted by the shear, i.e. it does not cause damping of 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. Armed with 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_{05}=2000$ the growth of shearing vortices in the boundless Keplerian shear attains a maximum value $g_{s}\sim 100$ at $k_{y}\sim 0.3$ expressed in units of $H^{-1}$, which is substantially larger in comparison with corresponding $G_{max}\sim 20$ obtained by MAN05 in an incompressible flow between the walls. The optimal growth declines slower as $k_{y}\to 0$, rather than as $k_{y}\to\infty$, which is in agreement with analytical estimates (71) and (73), see top panel in Fig. 1. In the case of a much larger Reynolds number typical for astrophysical discs turbulent energy reservoir can be, in principle, replenished by perturbations with $\lambda_{y}\gtrsim H$ provided that there is a subsequent delivery of their energy to SFH with negative $k_{x}<0$ on smaller azimuthal scales $\lambda_{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 et al. (2006) performed high resolution 3D numerical simulations in the large shearing box approximation with horizontal size $L>H$, however, it was not still enough to well resolve the region $k_{y}<1$, 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 enough elongated in azimuthal direction to survive for longer times. In light of our results, the high resolution simulations in box extended in 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 akin existing in density waves: its amplitude increases as $\tilde{k}_{x}\to 0$ and, vice versa, decreases as $\tilde{k}_{x}$ of SFH starts growing back. Also, as it is suggested by balanced solutions, the perturbation of enthalpy entering the large-scale vortex is proportional to the potential vorticity perturbation. By 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 reorganise a picture of non-linear interaction between SFHs, presumably, via the term $\nabla(\rho_{1}{\bf u})$ in the continuity equation. For example, the last term must be significant if it is evaluated for density perturbation generated by swinging leading spiral with $k_{y}\ll 1$ and velocity perturbation generated by swinging leading spiral with $k_{y}\gtrsim 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 (73), see Fig. 3. Thus, with the account of all perturbation scales as compared to $H$, it is confirmed that there is no such an abrupt depression of $G_{max}$, which takes place for perturbations $\lambda_{y}\ll 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 it is 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 raises faster than the amplitude of $\hat{W}$ as $q\to 2$. Additionally, as soon as $q\to 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 like $\propto k_{y}^{-4/3}$ exceeding the value of $G_{max}$ produced by streamwise streaks, cf. Fig. 2d by MAN05, at $k_{y}\lesssim 0.1$. At this point let us remark that Rincon et al. (2007) employed a sophisticated non-linear continuation methods 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 anti-cyclonic 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 scale-height, 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 $\alpha(q)$, see equation (28), where $\alpha$ 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, e.g., 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, since $k_{y}\sim\delta$ would correspond to SFH with $\lambda_{y}\sim r_{0}$, which is not local anymore. Therefore, a study of transient growth of large-scale shearing vortices would not be complete without an investigation of global azimuthal harmonics of perturbations $\propto\exp({\rm i}m\varphi)$ with $m\sim 1$. In §5 we additionally consider a global 2D dynamics of shearing vortices described by the set of equations (48-50) supplemented by global viscous terms (80-83). We provide a detailed construction of several global background models, introduce boundary conditions for both the basic and the adjoint variables and give brief notes on the numerical method. Note that here we leave out a problem of global density wave excitation setting $R_{b}\ll 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 scale height $H$, but additionally by radial profiles of (i) a rotational frequency, (ii) a background vorticity, (iii) a surface density and a viscosity, (iv) a shear rate itself, e.g., 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 particular disc model. Nevertheless, we manage to understand the behaviour of global transient spirals with $m\gg 1$ employing interpretation in terms of the local approach and 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\sim 1$ produce $G(t)$, which deviates from the analytical curve. However, they basically preserve the qualitative properties of spirals with $m\gg 1$. It turns out that global optimal vortex is always radially localised and, as one goes to longer timespans, 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 locations of initial optimal vortices $m=1$ obtained for different optimisation timespans in any of the disc models we used, see Fig. 7, 9, 11 and 13. In order to construct an analytical $G(t)$ we assume that (i) the optimisation time $t$ is equal to the time elapsed before the instant of swing $t_{s}$, (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 localisation of global spiral corresponding to this optimisation timespan. 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 estimation of $r_{sp}(t)$ and, subsequently, with 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 outskirts in comparison with similar spirals obtained in models with a constant viscosity and surface density. As soon as someone includes relativistic effects, cf. the models P1 and P2, this causes the suppression of $G_{max}$ down to 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 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 optimisation method applied at long timespans. Its profile is shown is Fig. 8 and Fig. 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 mode is located inside the flow and the increment attains a few percents 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 corotation point. Let us point out that this sort of modal growth competes with transient growth at Reynolds numbers equal to few thousands only. Since $G_{max}\propto R^{2/3}$, we expect that in astrophysical situation $R\sim 10^{10}\div 10^{13}$ 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 $\lambda_{y}\ll H$ or $\lambda_{y}\gtrsim 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 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 quasi-Keplerian regime. Thus, an interesting issue would be to find out whether the optimal666or, 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, $\lambda_{y}\gtrsim 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 is 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 generation of vortices of considerable amplitude, which are usually invoked to resolve the problems of angular momentum transport and planetesimal formation inside non-ionised parts of protoplanetary discs. Some of them 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 et al., 2013; Richard et al., 2016), whereas the 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 subcritical behaviour of stratified Keplerian flow it is important to note that a novel type of linearly growing vortical 3D perturbations appears due to the so called strato-rotational balance, see Tevzadze et al. (2003) and Tevzadze et al. (2008). The existence of these 3D vortices is closely related to special invariant of stratified fluid dynamics equal to scalar product of the entropy gradient and the vortensity. These vortices have been suggested as transient perturbations for bypass transition scenario in stratified discs, however, note that they degenerate into trivial solution within unstratified model considered in this paper. Also interestingly, that 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 may play in subcritical instability of stratified Keplerian flows, motivating to consider the corresponding 3D optimisation 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 optimisation problem in baratropic and vertically stratified backgrounds. Particularly, this would be necessary to do in the large-scale limit $\lambda\gg H$ in order to establish the link with the work done by Umurhan et al. (2006) and later Rebusco et al. (2009). Employing the asymptotic time-dependent model of geometrically thin viscous Keplerian disc, these authors obtained that the vertical sound waves exhibit non-modal growth being 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 et al. (1999) gets smaller as one considers the non-linear dynamics of shearing vortices with $k_{y}\lesssim 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 the usual phenomena in space and should be understood better in astrophysical science and, mainly, in theory of astrophysical discs.

We would like to finish this paper making link with a neighboring field of magnetohydrodynamic (MHD) turbulence in accretion discs. Actually, it became 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) (see also Nath & Mukhopadhyay (2015)) has shown that the established practice to oppose the subcritical HD transition to supercritical magnetohydrodynamic (MHD) transition to turbulence in astrophysical discs was misleading in a particular sense. It turned out that even in a spectrally unstable magnetised Keplerian shear the transient growth prevails the modal growth by orders of magnitude at certain short timespans. Thus, in fact, it is not yet clear what linear mechanism gives rise to a sustained MHD turbulence even in hot ionised accretion discs: either the magnetorotational instability, or the transient growth of magnetised 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 turbulisation of spectrally stable Keplerian background with zero net magnetic flux. The above extension of non-modal concepts onto 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 magnetised 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 remaining issues in both HD and MHD problems using a general non-modal approach.

Acknowledgements

We thank P. B. Ivanov for careful reading the manuscript and important remarks. DNR was supported by grant RSF 14-12-00146 for writing §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.

Appendix A Subcase of incompressible dynamics

A.1 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. On the contrary, here we are going to deal with an unbounded shear. As $a_{*}\to\infty$, the continuity equation (3) yields the following restriction on the velocity components

[TABLE]

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

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

[TABLE]

where the quantity $L$ is an auxiliary radial scale, which does not have any physical meaning. The coordinates (61) 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^{-1}$. Similarly to §2.1, we omit the primes and introduce the shearing radial wavenumber $\tilde{k}_{x}\equiv k_{x}+qk_{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:

[TABLE]

[TABLE]

provided that

[TABLE]

and it is assumed that $\hat{u}_{x},\hat{u}_{y}$ are expressed in units of $U\equiv\Omega_{0}L$. $\hat{W}$ is expressed in units of $U^{2}$. Correspondingly, $R\equiv UL/\nu$ is Reynolds number. Note that MAN05 used Reynolds number

[TABLE]

Using (62-64), we obtain equation for $\hat{u}_{x}$ which reads

[TABLE]

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

[TABLE]

The solution (67) is used to obtain the transient growth factor of incompressible perturbations for particular values of $k_{x},k_{y}$, $R$ and $q$, see §A.2 for further consideration.

A.2 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 simple expression for a single SFH:

[TABLE]

Substituting equation (67) into (68) yields the following growth factor of SFH:

[TABLE]

Given the particular values of $k_{y}$, $R$ and $q$, (69) has a global maximum over all $k_{x}$ and $t$ which we are going to denote here as

[TABLE]

since this quantity is analogous to the maximum optimal growth defined in §3, see equation (21).

Appendix B Analytical estimations

of optimal growth

By means of the analytical supplement to their numerical optimisation 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 its Fig. 3b). In Fig. 15 we take the Keplerian shear, $q=3/2$, $R_{05}=2000$ and plot equation (77) of MAN05 (dot-dashed line) together with our quantity (70) obtained numerically (solid line). It can also be checked that equation (77) of MAN05 taken for a boundless shear,

[TABLE]

virtually coincides with our numerical curve in Fig. 15. Equation (71) can be easily obtained from (69) setting $\tilde{k}_{x}=0$ (the swing instant of SFH) and taking the limit $|k_{x}|\gg k_{y}$. Thus, in general, $G_{max}\to\infty$ as $k_{y}\to 0$ in the absence of boundaries. Physically, this takes place due to a simple fact that for a tightly wound SFH with $|k_{x}|\gg k_{y}$ both viscous dissipation time-scale, $t_{\nu}$, 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 (69). As it is discussed, e.g., by RZ15, equation (71) can be approximately derived suggesting that $G_{max}$ corresponds to equality $t_{\nu}=t_{tg}$, which gives the time when the maximum growth occurs as

[TABLE]

Note that estimation (72) provides exactly the same value which one gets looking for a maximum of (69) in order to derive an approximate $G_{max}$ given by equation (71). Then, using the inviscid form of $g\approx(qt)^{2}$, we recover equation (71) up to the damping factor $\exp(-2/3)$.

Of course, the growth factor (69) and so $G_{max}$ do not depend on the scaling factor $L$. As we see from equation (71), the transient growth value is controlled by the ratio $R/k_{y}^{2}\sim(\lambda_{y}/l_{p})^{2}$, where $\lambda_{y}$ is characteristic azimuthal length of perturbations and $l_{p}$ is mean free path of gas particles which determine microscopic viscosity. In fact, $\lambda_{y}^{-1}\sim k_{y}L^{-1}$ of SFH exhibiting $G_{max}$ cannot be too low since $|k_{x}|\sim(qk_{y}R)^{1/3}$ also decreases: but $|k_{x}|L^{-1}$ cannot be less than $H^{-1}$, since in the last case we start dealing with patch of disc comparable to its thickness in the direction of shear what forces us to take into account the effects of compressibility, see e.g. the discussion about SFH swing interval in §3.2.1 of ZR14. Thus, we of course cannot reach an arbitrary high $G_{max}$ at the expense of $\lambda_{y}$ value. It turns out, see §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 $\lambda_{y}\gg H$ (i.e. $k_{y}\ll 1$) employing the balanced solutions to equations (7-9), see Heinemann & Papaloizou (2009). The inviscid result reads $g_{opt}\approx 4\kappa^{-4}(k_{y}qt)^{2}$, where $\kappa^{2}=2(2-q)$ is the dimensionless epicyclic frequency squared expressed in units of $\Omega_{0}^{2}$. Substituting $t_{max}$ given by equation (72) into this result we obtain

[TABLE]

where damping factor $\exp(-2/3)$ has been added since we expect a similar action of viscous forces on shearing vortices having both $k_{y}\gg 1$ and $k_{y}\ll 1$, see a detailed justification of this assumption in §4.1. As we see, one gets 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}\to 0$, see Fig. 15. At the same time, as it is seen in Fig. 15, one expects much larger transient growth of columnar perturbations (i.e. perturbations independent of vertical direction) in the domain of $\lambda_{y}\gtrsim 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$777or, equivalently, outside of 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 e.g. 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}$:

[TABLE]

[TABLE]

where $K\equiv k^{2}+\kappa^{2}$ 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 corresponding expressions e.g. of Bodo et al. (2005) or Tevzadze et al. (2003) for the case of baratropic perturbations. $I_{\nu}$ does not depend on time in the inviscid limit, however, in presence of non-zero viscosity this is not true anymore. To emphasise this fact, we add a subscript ’$\nu$’ to the notation of the potential vorticity perturbation. It is not difficult to construct equation for $I_{\nu}$ using the set of equations (7-9):

[TABLE]

where the second term in the right-hand side of equation (76), 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\gg 1$, $R_{b}\gg 1$. Thus, in the leading orders in $R^{-1}$ and $R_{b}^{-1}$ one can use an inviscid balanced solution for $\hat{W}$ (see equation (37) of RZ15) which gives that $|\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 in the right-hand side of equation (76). Similarly, in the leading orders in $R^{-1}$ and $R_{b}^{-1}$, it is sufficient to take the slow varying solution of inhomogeneous equation (75) in its inviscid form,

[TABLE]

but with $I_{\nu}$ decaying as described by equation (76). This is the decay of $I_{\nu}$ which yields an additional factor $\exp(-2/3)$ in equation (73) being evaluated at the $t_{s}=-k_{x}/(qk_{y})$ with $t_{s}$ additionally equated to (72). Correspondingly, $\hat{u}_{x}$ is given by

[TABLE]

see equation (35) by 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:

[TABLE]

The tightly wound vortical spirals are of primary interest here, since we are dealing with high $R$ and focus on $k_{y}\lesssim 1$ which leads to $t_{max}\gg 1$, see equation (72). Further, since $t_{max}$ approximately equals to $t_{s}$ of an optimal SFH exhibiting $G_{max}$, at the initial moment $|k_{x}|\gg k_{y}$ and the same is about $|\tilde{k}_{x}|$ until $t>0$ is not too close to the instant of swing.

Equation (79) implies that terms in the second square brackets in the right-hand sides of equations (74) and (75) 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_{b}^{-1}$ are absent in equation (76), we come to the conclusion that bulk viscosity cannot affect the dynamics of vortical SFH even if $R_{b}\sim 1$.

The opposite situation takes place with density waves. These are rapidly oscillating solutions to the homogeneous equations (27-29) of RZ15. That 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_{b}^{-1}\tilde{k}_{x}^{2}\dot{\hat{u}}_{x}$ becomes the leading one among all viscous terms in equation (74) and it exceeds the inviscid terms as soon as $R_{b}\lesssim|\tilde{k}_{x}|^{2}$.

Appendix D Supplementary material to the global problem

D.1 Viscous terms

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

[TABLE]

[TABLE]

[TABLE]

[TABLE]

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

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

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 the Appendix D.2. Also, similar to what we have in the local problem, the terms in equations (84), (85), (87) and (88) acquire opposite signs comparing to corresponding terms in equations (80-83). Additionally, the term reciprocal to term $\propto\delta h$ in equation (81) appears in the adjoint equation for $\delta\tilde{h}$, see equation (55). Note that the last term vanishes as one turns to local perturbations.

D.2 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 timespan.

D.2.1 Models N1 and P1

At first, let us consider the set of equations for projections of the displacement vector $\mathbf{\xi}$. It can be done using a general relationship between $\mathbf{\xi}$ and the velocity Lagrangian perturbation (Lynden-Bell & Ostriker, 1967, Eq. 16):

[TABLE]

[TABLE]

Since $v_{r}$ tends to infinity at $r_{i}$ (see equation (37)), $\xi_{r}$ and $\xi_{\varphi}$ must vanish at $r_{i}$. This implies that the boundary remains unperturbed, which simplifies the boundary conditions. Now, as $v_{r}\propto D^{-3}$ in the vicinity of the boundary, the components of the displacement vector must behave at least as

[TABLE]

[TABLE]

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, $\sigma_{ik}$, can be found, e.g., in Landau & Lifshitz (1987, §15). Since $\xi_{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 $\sigma_{ik}$:

[TABLE]

[TABLE]

Integration of equation (94) over $z$ yields

[TABLE]

In equation (LABEL:boundary1) all variables are taken at the inner boundary of disc. However, because of the relation (38) and limitations on the displacement vector (91), (92) equation (LABEL:boundary1) is satisfied for arbitrary $\delta v_{r}$ and $\delta v_{\varphi}$. Thus, equation (LABEL:boundary1) does not provide us with any restriction on perturbations at $r_{i}$.

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

[TABLE]

In equation (LABEL:boundary2) all variables are taken at the inner boundary of disc. Equation (LABEL:boundary2) allows us to formulate the necessary restrictions on $\delta v_{r}$ and $\delta v_{\varphi}$ at $r_{i}$, provided that it is supplied by the regularity condition for $\delta h=\delta p/\rho$ at the boundary. Indeed, equation (LABEL:boundary2) yields

[TABLE]

The last two terms vanish due to limitations (91), (92). Since $\nu$ diverges as $r\to r_{i}$, $\delta h$ is regular at $r_{i}$ as soon as

[TABLE]

[TABLE]

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

[TABLE]

[TABLE]

[TABLE]

D.2.2 Models N2 and P2

In this case we impose the no-slip boundary conditions at $r=r_{i}$:

[TABLE]

[TABLE]

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

[TABLE]

[TABLE]

at $r=r_{i}$. Note that for models N2 and P2 the boundary condition onto $\delta\tilde{h}$ is not required, since the advective terms are absent in this case, see explanation in §5.4.