Excitation of turbulence in accretion disks of binary stars by non-linear perturbations

TL;DR

This paper studies how nonlinear perturbations in accretion disks of binary stars can generate turbulence through shock formation and flow instability, with turbulence characteristics influenced by perturbation parameters.

Contribution

It investigates the nonlinear evolution of finite-amplitude perturbations and their role in turbulence excitation in accretion disks, highlighting shock formation and flow instability mechanisms.

Findings

Shock waves form at perturbation fronts causing density and velocity discontinuities.

Turbulence characterized by alpha-parameter up to ~0.1.

Flow instability leads to disk turbulization influenced by perturbation parameters.

Abstract

Accretion disks in binary systems can experience hydrodynamic impact at inner as well as outer edges. The first case is typical for protoplanetary disks around young T Tau stars. The second one is typical for circumstellar disks in close binaries. As a result of such an impact, perturbations with different scales and amplitudes are excited in the disk. We investigated the nonlinear evolution of perturbations of a finite, but small amplitude, at the background of sub-Keplerian flow. Nonlinear effects at the front of perturbations lead to the formation of a shock wave, namely the discontinuity of the density and radial velocity. At this, the tangential flow in the neighborhood of the shock becomes equivalent to the flow in in the boundary layer. Instability of the tangential flow further leads to turbulization of the disk. Characteristics of the turbulence depend on perturbation parameters, but alpha-parameter of Shakura-Sunyaev does not exceed ~0.1.

Equations

$$\frac{1}{{\rho }_0}\frac{\partial p_0}{\partial r}=-r{\Omega }_{\mathrm{K}}^2+\frac{{\lambda }_0^2}{r},$$

$${\Omega }_{\mathrm{K}}={\Omega }_{\mathrm{K}\ast }{\left(\frac{r}{r_{\ast }}\right)}^{-3/2}.$$

$$r{\Omega }_{\mathrm{K}}^2-\frac{{\lambda }_0^2}{r^3}=\chi \frac{H^2}{r^2}r{\Omega }_{\mathrm{K}}^2,$$

$${\rho }_0(r)={\rho }_{\ast }{\left(\frac{r}{r_{\ast }}\right)}^{-d},$$

$$p_0=p_{\ast }{\left(\frac{{\rho }_0}{{\rho }_{\ast }}\right)}^{\gamma }=p_{\ast }{\left(\frac{r}{r_{\ast }}\right)}^{-\gamma d},$$

$$c_0^2(r)=c_{\ast }^2{\left(\frac{r}{r_{\ast }}\right)}^{-(\gamma -1)d},$$

$$d{c}_{\ast }^2{\left(\frac{r}{r_{\ast }}\right)}^{-(\gamma -1)d+3}=\chi {\Omega }_{\mathrm{K}\ast }^2{H}^2(r).$$

$$\rho (r,z)={\rho }_0(r)\exp \left[-\frac{z^2}{2H^2(r)}\right],$$

$$H^2(r)=\frac{c_0^2(r)}{{\Omega }_{\mathrm{K}}^2(r)}=\frac{c_{\ast }^2}{{\Omega }_{\mathrm{K}\ast }^2}{\left(\frac{r}{r_{\ast }}\right)}^{-(\gamma -1)d+3}.$$

$$\chi =d.$$

$$\Sigma (r)=2{\int }_0^{\infty }dz\rho (r,z)=\sqrt{2\pi }\frac{{\rho }_{\ast }{c}_{\ast }}{{\Omega }_{\mathrm{K}\ast }}{\left(\frac{r}{r_{\ast }}\right)}^{[-(\gamma +1)d+3]/2}.$$

$$\gamma =5/3,d=3/2,(\gamma -1)d=1.$$

$${\rho }_0={\rho }_{\ast }{\left(\frac{r}{r_{\ast }}\right)}^{-3/2},$$

$$c_0=c_{\ast }{\left(\frac{r}{r_{\ast }}\right)}^{-1/2},$$

$${\lambda }_0={\left(1-\frac{3}{2}\frac{c_{\ast }^2}{r_{\ast }^2{\Omega }_{\mathrm{K}\ast }^2}\right)}^{1/2}{\left(\frac{r}{r_{\ast }}\right)}^{1/2}{r}_{\ast }^2{\Omega }_{\mathrm{K}\ast },$$

$$H=\frac{c_{\ast }}{{\Omega }_{\mathrm{K}\ast }}\frac{r}{r_{\ast }}.$$

$$\frac{\partial \rho }{\partial t}+\frac{1}{r}\frac{\partial (r\rho v)}{\partial r}=0,$$

$$\frac{\partial v}{\partial t}+v\frac{\partial v}{\partial r}=-\frac{1}{\rho }\frac{\partial p}{\partial r}-r{\Omega }_{\mathrm{K}}^2+\frac{{\lambda }^2}{r^3},$$

$$\frac{\partial \lambda }{\partial t}+v\frac{\partial \lambda }{\partial r}=0,$$

$$\lambda^{2}(t,r)=\lambda_{0}^{2}\bigl{(}r-\xi(t,r)\bigr{)}=\lambda_{0}^{2}(r)-\frac{d\lambda_{0}^{2}(r)}{dr}\,\xi(t,r)+\mathcal{O}\!\left(r\Omega_{\mathrm{K}}^{2}\frac{\xi^{2}}{r^{2}}\right)\;.$$

