© ESO 2018

1. Introduction

Rotating boundary layers are believed to exist in the vicinity of weakly magnetised stars accumulating the material from an accretion disc. A slowly spinning star surface causes transformation of the rotational energy of the disc into thermal energy inside the boundary layer. Thus, the boundary layer is expected to have a brightness comparable to that of the accretion disc and be the main source of radiation in a relatively higher frequency band than the disc itself. The structure of the boundary layer and the adjacent disc inner part must be determined together (see Bisnovatyi-Kogan 1994). A number of boundary layer models have been developed for various classes of objects. Usually they employ slim disc equations (Abramowicz et al. 1988) in order to match consistently the disc inner part to the boundary layer. Among them the accreting pre-main sequence stars were studied by Popham et al. (1993), whereas the accreting white dwarfs and neutron stars were studied by Narayan & Popham (1993) and Popham & Sunyaev (2001). In the latter case, the formation of a spreading layer is also possible in systems with high enough accretion rate, as was suggested by Inogamov & Sunyaev (1999).

Just as for the formation of an accretion disc, the cornerstone of boundary layer formation is the origin of effective shear viscosity providing enhanced angular momentum and mass transfer onto the star surface. Conventionally, it is parametrised by α, which is the kinematic viscosity coefficient scaled by the speed of sound and pressure radial scaleheight (see Shakura & Sunyaev 1988). This variant of viscosity prescription is similar to what is done in the disc models, where α is scaled by the disc thickness (see Shakura & Sunyaev 1973). However, the physical process veiled by an effective viscosity as well as magnitude of α remains a matter of debate, probably, to a greater extent, in the case of boundary layer rather than in the case of accretion disc. Turbulence is known to be a natural solution to this issue, but its simple (supercritical) variant is discarded by the centrifugal stability of the rotating shear flows, which represent both accretion discs and boundary layers. A “magic wand” working in hot magnetised discs, where angular velocity Ω decreases with distance to the rotation axis r, is a magnetorotational instability. However, it does not work if dΩ²/dr > 0, which is the case for the boundary layers. Such flows are locally linearly stable even in the presence of the magnetic field. In this context, an alternative mechanism of angular momentum transfer has been proposed by Belyaev et al. (2013) and more recently by Philippov et al. (2016). It was shown that global sonic instability of the flow associated with the presence of the star surface excites global non-axisymmetric acoustic modes, which are responsible for effective viscosity even in the absence of turbulence.

At the same time, the possibility that astrophysical boundary layers acquire the effective viscosity through the turbulence is not ruled out, since flows with shear rates including that of the order of |Ω| may be locally non-linearly unstable. Along with the hydrodynamical, non-linear stability of homogeneous Keplerian shear flow, which locally represents the accretion and protoplanetary discs, this is one of the major unresolved problems in astrophysical fluid dynamics. Below we consider both problems jointly. In order to be non-linearly unstable, such flows first require a candidate for transiently growing perturbations, which could be amplified by background shear by orders of magnitude for the Reynolds number R typical in astrophysical situations; second, in phase of high amplitude, those perturbations have to give a positive non-linear feedback to the basin of small perturbations subject to transient growth. Two of these conditions together introduce what is usually called a bypass scenario for transition to turbulence.

The bypass scenario is best developed in application to plane shear flows and, in particular, to the Couette flow (see Hamilton et al. (1995) and Waleffe (1997), who proposed a non-linear mechanism sustaining the turbulence). This mechanism is called the self-sustaining process (SSP) and is sketched, for example in Fig. 1 by Waleffe (1997). The transient growth of initially small streamwise rolls produces high amplitude streamwise streaks, the process usually called lift-up effect after the work by Ellingsen & Palm (1975). As streamwise streaks modify the background velocity yielding an inflexion point in a spanwise direction, the flow becomes linearly unstable with respect to non-axisymmetric modes, which in turn regenerate streamwise rolls via non-linear interactions with each other.

However, the situation becomes less clear in centrifugally stable rotating shear flows. Let us define rotation number, which characterises a dynamical contribution of rotation to shear, R_Ω = −2/q, where q = −(r/Ω)dΩ/dr is the local dimensionless shear rate. Steady, non-linear self-sustaining solutions obtained in the non-rotating case, for which R_Ω = 0, can be non-linearly continued into both the centrifugally unstable regime, −1 < R_Ω < 0 (q > 2), and the cyclonic regime with high shear rate, R_Ω ≪ 1 (−q ≫ 1) (see Rincon et al. 2007). Nevertheless, this operation cannot be extended either until the centrifugally stable anti-cyclonic (or quasi-Keplerian) regime with R_Ω ≤ −1 (0 < q < 2), or far into the cyclonic regime with R_Ω ≲ 1 (q ≲ −1). Presumably, this takes place because the lift-up mechanism ceases to work in both cases. On the other hand, direct numerical simulations (DNS) allow us to obtain turbulent solutions until, at least, R_Ω ≈ 0.3 (see Lesur & Longaretti 2005; LL05 hereafter). It is interesting to note that such a value of rotation number is far beyond the interval of small positive R_Ω corresponding to a substantial transient growth of streamwise rolls (see Sect. 5.1.1 of this paper). Another notable fact is that the quasi-Keplerian regime demonstrates a strong asymmetry in comparison with the cyclonic regime. Namely, DNS show dramatic stabilisation of the quasi-Keplerian shear, as one goes away from the Rayleigh line, R_Ω = −1, to smaller rotation numbers (cf. Figs. 4 and 7 of LL05; see also the earlier DNS of the quasi-Keplerian regime made by Hawley et al. 1999). We have already noted that when R_Ω = −1 exactly, the lift-up mechanism does not work, however, there is an inverse process of small streamwise streaks growing into high amplitude streamwise rolls, which is called the anti-lift-up effect (see Sect. 2.3 in Rincon et al. 2008 for more detail). The regime R_Ω = −1 is neutrally stable with respect to local modal perturbations, that is, the growth of streaks can be regarded as epicyclic oscillations with an infinite period. Because finite amplitude streamwise rolls of the anti-lift-up are quite different from finite amplitude streamwise streaks of the lift-up, Rincon et al. (2007) put forward several arguments that the turbulent solution should be represented by a bypass scenario other than SSP. In any case, since a bypass transition in the regime R_Ω = −1 should involve the anti-lift-up mechanism, it looks natural that LL05 observed a sharp increase of the transition Reynolds number R_T as one goes to R_Ω < −1, in contrast to the situation with the cyclonic regime. Indeed, the streak’s growth is strongly suppressed at R_Ω < −1, as it turns into epicyclic oscillations with nearly rotational frequency (see also Sect. 5.1.2 of this paper). LL05 observed a decay of initial turbulence in the ranges −∞ < R_Ω < −1.03 and 0.3 < R_Ω < ∞ for R ≲ 10⁵ within a local model of the shearing box with shearing boundary conditions. A bypass scenario proposed in this work naturally resolves the issues mentioned above.

There are a number of studies that examine properties of centrifugally stable flows either experimentally or employing DNS in the framework of the Taylor-Couette problem (the flow between the coaxial corotating cylinders, see Sect. 6 of the review by Grossmann et al. 2016). Among the recent ones are Burin & Czarnocki (2012), who observed the turbulisation of the cyclonic flow at R_T ∼ 10⁵ in the case of R_Ω ≈ 0.8, and Ostilla-Mónico et al. (2016), who considered the non-stationary dynamics in the same regime numerically excluding the influence of axial boundaries, which is unavoidable in the laboratory setup. At the same time, the quasi-Keplerian regime was investigated in the experiment by Schartman et al. (2012) subsequently continued by Edlund & Ji (2014), and numerically by Ostilla-Mónico et al. (2014) subsequently continued by Shi et al. (2017). All those studies demonstrated an outstanding stability of quasi-Keplerian flow with respect to finite amplitude perturbations from R ∼ 10⁵ up to higher than R ∼ 10⁶. Yet, this result should be accepted with a caution, because the presence of walls may prevent the onset of self-sustaining solutions inherent to the anti-cyclonic subcritical regime as was argued by Rincon et al. (2007). Additionally, let us emphasise that most of attempts to detect turbulence have been made at super-Keplerian shear rates, 3/2 < q < 2. For example, Edlund & Ji (2014) and Shi et al. (2017) used q = 1.8 and q = 5/3, respectively¹. Indeed, there is a tendency to consider the shear rates starting from q = 2, which stems from the established opinion that the closer quasi-Keplerian shear is to the Rayleigh line, R_Ω = −1 (q = 2), and the farther it is from the solid-body line, R_Ω → −∞ (q = 0), the less stable it becomes. However, an accompanying purpose of this work is to show that the most appropriate shear rate to test the transition to turbulence in the quasi-Keplerian regime may be found in the vicinity of q = 1/2 (see Sect. 7).

Two decades ago it was established that shear flow turbulence is a result of the subtle interplay between non-normal linear dynamics of perturbations and their non-linear interaction (see for example the review by Grossmann 2000). Although the subcritical transition is in essence a loss of non-linear stability, the phase of linear, non-modal transient growth is of much importance in the general design of shear flow turbulence. Among others, this was elucidated by Henningson (1996). It was argued that (i) transient growth as linear process is the only source of energy for subcritical turbulence, (ii) transient growth must be large, that is the linear growth factor must be much greater than unity, for the transition to occur: for particular examples, the plane Couette and Poiseille flows may become turbulent at the Reynolds numbers as low as, respectively, R_T ≃ 350 and 1000, when according to Butler & Farrell (1992) and Reddy & Henningson (1993) the maximum growth factor of streamwise rolls is, respectively, ≃150 and 200. Also, a number of low-dimensional models of plane shear flow employing various representations of the non-linear dynamics demonstrated qualitatively similar behaviour during the transition: substantial non-modal growth followed by the “bootstrapping” effect (Trefethen et al. 1993; see Baggett & Trefethen 1997). These authors explained threshold exponents for transition to turbulence in each case using simple arguments based on the secondary linear instability of the transient perturbation attained the high amplitude. Thus, it was shown that transient dynamics can always be distinguished in the whole loop of the sustenance of turbulence (see caption for Fig. 3 by Baggett & Trefethen 1997)).

The reasoning just mentioned makes the transient growth a pivotal process responsible for the generation of the specific high-amplitude inputs for the sequential non-linear feedback. Such high-amplitude inputs are the streamwise streaks in the case of plane shear flows, but what are they in the case of centrifugally stable homogeneous and boundless rotating shear flow? To answer this question in this paper, we (i) consider three-dimensional (3D) perturbations, which are able to exhibit large transient growth, (ii) perform a critical test, which indicates that high-amplitude perturbations generated during this transient growth have to do with the sustenance of turbulence. Issues (i) and (ii) are tackled solving, respectively, the linear problem (1)–(4) and the more general non-linear problem (53)–(55).

It is known that shear flow allows for the existence of the so-called Kelvin modes, or alternatively, shearing vortices exhibiting the transient growth via the swing amplification² (see Chagelishvili et al. 2003). The shearing vortices are spanwise-invariant, nearly streamwise perturbations contracted by a background motion. In rotating flow, they take the form of spirals, whereas the instant of swing corresponds to the transformation of the leading spirals into trailing spirals or vice versa (see e.g. Razdoburdin & Zhuravlev 2015). Conservation of a spanwise vorticity perturbation is the physical reason for the growth of their velocity and pressure perturbations. The swing amplification is a two-dimensional (2D) process, which takes place independently of whether the shear is plane or rotating: the Coriolis force does not affect the growth of shearing vortices in contrast to previously discussed rolls and streaks. However, shearing vortices grow considerably weaker rather than rolls and streaks, since they are much more suppressed by viscosity. The swing amplification gives the kinetic energy growth by a factor ∝ R^2/3, rather than ∝ R² as in the case of lift- and anti-lift-up. Mukhopadhyay et al. (2006) estimated that R ∼ 10⁶ is required to obtain the growth factor of shearing spirals comparable to the growth factor of rolls necessary to activate the SSP in the plane Couette flow. Remarkably, this claim is valid for both the plane and the rotating (no matter which cyclonicity) flows and does not change with the shear rate, since it can be shown that the maximum growth factor of shearing vortices G_max ≈ R^2/3exp ( − 2/3) provided that the time is measured in the inversed shear rate rather than in the inversed rotational frequency (cf. Eq. (83) by Afshordi et al. (2005) and Eq. (9) by Mukhopadhyay et al. 2005). Thus, the role of shearing vortices is negligible either at the Rayleigh line, or in absence of rotation, but as soon as R_Ω is either sufficiently smaller than −1, or sufficiently larger than 0, the growth of streaks and rolls vanishes and the shearing vortices become the kind of perturbations exhibiting the largest transient growth in the flow. This picture was confirmed for a variety of models by Ioannou & Kakouris (2001), Mukhopadhyay et al. (2005), Yecko (2004), and more recently by Maretzke et al. (2014), Zhuravlev & Razdoburdin (2014), and Razdoburdin & Zhuravlev (2017) who carried out studies of the transient dynamics of linear perturbations in the quasi-Keplerian regime using rigorous optimisation methods, which allow us to obtain the optimal perturbations exhibiting the largest possible transient growth.

For very high Reynolds numbers typical in astrophysical boundary layers and the adjacent accretion discs, the growth of shearing spirals becomes quite significant. However, its role in the transition to turbulence remains unclear. On the one hand, it seems plausible that the shearing spirals are involved in 2D models of subcritical turbulence in shear flows simulated by Umurhan & Regev (2004), Johnson & Gammie (2005), and Horton et al. (2010). Furthermore, Lithwick (2007) suggested a 2D variant of positive non-linear feedback sustaining the growth of shearing spirals. He showed that the coupling between shearing spirals at the instant of swing and small-amplitude axisymmetric vortices generates a new shearing spiral capable of transient growth. On the other hand, Lithwick (2009) revealed that the streamwise scale of a shearing spiral must exceed the disc scale height to make this process work in the 3D model, otherwise the shearing spiral is destroyed by resonant interaction with a 3D axisymmetric vortex (see also the simulations by Shen et al. 2006). Consequently, it seems reasonable that shearing spirals are only relevant to 2D turbulence, which may occur at scales above the accretion disc scale height: moreover, Razdoburdin & Zhuravlev (2017) have recently confirmed that shearing spirals keep the ability for significant transient growth even when their streamwise scale is comparable to the disc global radial scale. Also, the lack of the positive 3D feedback from shearing spirals at phase of swing is indirectly confirmed by the DNS of LL05, which show that R_T substantially depends on the shear rate for R_Ω ≲ 1 in the cyclonic case. In turn, this means that G_max strongly varies along the curve of R_T(q). Hence, positive non-linear feedback would come into play at a different amplitude of working perturbation $\propto \sqrt{G_{\max }}$ for different shear rates, which would require further clarification³. Additionally, Maretzke et al. (2014) revealed a small difference between the optimal perturbations in cyclonic and quasi-Keplerian regimes of the flow between the cylinders: in the cyclonic regime, the optimals are not shearing spirals but their slightly modified counterparts with a weak spanwise dependence. Coincidently or not, the non-linear instability of the cyclonic flow between cylinders is already observed, rather than that of the quasi-Keplerian flow. To summarise, we suppose that the shearing spirals considered as spanwise invariant perturbations seem unlikely to be an ingredient of the non-linear feedback in the bypass scenario of subcritical transition to 3D turbulence on scales small compared to the disc scale height.

In this work we consider analytically the linear local dynamics of quasi-2D shearing spirals in the rotating, homogeneous, boundless, and centrifugally stable shear flow. By a quasi-2D shearing spiral, we mean a spatial Fourier harmonics (SFH) of small perturbation, which initially has a zero spanwise velocity perturbation, but a small non-zero spanwise wavenumber. Such perturbations are not optimal, meaning they do not attain the largest possible growth factor measured at the instant of swing. Instead, during the swing amplification they generate essentially 3D vortices of a new type, which we refer to as cross-rolls, since those are aligned in the shearwise direction in contrast to ordinary streamwise rolls generated via the anti-lift-up mechanism on the Rayleigh line (see Fig. 9). We show that while the “planar” eddies constructed of shearwise and streamwise velocity perturbations decay as the shearing spiral is wound back by a background flow, the cross-rolls conserve their amplitude and fade due to viscous dissipation only. The growth factor of cross-rolls with respect to the growth factor of the shearing spiral at the instant of swing is less than unity and decreases when the shear rate alters from q = −∞ up to q = 2. Thus, the cross-rolls amplitude is larger in the cyclonic regime, but smaller in the quasi-Keplerian regime. Also, it vanishes both at the solid-body line and at the Rayleigh line. We note that the growth factor of cross-rolls can be arbitrarily high provided that the initial quasi-2D shearing spiral is sufficiently tightly wound. We suggest that the non-linear instability and corresponding self-sustained turbulent solutions in the regime sufficiently far from both the centrifugally unstable interval and the non-rotating case are provided by the positive non-linear feedback associated with either a secondary linear instability of the flow containing the finite amplitude cross-rolls, or the interaction of cross-rolls with each other. Either of these two processes may generate new small quasi-2D shearing spirals capable of substantial swing amplification. Evidence for that comes from our simulations of subcritical transition to turbulence in the cyclonic regime. It turns out that in the shearing box model with periodic boundary conditions the transition is facilitated in the box stretched along the rotation axis, which indicates that perturbations with small but non-zero spanwise wavenumber are involved in that process. Indeed, the numerical domain of finite size does not affect the 2D shearing spirals while it works as an ideal high-pass filter for harmonics with non-zero spanwise wavenumbers: it completely cuts off wavelengths larger than spanwise size of the numerical domain, whereas the first wavelength smaller than this size and satisfying the periodic boundary conditions is the best resolved by the numerical scheme. The inherent “flaw” of shearing box simulations becomes the useful tool in our case: the cross-rolls with streamwise size comparable to the box streamwise size, so the least susceptible to viscous damping, are forbidden in a cubic box but allowed in a tall box. This way, in the tall box we manage to find the transition up to q = −3 (R_Ω ≈ 0.67). Furthermore, we find that R_T obtained in the tall box simulations corresponds approximately to the constant maximum growth factor of cross-rolls in the range q ∼ −35 ÷ −3. Supposing that the transition to self-sustained turbulence in the quasi-Keplerian regime corresponds to the same threshold growth factor of cross-rolls, we predict the transition Reynolds number for 0 < q < 2 (see Fig. 8), which is R_T ≲ 10⁸ for the Keplerian shear rate. This is consistent with the results of LL05 concerning the asymmetry of R_T(q) in the vicinity of R_Ω = 0 and R_Ω = −1. Besides, it seems unrealistic to obtain the transition in the quasi-Keplerian regime considering the super-Keplerian shear rates. If our suggestion is true, the strategy of searching for hydrodynamical turbulence in the Keplerian flow should be changed: one should extend the turbulent solution from the cyclonic regime with high negative shear rate to the quasi-Keplerian regime across the solid-body line, R_Ω = +∞ →R_Ω = −∞, rather than start from a turbulent solution obtained on the Rayleigh line, which has a different origin and properties. We predict that the most favourable value of the shear rate to test the existence of hydrodynamical turbulence in the quasi-Keplerian regime is q = 1/2.

In the case of plane shear flow, R_Ω = 0, SFH with arbitrary non-zero spanwise wavenumber is described by an exact analytical solution obtained by Chagelishvili et al. (2016). At the limit of zero spanwise wavenumber, quasi-2D shearing spirals turn into SFH, which are usually called shearing spirals. From now on, we omit prefix the “quasi-2D”, so that by shearing spirals we mean all SFH with initially planar velocity perturbation, which may have a small non-zero spanwise wavenumber.

2. Linear perturbations in rotating shear flow

Local vortical 3D perturbations in viscous homogeneous rotating shear flow obey the following equations:

$$\begin{array}{c}\hfill (\frac{\partial }{\partial t}-q{\mathrm{\Omega }}_{0}x\frac{\partial }{\partial y})u_x-2{\mathrm{\Omega }}_0{u}_y=-\frac{1}{{\rho }_0}\frac{\partial p}{\partial x}+\nu {\mathrm{\nabla }}^2{u}_x,\end{array}$$(1)

$$\begin{array}{c}\hfill (\frac{\partial }{\partial t}-q{\mathrm{\Omega }}_{0}x\frac{\partial }{\partial y})u_y+(2-q){\mathrm{\Omega }}_0{u}_x=-\frac{1}{{\rho }_0}\frac{\partial p}{\partial y}+\nu {\mathrm{\nabla }}^2{u}_y,\end{array}$$(2)

$$\begin{array}{c}\hfill (\frac{\partial }{\partial t}-q{\mathrm{\Omega }}_{0}x\frac{\partial }{\partial y})u_z=-\frac{1}{{\rho }_0}\frac{\partial p}{\partial z}+\nu {\mathrm{\nabla }}^2{u}_z,\end{array}$$(3)

$$\begin{array}{c}\hfill \frac{\partial u_x}{\partial x}+\frac{\partial u_y}{\partial y}+\frac{\partial u_z}{\partial z}=0,\end{array}$$(4)

where u_x, u_y, and u_z are the Eulerian perturbations of velocity components, ρ₀ is a constant background density, and p is the Eulerian perturbation of pressure. By definition,

$$\begin{array}{c}\hfill {\mathrm{\nabla }}^2\equiv \frac{{\partial }^2}{{\partial }^{2}x}+\frac{{\partial }^2}{{\partial }^{2}y}+\frac{{\partial }^2}{{\partial }^{2}z}\end{array}$$