$$\frac{d{\lambda }_0^2}{dr}=\left(1-\frac{3}{2}\frac{c_{\ast }^2}{r_{\ast }^2{\Omega }_{\mathrm{K}\ast }^2}\right)r_{\ast }^3{\Omega }_{\mathrm{K}\ast }^2.$$

$$\frac{\partial {\rho }_1}{\partial t}+\frac{\partial ({\rho }_0{v}_1)}{\partial r}=0.$$

$$\frac{\partial v_1}{\partial t}+\frac{1}{{\rho }_0}\frac{\partial p_1}{\partial r}-\frac{d{p}_0}{dr}\frac{{\rho }_1}{{\rho }_0^2}+\frac{1}{r^3}\frac{d{\lambda }_0^2}{dr}\xi =0.$$

$$\frac{d{p}_0}{dr}=c_0^2\frac{d{\rho }_0}{dr},$$

$$p_1=c_0^2{\rho }_1,$$

$$\frac{d{c}_0^2}{dr}=(\gamma -1)\frac{c_0^2}{{\rho }_0}\frac{d{\rho }_0}{dr}.$$

$${\rho }_1=-\frac{\partial ({\rho }_0\xi )}{\partial r},$$

$$v_1=\frac{\partial \xi }{\partial t}.$$

$$\frac{{\partial }^2\xi }{\partial t^2}=\frac{c_0^2}{{\rho }_0}\frac{{\partial }^2({\rho }_0\xi )}{\partial r^2}+(\gamma -2)\frac{c_0^2}{{\rho }_0^2}\frac{\partial {\rho }_0}{\partial r}\frac{\partial ({\rho }_0\xi )}{\partial r}-\frac{1}{r^3}\frac{d{\lambda }_0^2}{dr}\xi .$$

$$\xi =r_{\ast }\Re \left\{C{e}^{\mp i\omega t}{\left(\frac{r}{r_{\ast }}\right)}^{7/4}{H}_{\nu }\left[\frac{2r_{\ast }\omega }{3c_{\ast }}{\left(\frac{r}{r_{\ast }}\right)}^{3/2}\right]\right\},$$

$$\nu =\frac{1}{6}{\left(16\frac{r_{\ast }^2{\Omega }_{\mathrm{K}\ast }^2}{c_{\ast }^2}-23\right)}^{1/2}.$$

$${\omega }^2>\left(1-\frac{3}{2}\frac{c_{\ast }^2}{r_{\ast }^2{\Omega }_{\mathrm{K}\ast }^2}\right){\Omega }_{\mathrm{K}}^2\approx {\Omega }_{\mathrm{K}}^2.$$

$$\frac{r}{r_{\ast }}\gtrsim {\left(\frac{{\Omega }_{\mathrm{K}\ast }}{\omega }\right)}^{2/3}.$$

$$\xi \approx |C|r_{\ast }{\left(\frac{3c_{\ast }}{\pi r_{\ast }\omega }\right)}^{1/2}\frac{r}{r_{\ast }}\cos \left[\frac{2r_{\ast }\omega }{3c_{\ast }}{\left(\frac{r}{r_{\ast }}\right)}^{3/2}\mp \omega t+(\dots )\right],$$

$${\rho }_1\approx |C|{\rho }_{\ast }{\left(\frac{3r_{\ast }\omega }{\pi c_{\ast }}\right)}^{1/2}\sin \left[\frac{2r_{\ast }\omega }{3c_{\ast }}{\left(\frac{r}{r_{\ast }}\right)}^{3/2}\mp \omega t+(\dots )\right],$$

Full text

PACS numbers:  47.27.-i, 97.10.Gz, 97.80.-d

1 Introduction

Turbulent viscosity is the most effective mechanism of the angular momentum transfer in astrophysical disks [1, 2, 3]. The turbulence is caused by development of instability of small perturbations. In accretion disks of nonmagnetic stars, as well as in the protoplanetary disks (in cold non-ionized gas), the turbulence can develop from instabilities of the hydrodynamic type only.

The stability of linear hydrodynamic perturbations in Keplerian disks is studied in numerous papers (see, for instance, references in [4]). Both radial and non-axisymmetric perturbations of the modal type turn out to be stable according to the Rayleigh criterion [5, 6]. The reason for this is spectral stability of the Keplerian flow. There are non-modal type perturbations or transient ones, which display a growth in the linear stage (see a review in [7]). However, the problem of the mechanism of their continuous generation in the Keplerian flow, as well as the problem of efficiency of removal of angular momentum, on large time scales needs further study.

There are at least two classes of astrophysical disks in which accretion flow experiences an external action: the disks around components of close binary stars and disks around young stellar systems (protoplanetary ones). In the objects of the first kind, the stream from the donor-star acts upon the outer part of the accretion disk, generating a tangential discontinuity and associated with the latter shock wave, the so-called „hot line“ [8, 9]. In the objects of the second kind act detached shocks from the components of the binary system [10] or shockless perturbations, if the star moves with a subsonic velocity [11]. The last affect the protoplanetary disk at the inner side. As a result, nonlinear waves of different amplitude and scale may be exited in the accretion disk. Although in a certain conditions the nonlinear waves may provide a conducive background for linear instability and the turbulence formation [12], it is not the case in the protoplanetary disks, since the nonlinear perturbations are too weak.