and kinematic viscosity ν is assumed to be a constant. Variables x, y, and z are Cartesian coordinates, which locally correspond to radial, azimuthal, and axial directions in boundary layer (or disc), respectively, whereas Ω₀ is angular velocity of fluid rotation at x = 0 corresponding to some radial distance r₀ ≫ x. Also, q is a constant shear rate that defines the background azimuthal velocity as v_y = −qΩ₀x. Equations (1)–(4) are vertically unstratified linearised versions of small shearing box equations (see A.3 Umurhan & Regev 2004). Throughout the text, the variables x, y, z are alternatively referred to as, respectively, shearwise, streamwise, and spanwise coordinates (or directions). The shear rate is positive (negative) in the case of anti-cyclonic (cyclonic) flow. The anti-cyclonic case is usually divided into Rayleigh-stable regime 0 < q ≤ 2, also referred to as the quasi-Keplerian regime with particular q = 3/2 corresponding to the Keplerian flow, and Rayleigh-unstable regime q > 2 subject to centrifugal linear instability. The solid-body line and the Rayleigh line are represented, respectively, by q = 0 and q = 2, whereas q → ±∞ are the limits corresponding to plane shear flow. In what follows, we are interested in the whole range −∞ < q < 2.

2.1. Dimensionless equations for shearing harmonics

We consider a boundless shear, which allows us to look for the partial solutions of Eqs. (1)–(4) in the form of SFH,

$$\begin{array}{c}\hfill f=\widehat{f}(k_x,k_y,k_z,t^{\prime })exp(\mathrm{i}{k}_x{x}^{\prime }+\mathrm{i}{k}_y{y}^{\prime }+\mathrm{i}{k}_z{z}^{\prime }),\end{array}$$(5)

where f is any of perturbation quantities and f̂ is its Fourier amplitude specified by the dimensionless wavenumbers (k_x, k_y, k_z). By Eq. (5) it is implied that the problem is considered with respect to the dimensionless comoving Cartesian coordinates

$$\begin{array}{c}\hfill x^{\prime }=x/L,y^{\prime }=(y+q{\mathrm{\Omega }}_{0}xt)/L,z^{\prime }=z/L,t^{\prime }={\mathrm{\Omega }}_{0}t,\end{array}$$(6)

where L is an auxiliary distance. Correspondingly, k_x, k_y, and k_z are expressed in units of L⁻¹. Let us assume from now on that k_y > 0 and L is chosen in such a way that k_y = 1. We omit the primes everywhere below and note that with respect to the original coordinates SFH acquires a varying radial wavenumber k̃_x ≡ k_x + qt. We arrive at the following dimensionless equations for SFH:

$$\begin{array}{c}\hfill \frac{\mathrm{d}{\widehat{u}}_x}{\mathrm{d}t}=2{\widehat{u}}_y-\mathrm{i}{\tilde{k}}_x\widehat{W}-R^{-1}(k^2+k_z^2){\widehat{u}}_x,\end{array}$$(7)

$$\begin{array}{c}\hfill \frac{\mathrm{d}{\widehat{u}}_y}{\mathrm{d}t}=-(2-q){\widehat{u}}_x-\mathrm{i}\widehat{W}-R^{-1}(k^2+k_z^2){\widehat{u}}_y,\end{array}$$(8)

$$\begin{array}{c}\hfill \frac{\mathrm{d}{\widehat{u}}_z}{\mathrm{d}t}=-\mathrm{i}{k}_z\widehat{W}-R^{-1}(k^2+k_z^2){\widehat{u}}_z,\end{array}$$(9)

$$\begin{array}{c}\hfill {\tilde{k}}_x{\widehat{u}}_x+{\widehat{u}}_y+k_z{\widehat{u}}_z=0,\end{array}$$(10)

where velocities are expressed in units of Ω₀L and pressure normalised by a constant background density, W ≡ p/ρ₀, is expressed in units of L²Ω₀². The non-zero kinematic viscosity yields a finite value of the Reynolds number

$$\begin{array}{c}\hfill R\equiv \frac{{\mathrm{\Omega }}_0{L}^2}{\nu }\end{array}$$(11)

and $k^2\equiv 1+{\tilde{k}}_x^2.$

The normalisation used above is somewhat different from the accepted one (see e.g. LL05 and Mukhopadhyay et al. 2005). Specifically, the time is measured in units of Ω₀⁻¹ rather than in units of (qΩ₀)⁻¹. We choose the former variant, since it allows us to cross continuously the solid-body line.

Next, as soon as we are dealing with vortical dynamics, it is convenient to operate with equations for û_x and a shearwise component of vorticity SFH defined as ω̂ ≡ k × û:

$$\begin{array}{c}\hfill {\widehat{\omega }}_x={\widehat{u}}_z-k_z{\widehat{u}}_y.\end{array}$$(12)

Equations (7)–(10) yield

$$\begin{array}{c}\hfill \frac{\mathrm{d}{\widehat{u}}_x}{\mathrm{d}t}=-2\frac{k_z}{k^2+k_z^2}{\widehat{\omega }}_x-2q\frac{{\tilde{k}}_x}{k^2+k_z^2}{\widehat{u}}_x-R^{-1}(k^2+k_z^2){\widehat{u}}_x,\end{array}$$(13)

$$\begin{array}{c}\hfill \frac{\mathrm{d}{\widehat{\omega }}_x}{\mathrm{d}t}=(2-q)k_z{\widehat{u}}_x-R^{-1}(k^2+k_z^2){\widehat{\omega }}_x.\end{array}$$(14)

3. Production of 3D vortices by shearing spirals: inviscid dynamics

We assume throughout this section that R → ∞.

3.1. Swing amplification

Here we revisit the well-known mechanism of swing amplification of vortices, which takes place in the xy-plane. Namely, in the case û_z = 0, k_z = 0, the problem (13)–(14) becomes one-dimensional and the evolution of the corresponding 2D SFH is described by a single equation

$$\begin{array}{c}\hfill \frac{\mathrm{d}{\widehat{u}}_x}{\mathrm{d}t}=-2q\frac{{\tilde{k}}_x}{k^2}{\widehat{u}}_x,\end{array}$$(15)

provided that the velocity is divergence-free,

$$\begin{array}{c}\hfill {\widehat{u}}_y=-{\tilde{k}}_x{\widehat{u}}_x.\end{array}$$(16)

The solution of Eq. (15) is

$$\begin{array}{c}\hfill {\widehat{u}}_x=\frac{k_0}{k^2},\end{array}$$(17)

where the normalisation k₀ ≡ k(t = 0) is obtained from the condition of doubled kinetic energy density of 2D SFH equals to unity:

$$\begin{array}{c}\hfill ({\widehat{u}}_x^2+{\widehat{u}}_y^2){|}_{t=0}=1.\end{array}$$

Accordingly, we have

$$\begin{array}{c}\hfill {\widehat{u}}_y=-{\tilde{k}}_x\frac{k_0}{k^2}·\end{array}$$(18)

The perturbation of pressure behaves in the following way:

$$\begin{array}{c}\hfill \widehat{W}=2\mathrm{i}\frac{k_0}{k^4}(k^2-q).\end{array}$$(19)

The solution (17)–(18) conserves the spanwise component of the vorticity (see Chagelishvili et al. 2003):

$$\begin{array}{c}\hfill {\widehat{\omega }}_z={\tilde{k}}_x{\widehat{u}}_y-{\widehat{u}}_x=-k_0,\end{array}$$(20)

which results in the transient growth of 2D SFH in the case of leading spirals, k_x < 0, for the quasi-Keplerian regime, q > 0, and in the case of trailing spirals, k_x > 0, for the cyclonic regime, q < 0. Indeed, at the instant of swing, t_s ≡ −k_x/q, the shearwise wavenumber vanishes, k̃_x = 0, and their growth factor defined as the ratio of current energy density of 2D SFH to its initial energy density reads

$$\begin{array}{c}\hfill g_s\equiv g{{|}_{t=t_s}\equiv ({\widehat{u}}_x^2+{\widehat{u}}_y^2)|}_{t=t_s}=\frac{k_0^2}{k^2(t=t_s)}=k_0^2.\end{array}$$(21)

For |k_x|≫1 g_s ≈ k_x², attaining an arbitrary high value for sufficiently wound spirals. The maximum possible g_s is limited by damping due to the microscopic viscosity of the fluid and will be evaluated in Sect. 4.

3.2. Dynamics of SFH with non-zero spanwise wavenumber

Let us consider SFH with small but non-zero |k_z|≪1⁴. Without loss of generality, we assume that k_z > 0 from now on. We also assume that initially SFH has a planar velocity field, that is, û_z(t = 0) = 0. Combining (10) and (12) we formulate the initial condition in the form

$$\begin{array}{c}\hfill {\widehat{u}}_x=k_0^{-1},{\widehat{\omega }}_x=k_z{k}_x{k}_0^{-1}.\end{array}$$(22)

It can be seen that Eq. (14) together with the initial condition Eq. (22) implies that ω̂_x ~ k_z, which makes the first term in the right-hand side (RHS) of Eq. (13) as small as ∼O(k_z²) comparing to the second term therein. Consequently, in order to obtain the solution for 3D SFH in the leading order in k_z ≪ 1, we drop the first term in the RHS of Eq. (13) going back to the 2D solution for û_x given by Eq. (15). Thus, for SFH weakly depending on the spanwise direction, û_x(t) is approximately given by Eq. (17).

Next, combining (10) and (12) we get an exact expression for û_z in terms of û_x and ω̂_xand :

$$\begin{array}{c}\hfill {\widehat{u}}_z=\frac{{\widehat{\omega }}_x-k_z{\tilde{k}}_x{\widehat{u}}_x}{1+k_z^2}·\end{array}$$(23)

Thereby, the term k_zû_z in the continuity Eq. (10) must be as small as ∼O(k_z²) in comparison with the other two terms there. Thus, we assume that relation (16) holds for 3D SFH and û_y is approximately given by its 2D expression, (18). Of course, this implies that pressure also approximately obeys the 2D solution (19).

As we see, spanwise dynamics of 3D SFH is separated from the swing amplification dynamics described in Sect. 3.1 and ω̂_x is obtained from Eq. (14) by integrating the expression (17) over the time,

$$\begin{array}{c}\hfill {\widehat{\omega }}_x=k_z{k}_0[\frac{2-q}{q}(arctan{\tilde{k}}_x-arctan{k}_x)+\frac{k_x}{k_0^2}]·\end{array}$$(24)

As t → ∞ and, respectively, |k̃_x| → ∞, the shearwise component of the vorticity perturbation tends to a constant non-zero value. We refer to this asymptote of SFH as a “plateau” of SFH. We conclude that along with the existing “planar” eddies, a set of other vortices associated with this non-zero limit of ω_x is generated as a “byproduct” of swing amplification. The latter will be referred to as “cross-rolls” in this work. The cross-rolls are directed across the stream in contrast to the streamwise rolls produced by the anti-lift-up mechanism in the flow on the Rayleigh line (see Fig. 9 and the details in Appendix A below). An important feature of the cross-rolls is that the spanwise velocity perturbation gains a non-zero value at t → ∞, whereas the velocity perturbation components associated with the planar eddies vanish like û_x ∼ O(t⁻²) and û_y ∼ O(t⁻¹). Equations (17) and (24) together with the relation (23) taken in the leading order in k_z ≪ 1 yield⁵