In hydrodynamics, the mechanism of formation of discontinuities in nonlinear waves is known [13, 6]. This phenomenon is due to the “breaking” of the wave profile. If the discontinuity is a shock wave, dissipative effects will lead to the damping of the perturbation with time. On the other hand, in the accretion flows the tangential velocity of a gas can many times larger than the sound velocity. Therefore, it can not be ruled out that, for certain perturbation parameters the rollover of the profile may cause a tangential discontinuity. In this case, the difference of energy at different sides of the discontinuity will transform into kinetic energy of the turbulence.

In the present paper, we investigate whether evolution of nonlinear waves in the accretion disk may lead to the appearance of unstable configurations and, as a result, to turbulence. we We perform an analysis of nonlinear perturbations in sub-Keplerian disks with inhomogeneous distribution of density and pressure. The paper is organized as follows: in the Section 2 we describe the structure of accretion flow in semi-detached and young binary systems. In Section 3 we present the model of an equilibrium sub-Keplerian disk and a briefly analyze linear perturbations in the latter. In the fourth Section we explore evolution of nonlinear radial perturbations in accretion disks. In Section 5 we apply our results to two simple models of accretion disks. Conclusions and deductions are presented in the last Section.

2 The structure of accretion flows in young binaries with circumstellar envelope and in close binaries

Accretion disks in binary systems have the distribution of angular momentum close to the Keplerian one. Velocity of tangential motion in them can be hundreds times larger than the sound velocity (in the central regions). Often, accretion disks experience hydrodynamic action from accretion flows, encounter with the streams of incoming gas or detached shock waves. Such effects, along to large tangential velocities, can be potential sources of turbulence. Below, we will consider the flow patterns for two classes of typical binary stars: young stars surrounded by a protoplanetary disk, and close binary systems.

The envelope of a T Tauri binary star is a remnant of the protostellar cloud from which the matter falls and forms a Keplerian protoplanetary disk. In the central part of the disk around the binary forms a gas-free region, a “gap”. The reason for the existence of the gap can be both tidal action (Lindblad resonances) of the binary [14] and the effects of the detached shock waves produced by the components of the system [10]. Typical gap width is of the order of several orbital separations (see, e.g., [15, 16]). Inside the gap, around each of the stars forms its own circumstellar accretion disk and in the case of the supersonic orbital motion of components (which is typical of such systems), detached shock waves appear in the front of the stars [17, 18, 19]. In the region of the inner disk boundary, cold and dense gas of the Keplerian flow is contiguous to the hot and rarefied gas of the gap. Keplerian velocity in this region can exceed the sound velocity by many dozen times. At the same time, the impact of the detached shock wave upon the disk is concentrated on a small part of the boundary. Since the impact has a shock nature, and it is quite intense, it is obvious that it will result in destruction of the Keplerian flow in the vicinity of the contact. Blurring of this region by differential rotation of the gas can lead to the formation of a turbulent layer over a time of the order of one revolution of the disk at the radius of the inner boundary.

In close binary systems, one of the components (the donor) fills its Roche lobe. This results in the flow of the matter onto the second component (accretor). Under the influence of the dissipative processes, a disk forms around the accretor. The former has a Keplerian angular momentum distribution. Typical is the case when the accretor is more massive and is a $\sim 1\leavevmode\nobreak\ M_{\odot}$ white dwarf, while the mass of the donor is lower. The radius of the accretion disk comprises the tenth of the orbital separation, while velocity of Keplerian motion even at the outer edge of the disk can exceed the sound speed be several tens of times. The region of the outer edge of the disk where the stream from the donor interacts with accretion disk, is a tangential discontinuity [8, 9, 20]. The step of tangential velocity at the discontinuity can be of the order of the Keplerian velocity at the edge of the disk. It is known that such discontinuities are unstable for many types of perturbations [21]. Because of the high flow velocities, the instability must have a drift character and manifest itself as a turbulent wake from the tangential discontinuity in the direction of rotation of the gas.

In both examples considered above, the turbulence is excited in the Keplerian flow: in the inner region of the protoplanetary disks and in the outer region of the circumstellar disk in semi-detached binary systems. The turbulence, however, engulfs only a small part of the disk flow. The problem of possible further propagation of turbulence into entire volume of the accretion disk is of particular interest, since if this is a case, existing nonlinear perturbations can be considered as the source of turbulence of nonmagnetic accretion disks.

The effect of detached shock waves upon internal boundary of the protoplanetary disk should not only lead to the destruction of the boundary Keplerian flow, but also to the propagation of the perturbations with different amplitudes and scales inside the disk. A similar situation should take place also in the disks of close binaries — slight perturbations of the rate of mass loss by the donor star or characteristics of the flow at the edge of the accretion disk (for example, due to the intrinsic precession of the accretion disk) can cause changes in the position and intensity of the tangential discontinuity, which, in turn, will disturb the accretion flow.

In the Introduction, we mentioned that small amplitude perturbations are either stable (for instant, radial ones are stable by the Rayleigh criterion [22]), or display only a weak instability (some nonaxisymmetric perturbations are unstable [7], but they have rather large time scale of growth). On the contrary, non-linear waves can participate in the universal scenario of generation of shock waves and tangential discontinuities which is realized in any environment with ordinary equation of state [13, 6]. At the conditions of fast shear flow these discontinuities will be unstable and will decay directly, generating turbulence.

A nonlinear wave is a propagating perturbation of density, pressure, and velocity, with a finite amplitude. Two factors are responsible for the formation of a discontinuity in a nonlinear wave: (i) in the areas of enhanced pressure, propagation velocity of perturbations is higher than in the low-pressure regions; (ii) velocity of perturbations is added with the velocity of the background flow. So, if consider the leading edge of a nonlinear wave, where the pressure and velocity vary from the maximum to the minimum values, the region of the maximum will catch up with the minimum region. This will lead to the steepening of the wave profile and to its “rollover” and formation of a discontinuous solution in a finite time. This pattern becomes somewhat more complicated if one takes into account the differential rotation of the disk matter. For instance, the wave solution for small radial perturbations exists only in the region of the Keplerian disk, where the frequency of the wave is not less than the angular velocity of the matter of the disk [23]. However, we will show below that the scenario of formation of discontinuities can take place in the rotating flows too.

Next, we will perform a brief analysis of the linear perturbations in the accretion disks and will consider the mechanism of formation of discontinuities in the nonlinear perturbations in such disks.

3 Linear perturbations in the accretion disk

3.1 Equilibrium disk

Let consider a stationary disk flow, symmetric respective to the $r=0$ axis and $z=0$ plane. The equation of radial balance of forces in the mid-plane of the disk:

[TABLE]

where $\rho_{0}$, $p_{0}$, $\lambda_{0}$ are, respectively, density, pressure, and angular momentum in the $z=0$ plane; $\Omega_{\mathrm{K}}$ is Keplerian angular velocity

[TABLE]

By „$\ast$“ we will mark the values of variables in the point $r_{\ast}$ at $z=0$. In particular, $\Omega_{\mathrm{K}\ast}=(GM/r_{\ast}^{3})^{1/2}$. The gas pressure usually decreases with the radius. Therefore the centrifugal acceleration is not exactly equal to the gravitational one, but is somewhat lower [24]. Let the disk to be sub-Keplerian, if the balance of gravitational and centrifugal force is violated by a small amount [25]:

[TABLE]

where $H$ is semi-thickness of the disk; $\chi$ — a dimansionless parameter. Let distribution of the matter to be isentropic with the exponent $\gamma$, distribution of density $\rho_{0}$ on radius is a power law with exponent $d$. Then, for the pressure, density and sound velocity one may write down

[TABLE]

where $c_{\ast}^{2}=\gamma p_{\ast}/\rho_{\ast}$. Then the equation of the radial balance of forces becomes

[TABLE]

Let assume that the vertical structure of the disk is isothermal. In a thin disk, $H\ll r$, the density is distributed as

[TABLE]

where semi-thickness of the disk is determined from the condition of the vertical balance of forces

[TABLE]

From this expression, using Eq. (7) we may derive the factor $\chi$ in the Eq. (3):

[TABLE]

For the accepted distribution $\rho(r,z)$ the dependence of the surface density on the radius will be a power law:

[TABLE]

In the realistic models of the accretion disks, surface density declines with radius, approximately, as $1/r^{0\dots 1}$ (see [25] and references therein). For instance, in the Shakura-Sunyaev disk model $\Sigma\propto r^{-3/5}$ [2]. In the paper [26] as a typical dependence is assumed $\sigma\propto r^{-1}$. In our model, we shall assume $\Sigma\propto r^{-1/2}$. Below, this model will allow to obtain more simple expressions for some of dependences, while remaining qualitatively consistent with more complete models.

The slope of the surface density distribution in combination with the value of the adiabatic exponent $\gamma$, determines the values of exponents in power-laws for remaining thermodynamical quantities: volumetric density, pressure, and sound velocity. Further in the present study, we will consider the gas as monatomic. Thus, we obtain

[TABLE]

Let write down final expressions for equilibrium distributions:

[TABLE]

The typical values of the gas velocity in the disks, according to numerical models, are equal to tens of sound velocities [8]. Thus, correction to the Keplerian momentum in expression (15) does not exceed one per cent.

3.2 Linear perturbations

Let consider the dynamics of gas flow in the mid-plane of the disk and confine ourselves to perturbations of the equilibrium configurations that (i) are symmetric about the axis $r=0$; (ii) do not depend on the vertical coordinate z; (iii) have a zero velocity in the direction of the z-axis. The system of equations that describes such a class of flows has the form:

[TABLE]

where $v$ is radial velocity; $\lambda$ — angular momentum per unit mass.

Let assume that the system of equations (17)–(19) describes a small perturbation of the equilibrium configuration, which was determined in the previous section. The perturbation itself we describe as a radial displacement of gas elements. Let denote by $\xi(t,r)$ the value of radial displacement of the particle whose position at the time $t$ is $r$, i.e., the radial coordinate of the nonmoved particle is $r-\xi(t,r)$. From the conservation of the angular momentum law in form (19) it follows that the moment of each element of the gas is conserved, i. e., we may write down