$$\begin{array}{c}\hfill {\widehat{u}}_z=k_z{k}_0[\frac{2-q}{q}(arctan{\stackrel{\sim }{k}}_x-arctan{k}_x)+\frac{k_x}{k_0^2}-\frac{{\tilde{k}}_x}{k^2}]·\end{array}$$(25)

Thus, Eq. (25) leads to the following asymptote of û_z at t → ∞,

$$\begin{array}{c}\hfill {\widehat{u}}_z\to \mathrm{sgn}(q)\pi k_z|k_x|\frac{2-q}{q}·\end{array}$$(26)

As soon as the shearing spiral is wound up by the shear after the instant of swing, the ratio |û_z/û_y| increases, which implies that the cross-section of cross-rolls becomes elongated in a spanwise direction. At the same time, cross-rolls shrink along their axes, as far as the shearwise wavenumber increases. This, in turn, makes a streamwise component of the vorticity perturbation, ω_y, grow linearly with time. Since

$$\begin{array}{c}\hfill {\widehat{\omega }}_y=k_z{\widehat{u}}_x-{\tilde{k}}_x{\widehat{u}}_z,\end{array}$$(27)

$$\begin{array}{c}\hfill {\widehat{\omega }}_y\to -\mathrm{sgn}(q)\pi k_z|k_x|(2-q)t.\end{array}$$(28)

Equation (26) demonstrates the basic features of cross-rolls production. Clearly, their amplitude is proportional to pressure perturbation spanwise gradient at the instant of swing, which is ∝ k_z|k_x| (cf. Eq. (19)) at k̃_x = 0. Also, it decreases as the shear rate goes from q = −∞ to q = 2 and vanishes exactly at q = 2. Those features will be illustrated below in Fig. 1, where the growth factors of shearing spirals with non-zero spanwise number will be plotted for different shear rates.

3.3. Restriction on value of spanwise wavenumber

Now we are a step away from estimating the maximum value of k_z = k_{z max}, which restricts the validity of the solution for 3D SFH obtained above in the leading order in k_z ≪ 1. At sufficiently large k_z, the swing amplification can no longer be separated from the spanwise dynamics, which starts to destruct the shearing spiral in the xy-plane. Moreover, comparing the first term in the RHS of Eq. (13) to the second term therein, we find that these terms always become of the same order at a sufficiently large time for any small value of k_z provided that the first term contains ω̂_x given by Eq. (24), while the second term contains û_x given by Eq. (17). This takes place since û_x decreases as O(t⁻²) long after the instant of swing. We must consequently provide a restriction on k_z estimating a higher order correction to the solution given in Sects. 3.1 and 3.2 at the selected period of time.

For that, let us assume that with the account of the first term in the RHS of Eq. (13), the 2D solution (17) acquires a small correction $k_z^2{\widehat{u}}_x\prime$. Substituting

$$\begin{array}{c}\hfill {\widehat{u}}_x=\frac{k_0}{k^2}+k_z^2{\widehat{u}}_x^{\prime }\end{array}$$

into Eq. (13), we derive the following equation for û′_x

$$\begin{array}{c}\hfill \frac{\mathrm{d}}{\mathrm{d}{\tilde{k}}_x}(k^2{\widehat{u}}_x^{\prime })=-2k_0[\frac{2-q}{q^2}(arctan{\tilde{k}}_x-arctan{k}_x)+\frac{k_x}{q{k}_0^2}-\frac{{\tilde{k}}_x}{k^4}]·\end{array}$$(29)

Given the zero initial condition û′_x(t = 0) = 0, we find that

$$\begin{array}{c}\hfill k^2{\widehat{u}}_x^{\prime }=-\frac{2k_x}{q{k}_0}({\tilde{k}}_x-k_x)+\frac{2k_0(2-q)}{q^2}[{\tilde{k}}_x(arctan{k}_x-\\ \hfill arctan{\tilde{k}}_x)+ln\left(\frac{k}{|k_0|}\right)]+\frac{1}{|k_0|}-\frac{|k_0|}{k^2}·\end{array}$$(30)

We estimate k_{z max} from the assumption that it corresponds to a constant ratio of the amplitude of $k_z^2{\widehat{u}}_x\prime$ and the amplitude of leading order û_x given by Eq. (17)

$$\begin{array}{c}\hfill k_{z\text{max}}^2\frac{|k^2{\widehat{u}}_x^{\prime }|}{|k_x|}{|}_{{\tilde{k}}_x={\tilde{k}}_{xmax}}=ϵ\ll 1,\end{array}$$(31)

where k̃_{x max} is taken at some time after the instant of swing and corresponds to a certain location of SFH on its plateau. We note that the higher |k̃_{x max}| is, the smaller is k_{z max} and the longer is the duration of our analytical approximation. Most naturally, k̃_{x max} can be evaluated at the latest extremum of û_z(t), which emerges as one incorporates viscous damping into the problem.

4. Taking viscous damping into account

Primarily, viscous force confines the swing amplification of SFH as the initial winding of the spiral and so the duration of 2D transient growth become limited by viscous dissipation. Thus, the largest value of g_s corresponds to some k_x = k_{x max}. Additionally, viscous force causes damping of vertical motions and the asymptote (26) is replaced by extremum of û_z at some time after the instant of swing. Below we identify this time with k̃_x = k̃_{x max} used in Eq. (31) (see Sect. 5).

To determine both k_{x, max} and k̃_{x max} rigorously, let us suppose that û_x, ω̂_x is the solution of Eqs. (13)–(14) obtained in the inviscid limit R → ∞. Then, the solution of the same equations taken with the account of non-zero viscosity R < ∞ is represented by

$$\begin{array}{c}\hfill {\widehat{u}}_{x,\nu }={\widehat{u}}_x{\mathrm{e}}^{-\gamma },\\ \hfill {\widehat{\omega }}_{x,\nu }={\widehat{\omega }}_x{\mathrm{e}}^{-\gamma },\end{array}$$

where

$$\begin{array}{c}\hfill \gamma =\frac{1}{\mathit{qR}}[(1+k_z^2)({\tilde{k}}_x-k_x)+\frac{1}{3}({\tilde{k}}_x^3-k_x^3)].\end{array}$$(32)

For simplicity below we neglect by $k_z^2$ in comparison with unity in (32) since k_z ≲ 1 as far as q ϵ ≲ 1, which is assessed from the condition (31).

The spanwise velocity perturbation is modified in the same way:

$$\begin{array}{c}\hfill {\widehat{u}}_{z,\nu }={\widehat{u}}_z{\mathrm{e}}^{-\gamma },\end{array}$$(33)

where û_z is given by Eq. (25). This allows one to obtain k_{x, max} and k̃_{x max} jointly as the roots of the set of equations

$$\begin{array}{c}\hfill \frac{\partial {\widehat{u}}_{z,\nu }}{\partial k_x}=0,\end{array}$$(34)

$$\begin{array}{c}\hfill \frac{\partial {\widehat{u}}_{z,\nu }}{\partial {\tilde{k}}_x}=0,\end{array}$$(35)

provided that the initial SFH is assumed to be the leading spiral (k_{x max} < 0 and k̃_{x max} > 0) in the quasi-Keplerian flow (q > 0) and, vice versa, the trailing spiral (k_{x max} > 0 and k̃_{x max} < 0) in cyclonic flow (q < 0).

4.1. The case of tightly wound SFH

In the limit of |k_{x max}|≫1 and |k̃_{x max}| ≫ 1, which is valid for sufficiently high R, the approximate relations

$$\begin{array}{c}\hfill arctan\{k_{x\text{max}},{\tilde{k}}_{x\text{max}}\}\approx \mathrm{sgn}(q)\pi /2-\{k_{x\text{max}}^{-1},{\tilde{k}}_{x\text{max}}^{-1}\}\end{array}$$

can be substituted into Eq. (25). Then, Eq. (34) simply recovers an approximate value of k_{x max}, which follows from the estimation of Razdoburdin & Zhuravlev (2017) of maximum duration of the swing amplification t_max (see their Eq. (B2)):

$$\begin{array}{c}\hfill k_{x\text{max}}\approx -q{t}_{\text{max}}=-\mathrm{sgn}(q){(R|q|)}^{1/3}.\end{array}$$(36)

Thus, a very high transient growth of 2D SFH ∝ R^2/3, and so the amplitude of the cross-rolls, can be reached taking into account that the Reynolds number may be as huge as 10¹⁰ in astrophysical boundary layers and discs.

Furthermore, Eq. (35) yields the estimation for k̃_{x max},

$$\begin{array}{c}\hfill {\tilde{k}}_{x\text{max}}\approx \mathrm{sgn}(q){\left(\frac{2|q|R}{(2-q)\pi }\right)}^{1/4},\end{array}$$(37)

provided that R^1/2 ≫ 1.

5. Maximum growth factor of cross-rolls

The growth of the cross-rolls can be measured by energy density stored in vertical motion. In this case, since the initial condition (22) sets SFH with the unit energy and the planar velocity field, the growth factor of the cross-rolls is defined as

$$\begin{array}{c}\hfill g_z\equiv {\widehat{u}}_{z,\nu }^2.\end{array}$$(38)

Then, the quantity

$$\begin{array}{c}\hfill g_{z\text{max}}\equiv {\widehat{u}}_{z,\nu }^2(k_{x\text{max}},{\tilde{k}}_{x\text{max}},k_{z\text{max}})\end{array}$$(39)

can be considered as an estimate of the largest possible amplification of cross-rolls generated by shearing spirals. We note that g_{z max}/ϵ is the function of R and q. One may assume that violation of our analytical approximation discussed in Sect. 3.3 corresponds to the breakdown of swing amplification being a unique mechanism generating the cross-rolls. Of course, an exact k_{z max} corresponding to the cross-rolls of the largest amplitude can be obtained solving the full set of Eqs. (13)–(14) without an expansion over k_z. Nevertheless, it is known that k_{z max} should be at least k_{z max} ≲ |q/(2 − q)|⁶, as the parameter β ≳ 1 introduced by Balbus & Hawley (2006) in their Eq. (39) substantially weakens the swing amplification and, consequently, the cross-rolls amplitude. In this study, we intend to draw main conclusions assuming that Eq. (39) gives a reasonable estimate of the largest cross-rolls generated by shearing spirals, whereas their exact identification from a strict solution of (13)–(14), as well as determination of the corresponding k_{z max}, is a subject for future studies with the use of an optimisation technique.

Below we determine the largest growth factor of cross-rolls using accurate and approximate procedures.

• i)

  An accurate determination of g_{z max}(q, R) is done by solving Eqs. (34)–(35) numerically, which subsequently provides us with k_{z max} from the condition (31). Then, g_{z max} is obtained from Eq. (33) with k_x = k_{x max}, k̃_x = k̃_{x max}, and k_z = k_{z max}.