[TABLE]

Here we took into account that $|\xi|\ll r$. In our model the derivative $d\lambda_{0}^{2}/dr$ does not depend on the coordinate and is equal to

[TABLE]

Let require that perturbations of interest have a small spatial scale of variability. Then we can neglect the „geometric“ term $\rho v/r$ in the equation of continuity, which is determined by the cylindrical geometry of the problem, and the advective term in (17) will take a Cartesian form. This approximation was used by many authors in the analysis of perturbations in the accretion disks (see, for instance, [27, 22]).

Let denote perturbed values by an index „$1$“. Then, after linearization of equation of continuity (17) over small perturbations, we get

[TABLE]

By means of (20) linearization of the Euler equation (18) will be

[TABLE]

Angular momentum conservation law was already taken into account explicitly when we employed expansion (20). The gradients of thermodynamical values for adiabatic equation of state have the form:

[TABLE]

To close the system (22), (23), let use the relation between dislocation, density and velocity [3]:

[TABLE]

It easy to show that, if these definitions are used, continuity equation (22) is satisfied automatically. Inserting (27) and (28) into the radial component of equations of motion equation (23) we get a closed equation for the displacement:

[TABLE]

The algorithm of solution of this equation is presented in the Appendix A. Solution has the form of a progressing wave:

[TABLE]

where $H_{\nu}$ is the Hankel function of the first kind of the order $\nu$; $C$ is a complex constant; $\Re\{\cdot\}$ — real part of its argument. The order is defined as

[TABLE]

It is necessary to note that the solution of Eq. (30) has a wave form only in the case, when the following inequality keeps (see Appendix A):

[TABLE]

A similar result was obtained in [23] by WKB-method. The last inequality may be rewritten as

[TABLE]

Thus, perturbation wave with a given frequency $\omega$ may progress in the outer part of the disk only, where $\omega>\Omega_{\mathrm{K}}$.

At large distance from the center dependence of radial displacement on the radius gets a simpler form (see Appendix A):

[TABLE]

where $(\dots)$ is an insignificant phase addition. Using (27) and (28) we can obtain expressions for perturbations of density and velocity:

[TABLE]

Let get an expression for perturbation of the angular momentum. Using expansion (20) one may write down

[TABLE]

Applying approximate relation (34), we get

[TABLE]

It is easy to estimate the typical space scale of the wave $\ell$ from oscillating part of solution (34). As the definition of $\ell$ let write down the equation for the phase:

[TABLE]

where $[\dots]$ stands for the term in square brackets in (34). Then

[TABLE]

Comparing this expression with condition (32) one may note that the perturbations have a wave form, if their „wavelength“ $\ell$ does not exceed semi-thickness of the disk. In other words, the waves can not propagate into the region $\ell>H$.

Linear perturbations in the rotating flows were considered by many authors. For instance, in the classical paper [23] by WKB-method, was obtained a dispersion relation for the waves:

[TABLE]

where $k$ is the parameter of the WKB-method, which is interpreted as the wave number. Such an interpretation implies that the sense of the value $k^{-1}$ is the length of the perturbation wave. This value is related to $\ell$ from (40) as

[TABLE]

In the depth of the wave zone ($\omega\gg\Omega_{\mathrm{K}}$) the values of both expressions comply with precision up to value of the order $\Omega_{\mathrm{K}}^{2}/\omega^{2}$. However, in the region $\omega\approx\Omega_{\mathrm{K}}$ definition of wave length by expression (42) provides a physically unacceptable result, since $k^{-1}\to\infty$. Below, as the length of perturbation wave we will bear in mind expression (40), which, in essence, is an estimate of the distance between neighboring maxima of the wave.

Description of the dynamics of small perturbations in the terms of radial displacement $\xi$ is completely equivalent to the usual approach to the analysis of linear waves in the accretion disks, examples of which can be found in [27, 28] and many other papers. This method was used repeatedly earlier, see, e.g., description of the general approach to the analysis of Lagrangian perturbations in [3]. In the present study, the use of the Lagrangian approach is preferable, since it allows to express more easily the initial conditions for the analysis of nonlinear perturbations, which is carried out in the next Section.

4 Formation of discontinuities

It is known that nonlinear perturbations propagating in the gas can evolve with time into shock waves [6, 13]. Usually, the proof of this statement is based on the solution in the form of a simple Riemann wave for one-dimensional gas dynamics [6]. In the polytropic flow, described by a simple wave, all gasdynamic quantities are expressed through one variable. As it can be seen from the expressions (35) and (36), in the case of a perturbation of inhomogeneous medium, density and velocity distributions depend on the coordinate differently and, therefore, they do not form a simple wave. However, it is possible to show that in this case the perturbation can also evolve into a shock wave. If there is a nonzero tangential velocity, on both sides of the shock wave forms a flow, equivalent to the flow in the boundary layer. To prove this assertion, we use the formulation of equations for one-dimensional gas flow presented in the book [6].

First, we note that we are interested in perturbations from the wave part of the accretion disk, with wavelength which does not exceed the semi-thickness of the disk (see the previous Section). Let us consider, for definitiveness, a perturbation in the form of a wave (35), (36). If nonlinear effects would be taken into account, the form of this perturbation will be distorted, but such characteristics as the amplitude and typical scale of variability („wavelength“ $\ell$) will be about the same as in the linear case, at least up to the beginning of the formation of discontinuities.

Let estimate the contribution of various terms to the linearized Euler equation (23). The amplitude of velocity perturbations we normalize by the Mach number $\mathcal{M}$: $|v_{1}|\sim\mathcal{M}c_{0}$. Then the amplitude of density perturbations can be estimated as $|\rho_{1}|\sim\mathcal{M}\rho_{0}$. During one wave period, every gas element oscillates along the direction of propagation of the perturbation with an amplitude of the order of $|v_{1}|\ell/c$. This value can be used as an estimate of the displacement amplitude $|\xi|$. It follows from these estimates that in the progressing wave, the ratio of pressure force and the sum of the gravitational and centrifugal forces is

[TABLE]

In the wave part of the disk $\ell<H$, therefore, in order to avoid complicating of further presentation, we neglect the influence of gravitational and centrifugal acceleration upon dynamics of the wave.

It may seem that, at the background of our approximations, all effects caused by differential rotation disappear. Indeed, we demanded that perturbations have sufficiently small scale of variability in comparison with the typical disk radius and the perturbation amplitude is small enough that it would be possible to interpret the sum of the gravitational and centrifugal accelerations in terms of the displacement of the gas element $\xi$. In additive, as it can be seen from the solution of the linearized problem (subsection 3.2), for sufficiently large $r$ the differential rotation law does not determine directly the dependence $\xi(t,r)$, but enters only the phase of the oscillating term. In reality, the role of differential rotation will be manifested in the fact that the distribution of the unperturbed density and the speed of sound are inhomogeneous in space and the law of this distribution is dictated, among other factors, by the law of differential rotation (see 3.1).

We reduced the gas flow equations to the one-dimensional Euler equation and continuity equation. It is shown in [6] that an arbitrary one-dimensional isentropic gas flow with the adiabatic index $\gamma$ may be described by a linear equation of the form

[TABLE]

where $v$ is the velocity; $w=\int dp/\rho=c^{2}/(\gamma-1)$ — enthalpy. Solution of this equation, the function $\chi(v,w)$, is defined as

[TABLE]

where $\varphi$ is velocity potential, i.e. $v=\partial\varphi/\partial r$. Time $t$ and coordinate $r$ are expressed as:

[TABLE]

We are not interested in the complete solution of the problem of nonlinear evolution of a perturbations with a given initial form, but only in the analysis of the possibility of formation of discontinuities in it. The location of the discontinuity can be defined as the point at which the plots of the velocity and/or enthalpy dependence on the coordinate become multivalued. In terms of the function $r(v,w)$ the position of discontinuity corresponds to the inflection point:

[TABLE]

while the derivatives are computed at the fixed instant of the time. It is easy to see that the r.h.s. of (47) describes dependence of $r(v,w)$ at $t=0$; we will denote it as $R(v,w)$. The first of expressions (48) together with (47) provides the time of formation of discontinuity in solution:

[TABLE]

The calculation of this derivative in the general case may turn out to be a non-trivial task. However, we can express the position of the element we can express position of a gas element $R$ кas a function of velocity, by inversion of the dependence of the velocity on coordinates at $t=0$. Let take as an initial perturbation configuration a section of the profile of a linear wave (36):

[TABLE]

where $\mathcal{M}_{\ast}$ is Mach number at $r_{\ast}$. Without limiting generality, let select $r_{\ast}$ as a point, where velocity becomes zero and declines in the direction of progression of the wave. It is evident that this point is also inflection point for dependence $r(v)$. Let consider the wave in a small vicinity of this point. Let introduce $x=r-r_{\ast}$; it will not be an additiveal restriction to require satisfaction of condition $|x|\ll r_{\ast}$:

[TABLE]

Now, with precision down to small quantities of the second order, we may write down:

[TABLE]

Inserting this expression into (47) and applying the first of conditions (48) (the second one is fulfilled automatically), we obtain an expression for the time of discontinuity formation:

[TABLE]

where $\mathcal{M}_{\ast}$ and $\omega$ should be interpreted as parameters of perturbation in the point $r_{\ast}$.

It is easy to find the point of formation of discontinuity $r_{\mathrm{sh}}$, if we recall that the local phase velocity of small-scale perturbations is equal to the sound velocity. In the general case of an arbitrary direction of progression of the wave we obtain

[TABLE]

The sign in the front of r.h.s. depends on direction of the wave propagation — the outer edge of the disk („$+$“) or to the inner one („$-$“). Inserting (53) and expression for the sound velocity $c_{0}=c_{\ast}\,(r/r_{\ast})^{-1/2}$ into (54), we obtain

[TABLE]

Note, for the waves progressing to the center of the accretion disk, the formula (55) is valid in the wave region, i.e., as long as condition $\omega>\Omega_{\mathrm{K}}(r_{\mathrm{sh}})$ is valid (see subsection 3.2). The time it takes for perturbation to reach the border of the wave zone may be estimated from the latter condition and expression (54):