• ii)

  In the case of sufficiently high R_{gz max} can be estimated analytically. Indeed, as the conditions k_{x max} ≫ 1 and k̃_{x max} ≫ 1 are fulfilled, Eq. (30) simplifies to

  $$\begin{array}{c}\hfill k^2{\widehat{u}}_x^{\prime }=-\mathrm{sgn}(q)\frac{2(2-q)}{q^2}|k_x|{\tilde{k}}_x\end{array}$$

  and yields the following estimation for k_{z max}

  $$\begin{array}{c}\hfill k_{z\text{max}}\approx {\left(\frac{{|q|}^7}{2^5{\pi }^3{(2-q)}^{3}R}\right)}^{1/8}{ϵ}^{1/2}.\end{array}$$(40)

  Finally, in order to estimate the largest growth factor of cross-rolls in this case, one must substitute Eqs. (36), (37), and (40) into

  $$\begin{array}{c}\hfill {\widehat{u}}_{z,\nu }\approx \pi \mathrm{sgn}(q)k_z|k_x|\frac{2-q}{q}{\mathrm{e}}^{-\frac{{\tilde{k}}_x^3-k_x^3}{3qR}}.\end{array}$$(41)

Now, setting g_{z max}/ϵ obtained according to the procedure (i) or (ii) to a constant value, we find R(q) as a root of Eq. (39). This also implies that k_{x max}(q), k̃_{x max}(q), and $k_{z\text{max}}(q)/\sqrt{ϵ}$ have been determined, so that we are ready to plot the growth factors of the particular SFH corresponding to constant g_{z max}/ϵ for several values of q (see Fig. 1). These growth factors have to be regarded as functions of time or, equivalently, k̃_x: so g(k̃_x) is defined by Eq. (21) and g_z(k̃_x) is defined by Eq. (38). As one can see on the top panel in Fig. 1, all curves of g_z(k̃_x) attain the same maximum value g_{z max}/ϵ, which occurs after the instant of swing at negative and positive k̃_x = k̃_{x max} for cyclonic and quasi-Keplerian regimes, respectively. If it were not for viscous damping, g_z would tend to the plateau, which has been found previously for the inviscid solution (see Eq. (26)). While comparing two panels in Fig. 1, we see that the ratio (g_s/g_{z max})ϵ, where g_s is a factor of 2D swing amplification of SFH, dramatically (and monotonically) increases as one goes from large cyclonic shear rates across the solid-body line up to the marginal case of the Rayleigh line q = 2. For example, at the Keplerian shear rate this ratio attains almost 10². Thus, the cross-rolls of the same intensity are generated by SFH with much larger |k_x| as one approaches the Rayleigh line. The larger is |k_x|, the tighter the winding of the initial shearing spiral is, the more it is susceptible to viscous dissipation at a typical time ∝ |k_x|⁻². Consequently, the closer one approaches the Rayleigh line, as starting from q = −∞, the higher the Reynolds number required in order to generate cross-rolls of the same intensity. This is reflected in the fit of our simulations of the subcritical transition to turbulence in the cyclonic regime (see Sect. 6 and solid curve on the top panel in Fig. 6) as well as in our prediction of the subcritical transition to turbulence in the quasi-Keplerian regime (see Sect. 7 and curve in Fig. 8).

Next, Fig. 2 illustrates the immediate cause of the suppression of the cross-rolls generated in the low-shear (|q|∼1) cyclonic regime and even more in the quasi-Keplerian regime in comparison with the high-shear (|q|≫1) cyclonic regime. There are the profiles of pressure perturbation amplitude as a function of k̃_x plotted for shearing spirals shown in Fig. 1. For SFH these profiles also represent the behaviour of spanwise acceleration of fluid elements in the perturbed flow. At first, as the shear rate goes from −∞ up to 1/2, the net acceleration impulse in the spanwise direction decreases just because |Ŵ| in the vicinity of the instant of swing attains a smaller value. Next, at q = 1/2 the point |Ŵ|(k̃_x = 0) transforms from maximum to minimum and, further, for q > 1, |Ŵ| changes sign in the vicinity of the instant of swing. Eventually, as q → 2, the net acceleration of fluid elements in positive and negative spanwise directions balance each other, so that one is left with planar velocity perturbation as the shearing spiral is wound up by the shear after the instant of swing. Certainly, an additional dip in growth of the cross-rolls, which must occur in the vicinity of the solid-body line q = 0, is imposed on the described picture. Therefore, there must exist a local maximum of growth of the cross-rolls generated by the initial SFH with constant k_x and k_z in the quasi-Keplerian regime. In Sect. 7 we show that this local maximum comes out through the local minimum of the transition Reynolds number occurring near the q = 1/2.

5.1. Cross-rolls in comparison with rolls and streaks

Let us estimate the growth factor of rolls generated by the anti-lift-up effect and streaks generated by the lift-up effect in centrifugally stable rotating flow. This is to be done in order to evaluate the role of those known mechanisms relative to the generation of cross-rolls discussed above.

Both lift-up and anti-lift-up mechanisms describe the transient growth of axisymmetric perturbations. Consequently, we set k_y = 0 in Eqs. (7)–(10) and obtain the set of equations for û_x and ω̂_x

$$\begin{array}{c}\hfill \frac{\mathrm{d}{\widehat{u}}_x}{\mathrm{d}t}=-2\frac{k_z}{k_x^2+k_z^2}{\widehat{\omega }}_x,\end{array}$$(42)

$$\begin{array}{c}\hfill \frac{\mathrm{d}{\widehat{\omega }}_x}{\mathrm{d}t}=(2-q)k_z{u}_x\end{array}$$(43)

in the inviscid limit R → ∞, which yields the equation for û_x

$$\begin{array}{c}\hfill \frac{{\mathrm{d}}^2{\widehat{u}}_x}{{\mathrm{d}}^{2}t}+{\sigma }^2{\widehat{u}}_x=0,\end{array}$$(44)

where

$$\begin{array}{c}\hfill {\sigma }^2\equiv \frac{2(2-q)}{1+{\alpha }^2}\end{array}$$

with α ≡ k_x/k_z.

The solution of Eq. (44) reads

$$\begin{array}{c}\hfill {\widehat{u}}_x=Asin(\sigma t+\phi ),\end{array}$$(45)

where amplitude A and phase φ are to be determined from initial conditions.

Further, ω̂_x is obtained with the help of Eq. (43)

$$\begin{array}{c}\hfill {\widehat{\omega }}_x=-A(2-q)k_z{\sigma }^{-1}cos(\sigma t+\phi ),\end{array}$$(46)

whereas

$$\begin{array}{c}\hfill {\widehat{u}}_y=A(2-q){\sigma }^{-1}cos(\sigma t+\phi ),\end{array}$$(47)

since ω̂_x = k_zû_y for axisymmetric perturbations. Finally,

$$\begin{array}{c}\hfill {\widehat{u}}_z=-\alpha {\widehat{u}}_x=-A\alpha sin(\sigma t+\phi ).\end{array}$$(48)

5.1.1. Growth of rolls

The lift-up mechanism provides the transient growth of streamwise rolls, that is, initial vortices with û_y = 0. This implies that ${\widehat{u}}_x=\sqrt{1/(1+{\alpha }^2)}$ and ${\widehat{u}}_z=-\sqrt{{\alpha }^2/(1+{\alpha }^2)}$ for perturbations with initial ${\widehat{u}}_x^2+{\widehat{u}}_z^2=1$. Thus, $A=\sqrt{1/(1+{\alpha }^2)}$ and φ = π/2.

Consequently, the growth factor of rolls defined as ḡ ≡ g + g_z reads

$$\begin{array}{c}\hfill \overline{g}=\frac{2-q}{2}{sin}^2\sigma t+{cos}^2\sigma t.\end{array}$$(49)

The corresponding maximum growth factor reads

$$\begin{array}{c}\hfill {\overline{g}}_{\text{max}}=\frac{2-q}{2}\end{array}$$(50)

and becomes unlimited as q → −∞, that is, in the non-rotating shear flow.

5.1.2. Growth of streaks

The anti-lift-up mechanism provides the transient growth of streamwise streaks, that is, initial vortices with û_x = 0, û_y = 1, û_z = 0. Thus, in this case $A=\sqrt{2/\left[(2-q)({\alpha }^2+1)\right]}$, φ = 0.

The growth factor of streamwise streaks is

$$\begin{array}{c}\hfill \overline{g}=\frac{2}{2-q}{sin}^2\sigma t+{cos}^2\sigma t.\end{array}$$(51)

The corresponding maximum growth factor reads

$$\begin{array}{c}\hfill {\overline{g}}_{\text{max}}=\frac{2}{2-q}\end{array}$$(52)

and becomes unlimited as q → 2, on the Rayleigh line.

Hence, the maximum transient growth factor attained through either the lift-up or anti-lift-up mechanism is independent of SFH wavenumbers. Since there is no low limit on their absolute values in the boundless shear, either the rolls or the streaks can be amplified up to ḡ_max given, respectively, by Eqs. (50) and (52) at any R. At the same time, the cross-rolls cannot be amplified stronger than the value estimated by definition (39) with û_{z, v} substituted from Eq. (33). One can see that additionally to a constraint coming from the damping of the swing amplification, which is $g_{\text{max}}\approx k_{x\hspace{0.17em}\text{max}}^2$, the growth of the cross-rolls is reduced by a factor of $\sim ((2-q)/q)\hspace{0.17em}{k}_z^2$ comparing to ḡ_max. We conclude that strictly in the absence of rotation or, alternatively, on the Rayleigh line, cross-rolls are generally weaker than streaks or, respectively, rolls. However, when −∞ < q < 2, the transient growth of both rolls and streaks is highly suppressed. For example, the growth of rolls and streaks drops below ḡ_max ∼ 20, which is equal to g_{z max} represented in Fig. 1 at ϵ = 0.1, for q ≳ −40 and q ≲ 1.9, respectively. That is why in both the cyclonic and the quasi-Keplerian regime of homogeneous boundless shear flow, cross-rolls are the only high amplitude 3D perturbations that can be generated via the transient growth mechanisms.

6. Transition to turbulence: Cyclonic regime

In this section we consider 3D non-linear dynamics of perturbations in the cyclonic regime of uniform and boundless rotating shear flow. Our goal is to recognize the subcritical transition to turbulence towards smaller |q| compared to the previous studies and to obtain the transition Reynolds number, R_T. This is done by performing hydrodynamical numerical simulations of the evolution of initial finite amplitude perturbations in the shearing box approximation. Importantly, we finally check how the cross-rolls maximum growth factor estimated above behaves along the dependence R_T(q).

6.1. Model for non-linear dynamics of perturbations and numerical details

Numerical simulations were carried out using the open-source code ATHENA for solving continuity and Euler equations (see Stone et al. (2008) and Stone & Gardiner (2010)),

$$\begin{array}{c}\hfill \frac{\partial \rho }{\partial t}+\mathrm{\nabla }·(\rho \mathbf{v})=0,\end{array}$$(53)