[TABLE]

If the discontinuity did not form during this time, then, following energy conservation law, the wave should reflect from the wave region boundary and head to the outer part of the disk. In this case, the wave has one more chance to form the discontinuity. However, the very reflection of the wave from the boundary $\ell=H$ can not be described in terms of the present model.

It is evident that all above mentioned is valid also for dependence of perturbation of enthalpy on radius. At this, from the form of functions (35) and (36) it is clear that inflection points for $r(v)$ and $r(w)$ coincide.

It is important to note that the angular momentum distribution in the perturb flow will behave differently. Indeed, expressions (35), (36) and (38) show that, if the dependence of density and radial velocity passes through zero (the front of the perturbation wave), the angular momentum passes through the extremum. This means that in the plot of the dependence of the angular momentum on the radius does not appear discontinuity. However, since the angular momentum is transferred by the gas as a passive scalar, see (19), at the instant of formation of the discontinuity of density and radial velocity, the distribution of angular momentum will get the shape of a caustic with a peak located at the point $r_{\mathrm{sh}}$. Schematically this is shown in Fig. 1. As it is seen, the amplitude of the density and radial velocity step is determined by the change of the corresponding value over the scale of half wavelength, $\ell/2$, of the initial perturbation. The change of tangential velocity on the caustic is determined mainly by the difference of the values of the tangential velocity in the original distribution over the scale $\ell/4$.

5 Application to accretion disks

In the previous Section we obtained the estimates of the formation time of discontinuities in a rotating flow. As was shown in [6], discontinuities of density and velocity, formed as a result of rollover of the wave profile are asymptotically damped because of microscopic viscosity. However, in the case of a rotating flow, the picture of the evolution of a perturbation after formation of a discontinuity may change qualitatively.

We noted above that at the time of the formation of the discontinuous solution the tangential velocity profile takes the form of a caustic. In the neighborhood of the caustic tangential velocity will experience a strong change. Let $\delta\lambda$ be the difference of the angular moment in the initial tangential velocity distribution over the scale $\ell/4$:

[TABLE]

Here, we neglected a small deviation from Keplerian distribution of angular momentum in the expression (15). The change of tangential velocity at the caustic is equal to $\delta\lambda/r_{\mathrm{sh}}$ or, in the units of sound velocity,

[TABLE]

In the reference frame connected with the surface of discontinuity, configuration of the flow resembles a laminar boundary layer of the „planar jet“. Indeed, the typical width of the caustic is much smaller than the wavelength $\ell$, which in the wave zone can not exceed semi-thickness of the disk. At the same time, in the tangential direction the flow is fairly homogeneous. It is well known that the flow of this type becomes unstable if the Reynolds number $\mathsf{Re}\approx 3.7$ (see, for example, [29]). With increase of the Reynolds number, the range of the wave numbers of perturbations which cause the instability of the boundary layer expands without limits. For a typical accretion disk, one can estimate that $\mathsf{Re}\sim 10^{10}$ [7]. Even if we determine the Reynolds number from the value of the tangential velocity difference (58), its estimate decreases by two or three orders of magnitude only. This allows to claim that the region of the caustic in the distribution of the tangential velocity is unstable with respect to perturbations of practically any scale. We note that the region of the caustics is located at both sides of the density and radial velocity jump. Since the jump is a shock wave, it propagates in the gas with supersonic velocity. By virtue of this, small-scale perturbations can propagate behind the shock only [6], and, therefore, the instability can develop only there. Given that accretion disks are typically have huge Reynolds numbers, instability of the distribution of tangential velocity in the caustic may well become a source of turbulence in the disk. Then the quantity $\alpha$ (58) should be interpreted as an estimate of the Shakura-Sunyaev coefficient [2].

Expression (58) is written down for the waves that move from the center (sign „$+$“) and to the center (sign „$-$“). If the perturbation propagates toward the outer region of the disk, the value of $\alpha$ increases with increasing amplitude of the initial perturbation $\mathcal{M}_{\ast}$, but it can not exceed $1/8$ in principle due to condition (32). Let us consider as an example a protoplanetary disk around binary system with the total mass $M_{\mathrm{tot}}=2M_{\odot}$, orbital separation $A=10R_{\odot}$ and the inner disk radius $r_{\mathrm{in}}=2A$. We set the temperature of the disk equal to $10^{4}$ К. Let assume that the wave moves from the inner boundary of the disk to the outer one (in this case it is necessary to assume $r_{\ast}=r_{\mathrm{in}}$). Figure 2 shows the plots of the dependence on time of the radius of formation of the discontinuity and the amplitude of the jump of tangential velocity on the spatial scale and Mach number of the initial perturbation. It is seen that for $\mathcal{M}_{\ast}\lesssim 10^{-2}$ the value of $\alpha$ weakly depends on the initial wavelength of perturbation, but increases with the amplitude as $\alpha\sim\mathcal{M}_{\ast}$. The radius, at which the discontinuity is formed and the time of its formation rapidly grow with increasing wavelength, i.e., long-wave perturbations of small amplitude can reach the outer boundary of the disk, without having had time to form a discontinuity. Perturbations of larger amplitudes, naturally, quite quickly form high intensity jumps, where $\alpha\lesssim 1/8$.

Scenario of formation of discontinuities in the case when the wave comes from the outer side of the disk, is somewhat different. Let consider as an example a circumstellar disk around a white dwarf with the mass $M=1M_{\odot}$. Let set the inner radius of the disk equal to $r_{\mathrm{in}}=0.03R_{\odot}$ and the outer one to $r_{\mathrm{out}}=0.3R_{\odot}$, the temperature of the disk to $10^{4}$ К; $r_{\ast}=r_{\mathrm{out}}$. The parameters of the discontinuities that form in such a disk are shown in Fig. 3. Unlike the first case, here the time of propagation of the wave is limited by the time it intersects the wave zone (a situations when the discontinuity forms strictly at the boundary of the wave zone are marked in Fig. 3 by asterisks). As the wave propagates into the disc, sound velocity increases. This explains a large gradient of the dependences $r_{\mathrm{sh}}(\mathcal{M}_{\ast})$ and $\alpha(\mathcal{M}_{\ast})$ in the neighborhood of the boundary of the wave zone. The value of $\alpha$ at this point attains its maximum value $1/8$. Perturbations of a small initial amplitude, as a rule, lead to higher values of the coefficient $\alpha$ than in the first case, up to two orders of magnitude. Wave region of small-scale perturbations can cover entire disk. If the amplitude of such perturbations is small, they can reach the inner boundary of the disk, before they form a discontinuity.

6 Conclusion

In this paper we studied the evolution of nonlinear perturbations in a single class of sub- Keplerian accretion disks. It was assumed that disks can be exposed to external action at inner and outer borders. The first case corresponds to protoplanetary disks around binary systems of young stars of T Tau type; the source of impact can be detached shock waves. The second case can be realized in accretion disks of close binary stars, where perturbations can be generated by the region of interaction of the accretion flow and the disk.

It is found that at the front of a radial longitudinal perturbation wave a discontinuity of the density and radial velocity forms with time, which is a shock wave. Distribution of the tangential velocity in this case takes the form of a caustic. Position of the peak ofthe latter corresponds to the position of the shock wave. Formation of the shock wave occurs due to the well-known mechanism of the „breaking“ of the profile of the nonlinear perturbation. The value of the jump of the tangential velocity at the caustic depends on the amplitude and length of the perturbation wave. In some cases it can reach about $0.1$ of the sound velocity. Such a configuration is unstable and disintegrates. The natural consequence of the decay is turbulization of the accretion disk. The estimates of the Shakura-Syunyaev parameter for the accretion disks — $\alpha<0.1$ — in both considered cases indicate a substantially subsonic turbulence.

7 Acknowledgments

E.P. Kurbatov was supported by the Russian Foundation for Basic research (contracts No. 14-29-06059 and 15-02-06365).

8 Appendix A

After substitution $\eta=\rho_{0}\xi$, Eq. (29) for displacement of the gas elements becomes:

[TABLE]

We do not provide boundary conditions for this equation, but seek solution in the form of propagating perturbance. Therefore, we may apply the standard method of separation of variables:

[TABLE]

Then, denoting by $-\omega^{2}$ integration constant, we get (dot and prime mark denote differentiation over time and spatial coordinate, respectively)

[TABLE]

Let write down solution for the temporal part as

[TABLE]

For the spatial part we get an equation

[TABLE]

Let insert to the last Eq. radial profiles of the density and sound velocity obtained in subsection 3.1. Then, after substitution of the variables, $x=r/r_{\ast}$, $R(r)=X(x)$, and introduction of denotations:

[TABLE]

we obtain:

[TABLE]

Solution of this equation will be expressed via Hankel function of the first type:111Hankel functions are related to the Bessel functions of the first and second type: H $H_{\nu}=J_{\nu}+iY_{\nu}$. [30, с. 169]:

[TABLE]

where $C$ is a complex constant; the index of Hankel function is equal to

[TABLE]

For larger arguments of Hankel function, $\mu x^{3/2}\gg 1$, one has $H_{\nu}(z)\approx(\pi z/2)^{-1/2}\operatorname{exp}[iz+i(\dots)]$, and then

[TABLE]

where $i(\dots)$ — insignificant phase addition.

Note, if solution has to have wave character, it is necesarry to satisfy the following inequality:

[TABLE]

or

[TABLE]

For low $x$ solution (67) very rapidly aperiodically approaches zero. Figure 4 shows approximately the dependence $X(x)$.

Exact solution for linear perturbation may be obtained also for full set up of the problem, without neglect of the „geometrical“ term in the continuity equation (17). In this case, the value of displacement and density will be related as

[TABLE]

After substitution $\eta=r\rho_{0}\xi$ we obtain an equation analogous to (59):

[TABLE]

From this equation, by means of separation of variables and substitutions listed above, it is possible to get an equation for $X$, which differs from the last Eq. by the sign in the front of the second term only:

[TABLE]

Its solution [30]:

[TABLE]

where the index of hankel function is

[TABLE]

This solution does not differ essentially from Eq. (67), obtained in quasi-Cartesian approximation.