$$\begin{array}{cc}\hfill \frac{\partial (\rho \mathbf{v})}{\partial t}+\mathrm{\nabla }·(\rho \mathbf{v}\mathbf{v}& +P)=2q\rho {\mathrm{\Omega }}_0^2\mathbf{x}-2{\mathrm{\Omega }}_0\times \rho \mathbf{v}+\hfill \\ \hfill & \mathrm{\nabla }·\left(\rho {\nu }_0\mathrm{\nabla }\mathbf{v}\right)+\mathrm{\nabla }·(\frac{1}{3}\rho \nu \mathrm{\nabla }·\mathbf{v}),\hfill \end{array}$$(54)

supplemented by

$$\begin{array}{c}\hfill P=c_s^2\rho ,\end{array}$$(55)

where x, y, z, and t are the Cartesian coordinates and time used in Eqs. (1)–(4); v, P, and ρ are, respectively, velocity, pressure, and density of the perturbed flow. In Eqs. (53)–(55) kinematic viscosity ν and isothermal speed of sound c_s are constants and the full non-linear dynamics of perturbations is considered in a box of size L_x, L_y, and L_z along x, y, and z directions, respectively. The usual periodic and shearing-periodic boundary conditions are imposed, respectively, along y, z, and x directions at each side of the box (see e.g. Hawley et al. 1995). We note that everywhere below L_x = L_y = 1 but L_z ≥ L_x, L_y.

In all our simulations

$$\begin{array}{c}\hfill \text{max}\{L_{x,y,z}\}\lesssim H,\end{array}$$(56)

where H ≡ c_s/Ω₀. According to hydrodynamical equilibrium in accretion disc boundary layer, H is a typical scaleheight of the flow. The condition (56) ensures that we work with vortical dynamics of perturbations as long as vortical initial perturbations are imposed.

It is suitable to parametrize the problem by two dimensionless numbers: Mach number

$$\begin{array}{c}\hfill M=\frac{|q|{\mathrm{\Omega }}_0{L}_y}{c_{\mathrm{s}}}\end{array}$$(57)

and Reynolds number

$$\begin{array}{c}\hfill R_{\mathrm{nl}}=\frac{{\mathrm{\Omega }}_0{L}_y^2}{\nu }.\end{array}$$(58)

We note that R_nl is different from R used above to describe the dynamics of the cross-rolls in the boundless flow. The relation between both Reynolds numbers can be derived along the following line. Suppose that there is SFH with dimensional shearwise wavenumber K_y. In the box it takes values

$$\begin{array}{c}\hfill K_y=\frac{2\pi }{L_y}n,\end{array}$$(59)

where n is a natural number. On the other hand,

$$\begin{array}{c}\hfill K_y=L^{-1},\end{array}$$

as far as in the analytical linear model of transient growth k_y = 1. Using the definitions (11) and (58), we find

$$\begin{array}{c}\hfill \frac{R_{\mathrm{nl}}}{R}=4{\pi }^2{n}^2\end{array}$$(60)

for SFH with shearwise wavelength n times shorter than the shearwise box size.

In all simulations, the density of the unperturbed flow is set to ρ₀ = 1 and the box rotation frequency Ω₀ = 0.001. Setting q, M, and R_nl allows us to specify an unperturbed value of c_s and value of ν, as well as pressure in the unperturbed flow, P₀, through the equation of state (55). Numerical resolution is chosen to keep a cubic form of a grid-cell, that is, the number of nodes in the x and y directions is equal to each other, N_x = N_y ≡ N, and the number of nodes in spanwise direction is N_z = l_zN, where the box ratio l_z ≡ L_z/L_y. We used N = 64 once in order to compare our results with those of LL05 for the case of q = −10. In all other simulations, N = 128.

6.2. Indication of transition to turbulence

In order to recognize turbulence, we employ a dimensionless angular momentum flux

$$\begin{array}{c}\hfill \alpha =-\frac{<\rho v_x{v}_y^{\prime }>}{{\rho }_0{c}_s^2}\end{array}$$(61)

arising due to velocity perturbation to the unperturbed value v_y = −qΩ₀x. In Eq. (61) v_x and v_y′ we introduce perturbations of shearwise and streamwise velocity components, respectively, and the brackets denote spatial averaging.

Angular momentum flux arising due to microscopic viscosity in unperturbed flow is

$$\begin{array}{c}\hfill {\alpha }_{\nu }=\frac{{\nu }_0{\rho }_{0}rd\mathrm{\Omega }/dr}{P_0}=-\frac{M^2}{\mathit{qR}}·\end{array}$$(62)

We compare α(t) with α_ν in the course of the non-linear evolution of perturbations and interpret the instant enhanced transport of angular momentum α(t)> α_ν as the existence of turbulence at time t.

Further, in order to check the transition to turbulence in the flow for given R_nl, the time-averaged variant of angular momentum flux is also used:

$$\begin{array}{c}\hfill \tilde{\alpha }\equiv \frac{1}{t_2-t_1}\underset{t_1}{\overset{t_2}{\int }}\alpha \mathrm{d}t,\end{array}$$(63)

where |q|t₁ = 500 Ω₀⁻¹ and |q|t₂ = 1000 Ω₀⁻¹. The transition to turbulence is defined by the condition

$$\begin{array}{c}\hfill \tilde{\alpha }>{\alpha }_{\nu },\end{array}$$(64)

which implies that a self-sustained turbulent state exists for a sufficiently long time.

It was found that the evolution of initial perturbations crucially depends on R_nl. If R_nl exceeds some threshold value, which we refer to as transition Reynolds number, R_T, initial perturbations give rise to a self-sustained non-stationary solution. In the opposite case, R_nl < R_T, initial perturbations fade out (see Fig. 3). In order to determine R_T, we perform a set of simulations on a uniform grid of Reynolds number values. The transition Reynolds number is defined as R_T = (R_turb + R_damp)/2, where R_turb equals the lowest R_nl at which the flow demonstrates the transition to turbulence and R_damp equals to the highest R_nl at which we do not observe a long-lasting turbulent state of the flow and, consequently, the condition (64) is not satisfied. As we see in Fig. 4, the pass through R_T is associated with a significant jump of α̃.

6.3. Initial conditions

At first, we check whether the subcritical transition to turbulence we observe in our numerical experiment is in accordance with the results of LL05. Performing simulations in a cubic box, they found R_T = 1200 at q = −10 for our definition of R_nl (LL05 measure time in units of the inversed shear rate, which makes their Reynolds number R_nl multiplied by q).

Generally, this is a feature of the subcritical turbulisation that R_T depends on the shape and amplitude of the initial perturbations. For example, this was shown experimentally by Darbyshire & Mullin (1995) and later numerically by Faisst & Eckhardt (2004) for pipe flow. One may assume that for any particular values of q and M there are optimal initial perturbations, which provide the lowest possible value of R_T. These optimals can be found by a variational approach to the non-linear Cauchy problem (see Cherubini et al. 2010 and Pringle & Kerswell 2010 for such an example). Unfortunately, this method is not easy to implement and it is time consuming. In this work we restrict a comparison with LL05 to a single variant of initial perturbations. Also, we specify the same resolution as LL05, N_x = N_y = N_z = N ≡ 64 using the cubic box.

We choose an initial condition in the form of the vortex

$$\begin{array}{c}\hfill v_x=0,\end{array}$$(65)

$$\begin{array}{c}\hfill v_y^{\prime }=-\frac{A{c}_{\mathrm{s}}\sqrt{2}{K}_z}{\sqrt{K_z^2+K_y^2}}sin(K_{y}y+K_{z}z),\end{array}$$(66)

$$\begin{array}{c}\hfill v_z=\frac{A{c}_{\mathrm{s}}\sqrt{2}{K}_y}{\sqrt{K_z^2+K_y^2}}sin(K_{y}y+K_{z}z),\end{array}$$(67)

where K_y and K_z are shearwise and spanwise dimensional wavenumbers of vortex, and the dimensionless constant A specifies the initial specific kinetic energy of perturbations, which equals A²/2.

It was revealed that for l_z = 1, A = 1, and K_y = K_z = 2π the initial vortex (65)–(67) leads to self-sustained turbulence at R_turb = 1400 and damps at R_damp = 1200 (see Figs. 3 and 4. Thus, R_T = 1300, which is close to the value obtained by LL05. We additionally checked that both the decrease and the increase of the amplitude of initial vortex, respectively, up to A = 0.5 and A = 2 lead to the increase of the transition Reynolds number, respectively, up to R_T = 2600 and to R_T = 1500. Until the end of this section, we use the initial perturbations in the form (65)–(67) with A = 1, K_y = 2π, and K_z = 2π/l_z.

6.4. Results

The picture of the subcritical transition to turbulence we extract from our simulations is presented in Table 1. For the first time (with regards to the model of unbounded uniform hydrodynamical flow considered here), we use not only the cubic box, but also a tall box with the box ratio l_z > 1. The box acts as a filter with respect to perturbations of different wavelengths. It prevents the existence of both small-scale perturbations with a characteristic scale comparable and less than the size of a numerical grid-cell, and large-scale perturbations with a characteristic scale larger than the size of the box in any direction. The tall box allows perturbation harmonics with K_z/K_y = l_z⁻¹ < 1 to take part in the non-linear evolution of the perturbed flow and, probably, in the subcritical transition to turbulence. In particular, this applies to shearing spirals that effectively generate the cross-rolls. Indeed, employing Eq. (40) one estimates K_z/K_y < k_{z max}/ϵ^1/2 ≈ 0.8 for q = −10 and R_T = 1300. This ratio becomes smaller as one shifts to a smaller shear rate and higher Reynolds number. For example, at q = −6.67 and R_T = 10⁴, which is the smallest shear rate examined by LL05 for the subcritical transition, its value is less than approximately 0.5. At the same time, the small size of the grid-cell is always required in order for the tightly wound spirals exhibiting large transient growth and generating high-amplitude cross-rolls to be involved into simulations.

The subcritical transition in the cubic box has been observed up to q = −4 (see Table 1 and Fig. 6). This is a higher value compared to q = −6.67 reached by LL05, which is due to a higher numerical resolution used in our simulations. Indeed, assuming that turbulisation, at least indirectly, is provided by the transient growth of shearing spirals and that the smallest shearwise scalelength of perturbations resolved in the numerical scheme comprises four grid-cells, one can estimate the largest achievable R_nl provided that the shearing spirals subject to the largest swing amplification are resolved by the numerical scheme. The latter is done by Eq. (36),

$$\begin{array}{c}\hfill R_{\mathrm{nl}}\approx 4{\pi }^2{\left(\frac{N}{4}\right)}^3{|q|}^{-1}.\end{array}$$(68)

For N = 128 and q = −4 Eq. (68) yields R_nl ≈ 330 000. This is close to the R_T obtained for q = −4 (see Table 1), which naturally explains why we did not see the subcritical transition at smaller |q|.

However, once we fix q = −4 and perform a set of simulations in the tall box, we find that R_T falls by a factor of ∼20 as l_z increases from 1.0 up to 1.5 ÷ 2 (see Fig. 5). We check that for l_z > 2.0, R_T tends to horizontal asymptote. Further, we managed to find the subcritical transition at even higher q = −3 in the box with the same resolution N = 128 but with the box ratio l_z = 2. We regard this result as evidence for the crucial role of harmonics of perturbations with certain 0 < (K_z/K_y)_T < 1⁷. It is plausible that R_T(l_z) acquires an asymptotic value as soon as l_z⁻¹ ≤ (K_z/K_y)_T. Such a harmonics already exists in the cubic box, however, it corresponds to n > 1 (see Eq. (59)), and, consequently, to the effective R smaller by a factor of n² (see Eq. (60). Harmonics of a columnar shape corresponding to the highest effective R (i.e. to the largest shearwise wavelength equal to L_y), which implies the highest transient growth, emerges in the tall box only.

The question arises, what can additionally testify in favour of cross-rolls as moderators of subcritical transition? Trying to address this issue, we obtain R_T in the tall box with l_z = 2.0 for various q (see Table 1 and Fig. 6). It was additionally checked that the box ratio l_z sufficient for R_T(l_z) to approach the horizontal asymptote becomes smaller as one goes to higher |q|. This guarantees that subcritical transition observed in the tall box with l_z = 2.0 occurred in the presence of harmonics with (K_z/K_y)_T in the range of the shear rate from q = −15 up to q = −3 covered by simulations. It was revealed that the difference between R_T obtained in cubic and tall boxes increases as |q| becomes smaller. The tall box result R_T = 16 250 at q = −4 allows us to evaluate the maximum growth factor of cross-rolls according to the procedure (i) formulated in Sect. 5,

$$\begin{array}{c}\hfill g_{z\text{max}}/ϵ\approx 150.\end{array}$$(69)

Previous inspection of the maximum optimal growth factor, G_max, corresponding to transition to turbulence at R_Ω = 0 and R_Ω = −1 in the framework of the shearing sheet model (see e.g. Mukhopadhyay et al. 2005) has shown that G_max ∼ 150 in those cases, which coincides with Eq. (69) up to the smallness factor. Using procedure (i) further in order to look for the root of Eq. (39), we obtain the curve of constant (69) on the plane (R_T, q) (see the solid line in Fig. 6). Additionally, the dashed line in Fig. 6 represents the corresponding approximate R_nl(q) obtained according to procedure (ii) formulated in Sect. 5. A notable accordance is found between the growth of cross-rolls and the subcritical transition: namely, the largest growth factor of cross-rolls, as estimated in the main order in small k_z, remains almost constant along the transitional path, R_T(q), deduced from our non-linear stability analysis. Below we will refer to Eq. (69) as the threshold value of g_{z max}/ϵ. If there is a common non-linear mechanism sustaining the subcritical turbulence in the uniform unbounded rotating shear flow of the cyclonic type, then it is natural to assume that it is triggered at the threshold value of the transient growth factor of relevant perturbations, at most weakly changing with the shear rate. Here we find such evidence in favour of the cross-rolls and their maximum growth factor, g_{z max}, derived in the limit of small k_z. It can be suggested that the subcritical transition, which occurs at asymptotic value of R_T(lz) in the tall box simulations, takes place due to a positive non-linear feedback from the harmonics of cross-rolls with k_z ≃ (K_z/K_y)_T. However, the ability of cross-rolls with g_z ≳ g_{z max} to generate secondary shearing spirals subject to swing amplification sufficient to sustain turbulence remains to be shown.

We note that as q ≲ −5, the approximate analytical curve in Fig. 6 is no longer in accordance with the curve obtained accurately. The reason for that is shown in Fig. 7, where k_{x max}, k̃_{x max}, and $k_{z\text{max}}/\sqrt{ϵ}$ along the corresponding curves of R_nl(q) are plotted. Clearly, the limit of tightly wound spiral, |k_{x max}|, |k̃_{x max}| ≫ 1, is valid as far as |q|∼1 and breaks when q ≲ −5. For |q|≫1 R_nl(q) corresponding to constant (69) becomes so small that spirals are no longer tightly wound (see Eqs. (36) and (37)). At the same time, we find that the first point of the subcritical transition obtained by LL05 in the presence of rotation, R_Ω > 0, which corresponds to q ≈ −33, is also in reasonable accordance with the threshold value (69). Along with the comparison of rolls and cross-rolls growth factors with each other at |q|≫1 (see the end of Sect. 5.1), this argues in favour of the assumption that, at least, when q ≳ −35, SSP is already replaced by a different self-sustaining solution incorporating the growth of cross-rolls.

The dependence k_{z max}(q) obtained along the cross-rolls threshold growth factor (69) (see Fig. 7) confirms that the latter is exhibited by shearing spirals with large spanwise length scale. As long as |q| becomes smaller, the destructive role of the Coriolis force in the evolution of SFH with non-zero k_z increases as compared to the pressure gradients enabling their swing amplification. That is why k_{z max} decreases as one approaches the line of rigid rotation, which is seen in Fig. 7 (see also Eq. (40)). In turn, this naturally explains why the difference between R_T obtained in cubic and tall boxes increases for smaller |q|: SFH with n = 1 (see Eq. (59)), corresponding to the growth factor (69) is progressively suppressed by periodic boundary conditions imposed in the cubic box in the spanwise direction. We note that SFH that generate the cross-rolls of the largest amplitude will be highly columnar at |q| ≈1, as far as $k_{z\text{max}}/\sqrt{ϵ}$ is estimated to become less than 0.1. If subcritical turbulence in the cyclonic regime is sustained by the non-linear feedback from the cross-rolls, it is worth checking its existence at |q|∼1 via the DNS in a very tall box.

7. Transition to turbulence: Towards the Keplerian shear across the solid-body line

Let us suppose that cross-rolls are responsible for the subcritical transition in the whole range of centrifugally stable rotating shear flows. If so, one can assume that the corresponding threshold value of the cross-rolls largest growth factor, g_{z max}, found in the cyclonic regime (see the previous section) remains independent of the shear rate further in the quasi-Keplerian regime. Thus, we suppose that equality (69) holds for R_nl = R_T across the solid body line q = 0 up to the Rayleigh line q = 2. Since the shearing spirals, which produce the threshold growth factor, become increasingly wound up as one goes from q = −∞ to q = 2 (see the bottom panel in Fig. 1 and the top panel in Fig. 7), it is sufficient to use an approximate evaluation of g_{z max} provided by the procedure (ii) described in Sect. 5. So R_nl(q) corresponding to a constant threshold growth factor (69) extends to the quasi-Keplerian regime as shown by the curve in Fig. 8. This curve may be considered as the prediction of the subcritical transition to turbulence in homogeneous boundless quasi-Keplerian flow. Generally, the turbulence provided by the non-linear feedback from the cross-rolls is expected to occur at much higher Reynolds number ∼10⁷ ÷ 10⁹ in the quasi-Keplerian regime when compared to the cyclonic regime. On the one hand, this is consistent with negative results from all existing studies of the non-linear instability of the quasi-Keplerian regime, whereas on the other, even such a huge Reynolds number is feasible in astrophysical discs (see e.g. Balbus 2003). Equations (40) and (26) show that such an asymmetry of growth of the cross-rolls in the cyclonic and quasi-Keplerian regimes is controlled by the factor

$$\begin{array}{c}\hfill g_{z\text{max}}\propto {(2-q)}^{5/4},\end{array}$$

which stems from the suppression of fluid vertical motion in the course of the swing amplification of SFH (see comments to Fig. 2). Consequently, we predict an asymmetry of R_T(q) with respect to the line of rigid rotation: for the same absolute value of the shear rate, |q|, the subcritical transition in the quasi-Keplerian regime should occur at the substantially higher R_T rather than in the cyclonic regime.

Let us note that the enhanced non-linear stability of quasi-Keplerian flow follows from our DNS in the cyclonic regime together with the results of LL05 and Rincon et al. (2007) discussed in Introduction. Indeed, the curves on top panel in Fig. 6 indicate that R_T ≈ 3 × 10⁵ at q = −1.5. At the same time, from extrapolation of DNS by LL05 carried out at R_Ω < −1, see expressions in their Fig. 7, it is expected that R_T ≳ 10⁹ at q = 1.5. This independently shows an asymmetry of R_T(q) with respect to change of sign of q. However, the swing amplification itself produces maximum growth factor $\sim \hspace{0.17em}{k}_{z\hspace{0.17em}\text{max}}^2$ , which would yield a symmetric R_T ∝ |q|⁻¹ with respect to change from cyclonic to quasi-Keplerian shear rate.

It is noteworthy that prediction of R_T in Fig. 8 combined with extrapolations of DNS suggested by LL05 for quasi-Keplerian regime gives the most stable rotating flow located at super-Keplerian shear rate q ∼ 1.7 ÷ 1.8.

Further, since in the vicinity of the solid-body line R_T naturally tends to infinity (cf. Eqs. (41) and (36), (40) at |q|→0), a minimal R_T is predicted for the quasi-Keplerian regime. The corresponding most “favourable” value of q can be immediately estimated in the limit of tightly wound SFH. Using Eq. (41) in combination with Eqs. (36) and (40) and the condition R ≫ 1, we find that constant g_{z max} corresponds to

$$\begin{array}{c}\hfill {(R_{\mathrm{T}}q)}^{2/3}{(2-q)}^2=\mathrm{const},\end{array}$$

which yields minimal R_T at q = 1/2. So, R_T(q = 1/2) exceeds 10⁷, while the transition at the Keplerian shear, q = 3/2, may occur at R approaching ∼10⁸. Hence, the super-Keplerian flow 3/2 < q < 2 is predicted to be the most non-linearly stable throughout the whole range of centrifugally stable rotating shear flows. Notably, most of the efforts to detect turbulence in the quasi-Keplerian regime in laboratory experiments and DNS have been made at super-Keplerian shear rates. In contrast, our results suggest that one should extend the turbulent solutions obtained so far in the cyclonic regime to smaller |q| and then across the solid-body line, finally focusing the efforts on sub-Keplerian rotation. The effective numerical viscosity existing in DNS is approximately proportional to the distance between the nodes of the numerical grid. If so, the maximum Reynolds number accessible in DNS is ∝ N². Then, assuming that R_T ≈ 40 000 at q = −3, see Table 1, is the maximum accessible R_nl for N = 128 in tall box, we suppose that in order to be able to check hypothetical transition to hydrodynamic turbulence in quasi-Keplerian flow, at least for “favourable” q = 1/2, N ≳ 2000 is required. This is quite a high DNS resolution, however, it is not unreachable for the most powerful modern supercomputers.

8. Summary

Angular momentum transport in boundary layers around the accreting weakly magnetised stars may be produced by turbulent fluctuations. Here we suggest that subcritical turbulence in boundless and homogeneous rotating shear flow of the cyclonic type, which is a local representative of the boundary layer, is supplied with energy by nearly optimal shearing spirals that generate essentially 3D vortices of high amplitude as a “byproduct” of swing amplification. We call those 3D vortices “cross-rolls” in order to distinguish them from the known streamwise rolls generated via the anti-lift-up mechanism in a flow on the Rayleigh line. Such turbulence must be sustained by a novel type of positive non-linear feedback from the cross-rolls, which would recover the basin of growing shearing spirals. Furthermore, we suppose that this scenario also works in the quasi-Keplerian regime. These suggestions are supported by the following evidence.

• i)

  As well as known spanwise invariant (columnar) shearing spirals, shearing spirals with small but non-zero spanwise wavenumber are subject to swing amplification, which causes their transient growth almost up to the highest possible values corresponding to a growth factor much greater that unity. Thus, these vortices satisfy a necessary condition for the subcritical transition in both the cyclonic and quasi-Keplerian regimes. At the same time, it is known that all shearing harmonics with k_z ≳ 1 do not exhibit considerable transient growth since they drastically degenerate into epicyclic oscillations.

• ii)

  The substantial decrease of R_T and the existence of an asymptotic value of R_T(l_z) for large l_z found in the tall box simulations (see Fig. 5) indicates that these perturbations are crucial for the transition and sustenance of subcritical turbulence, at least in the cyclonic regime. This is the most striking evidence that the cross-rolls could be ingredients of positive non-linear feedback completing the bypass scenario in centrifugally stable rotating shear flows.

• iii)

  It is revealed that the maximum growth factor of cross-rolls decreases as the shear rate alters from q = −∞ through the solid-body line q = 0 up to q = 2 (see Fig. 1). Moreover, its constant value accords with the curve of R_T(q < 0) obtained in tall box simulations (see Fig. 6). We regard this value as the threshold growth factor, which triggers the non-linear sustenance of subcritical turbulence in the cyclonic regime. On the contrary, the spanwise invariant shearing spirals with k_z = 0 do not produce any 3D vortices and acquire a growth factor symmetric with respect to change between the cyclonic and quasi-Keplerian regimes at the same |q|. However, DNS of turbulence in the cyclonic regime performed by LL05 and in this work far beyond the zone of the lift-up effect show a strong variation of the 2D growth factor along the obtained R_T(q < 0).

• iv)

  The disappearance of the cross-rolls in the vicinity of q = 2 naturally explains the asymmetry of the transition to turbulence revealed in DNS close to R_Ω = −1 versus R_Ω = 0: the bypass scenario based on the transient growth of shearing spirals with small but non-zero k_z for q < 0 replaces SSP, retaining the non-linear instability of the cyclonic flow immediately after the lift-up ceases to operate, whereas a similar situation with the anti-lift-up mechanism is not realised in the case of quasi-Keplerian rotation.

• v)

  The cross-rolls are the only high-amplitude essentially 3D vortices that can be generated via the transient growth mechanism in the whole range of centrifugally stable shear rates, provided that the Reynolds number is sufficiently high. They create the vorticity (namely, its shearwise and streamwise components), which must be of crucial importance for 3D turbulence.

• vi)

  The transition to turbulence in the tall box is observed until smaller |q| than that in the cubic box at the same numerical resolution is obtained (see. Fig. 6). Generally, a sensitivity of the transition to the shape of a box indicates the anisotropic nature of subcritical turbulence. In future work, the properties of subcritical turbulence including details of the non-linear sustaining process, as well as the precise role of the cross-rolls in it, should be studied employing, in particular, 3D Fourier analysis in k-space (see Mamatsashvili et al. (2016) and Gogichaishvili et al. (2017) for an examples of such studies in application to non-rotating shear flow and Keplerian flow with an azimuthal magnetic field).

Shearing spirals considered analytically in the leading order in small k_z allow us to draw the following conclusions.

• i)

  The dynamics of perturbations in the plane orthogonal to the spanwise direction remain approximately the same as for spanwise invariant (columnar) shearing spirals, which is the 2D swing amplification process.

• ii)

  Spanwise acceleration of fluid elements is determined by the swing amplification.

• iii)

  Spanwise velocity perturbation increases while the shearing spiral swings from trailing to leading in the cyclonic flow and vice versa in the quasi-Keplerian flow. After that, the spiral is wound by the shear, which causes the decay of streamwise and shearwise velocity perturbations, whereas the spanwise velocity perturbation tends to an asymptotic inviscid non-zero value proportional to the net impulse given to fluid elements in spanwise direction. The latter depends on the shear rate and vanishes as q → 2. Viscous damping results in a maximum value of the initial winding of shearing spiral and additionally in the maximum of spanwise velocity, which appears after the instant of swing (see Eqs. (36) and (37) giving the locations of those extrema in the limit of tightly wound spirals).

• iv)

  Assuming that the cross-rolls of the largest amplitude sustain turbulence throughout the range of centrifugally stable rotating shear flows and the threshold value of their growth factor found in the cyclonic regime (see Eq. (69)) remains the same, we obtain a tentative prediction of R_T(q) for the quasi-Keplerian regime, where hydrodynamical turbulence is not yet discovered. In particular, we find that the lowest R_T ∼ 10⁷ in the quasi-Keplerian regime must be attained at a sub-Keplerian shear rate q ≈ 1/2.

• v)

  The threshold value of the cross-rolls growth factor (69) exceeds the maximum growth factors of streaks and rolls generated via lift-up and anti-lift-up up to the shear rates, very close, respectively, to non-rotating flow and to flow on the Rayleigh line (see Sect. 5.1).

The cross-rolls’ appearance is shown in the top panels in Fig. 9 (see also Appendix A for the explanation of this representation). The growth factor of the corresponding shearing spiral is plotted by a solid curve in Fig. 1. For comparison, the rolls generated from axisymmetric streaks via anti-lift-up with the same k_x and k_z but in the flow close to the Rayleigh line are shown on the bottom panels in Fig. 9. The former and the latter ones are aligned in the shearwise and the streamwise directions, respectively. Furthermore, in contrast to rolls, the cross-rolls are essentially 3D vortices, since they are of a finite length along their axes of rotation, being contracted by the background shear ∝ |k̃_x|⁻¹. In combination with the opposite rotation of the adjacent cross-rolls, this causes monotonic growth of ω_y ∝ |k̃_x| in areas between the adjacent cross-rolls, in the vicinity of for example point (x = λ_x/4, y = 0, z = λ_z/4). As one can see, the cross-rolls coexist with “planar” eddies fading after the instant of swing. Such a complicated structure of the perturbation velocity field requires a new concept of non-linear feedback, which may be completely different from both SSP and a scheme associated with streamwise rolls in the flow on the Rayleigh line. The latter itself remains unknown.

The prediction of R_T(q) obtained in this work (see Fig. 8) makes it clear that, at least for boundless homogeneous flows, the issue of the transition to turbulence in the Keplerian regime should be attacked from the side of the cyclonic regime. This implies the crossing of the solid-body line, R_Ω = +∞ →R_Ω = −∞, which is not evident while looking at the axis of R_Ω. However, we suppose that R_T ∼ 10⁷ at the most “favourable” sub-Keplerian shear rate q ∼ 1/2, which presumably requires a resolution of N ≳ 2000 for the spectral Fourier code to be able to test the non-linear instability of the quasi-Keplerian regime. Still, such an extraordinary non-linear stability of Keplerian flow may be defeated by the huge spatial scales of real astrophysical discs, leading to R = 10¹⁰ ÷ 10¹³.

---

Acknowledgments

We kindly thank Professor K. Postnov and Professor P. Ivanov for genuine interest in the subject of this study and for a number of discussions during the preparation of the manuscript. The research is carried out using the equipment of the shared research facilities of HPC computing resources at Lomonosov Moscow State University (see paper Sadovnichy et al. 2013 for its detailed description). Equipment for the reported study was granted in part by the M. V. Lomonosov Moscow State University Programme of Development. DNR was supported by grant RSF 14-12-00146 when writing Sect. 6 of this paper. VVZ was supported in part by RFBR grant 15-02-08476 and by programme 7 of the Presidium of Russian Academy of Sciences.

References

Appendix A: Representation of cross-rolls

In order to plot the velocity perturbation in Fig. 9, we take the imaginary part of SFH given by Eq. (5), where the Fourier amplitudes of the velocity perturbations, û_x,y,z, are given by the analytic solution (17,18, 25). Explicitly, we have

$$\begin{array}{c}\hfill \mathfrak{I}\left[\begin{array}{c}{u}_x\\ u_y\\ u_z\end{array}\right]=\left(\begin{array}{c}{\widehat{u}}_x\\ {\widehat{u}}_y\\ {\widehat{u}}_z\end{array}\right)sin({\tilde{k}}_{x}x+y+k_{z}z)={\mathbf{C}}_1+{\mathbf{C}}_2+{\mathbf{C}}_3+{\mathbf{C}}_4,\end{array}$$(A.1.)

where

$$\begin{array}{c}\hfill {\mathbf{C}}_2({\tilde{k}}_{x}x,y,k_{z}z)={\mathbf{C}}_1({\tilde{k}}_{x}x+\pi /2,y-\pi /2,k_{z}z),\\ \hfill {\mathbf{C}}_3({\tilde{k}}_{x}x,y,k_{z}z)={\mathbf{C}}_1({\tilde{k}}_{x}x+\pi /2,y,k_{z}z-\pi /2),\\ \hfill {\mathbf{C}}_4({\tilde{k}}_{x}x,y,k_{z}z)={\mathbf{C}}_1({\tilde{k}}_{x}x,y+\pi /2,k_{z}z-\pi /2)\end{array}$$

Therefore, SFH is nothing but a combination of the single velocity pattern, C₁, with its duplicates translated along the coordinate axes on one-fourth of the corresponding wavelength. This velocity pattern has the form

$$\begin{array}{c}\hfill {\mathbf{C}}_1({\tilde{k}}_{x}x,y,k_{z}z)=\left(\begin{array}{c}{\widehat{u}}_{x}sin({\tilde{k}}_{x}x)cos(y)cos(k_{z}z)\\ {\widehat{u}}_{y}cos({\tilde{k}}_{x}x)sin(y)cos(k_{z}z)\\ {\widehat{u}}_{z}cos({\tilde{k}}_{x}x)cos(y)sin(k_{z}z)\end{array}\right)\end{array}$$(A.2.)

and is represented in Fig. 9 by the three mutually perpendicular cross-sections having the common point {x = 0, y = 0, z = 0}.

Appendix B: Effect of finite numerical resolution

In order to reveal the influence of the numerical viscosity on the transition to turbulence at terminally high R_nl, we fix q = −4 and perform a set of additional simulations with exactly the same setup as described in Sect. 6. The results are shown in table B.1.

The value of R_T obtained in cubic box substantially drops as we increase resolution up to N = 192. This is to be expected, since at N = 128 R_T is close to the restriction (68). However, an even higher resolution, N = 256, does not lead to further decrease of R_T: we conclude that convergence is achieved in this case. At the same time, R_T in tall box is only weakly affected by the increase of resolution. Thus, there is a difference between the transition Reynolds numbers in cubic and tall boxes, which cannot be attributed to effect of finite resolution.

All Tables

All Figures