1 Introduction.

Disk Accretion Driven by Spiral Shocks


Spiral density waves are known to exist in many astrophysical disks, potentially affecting disk structure and evolution. We conduct a numerical study of the effects produced by a density wave, evolving into a shock, on the characteristics of the underlying disk. We measure the deposition of angular momentum in the disk by spiral shocks of different strength and verify the analytical prediction of Rafikov (2016) for the behavior of this quantity, using shock amplitude (which is potentially observable) as the input variable. Good agreement between the theory and numerics is found as we vary shock amplitude (including highly nonlinear shocks), disk aspect ratio, equation of state, radial profiles of the background density and temperature, and pattern speed of the wave. We show that high numerical resolution is required to properly capture shock-driven transport, especially at low wave amplitudes. We also demonstrate that relating local mass accretion rate to shock dissipation in rapidly evolving disks requires accounting for the time-dependent contribution to the angular momentum budget, caused by the time dependence of the radial pressure support. We provide a simple analytical prescription for the behavior of this contribution and demonstrate its excellent agreement with the simulation results. Using these findings we formulate a theoretical framework for studying one-dimensional (in radius) evolution of the shock-mediated accretion disks, which can be applied to a variety of astrophysical systems.

Subject headings:
accretion, accretion disks — protoplanetary disks — planet-disk interactions — hydrodynamics — shock waves

1. Introduction.

Evolution of astrophysical disks is caused by the redistribution, gain or loss of angular momentum of the disk fluid, ultimately resulting in accretion onto the central object. Angular momentum of the fluid elements can change as a result of internal or external stresses. The former usually involve turbulence produced by operation of some instability in the disk. The most common example of such an instability is the magneto-rotational instability (MRI, Velikhov 1959; Chandrasekhar 1960; Balbus & Hawley 1991), but it may fail to operate in weakly ionized protoplanetary disks (Turner et al., 2014). Other processes, such as the gravitoturbulence (Gammie, 2001; Rafikov, 2015), vertical shear instability (Urpin & Brandenburg, 1998; Stoll & Kley, 2014), etc. have been proposed to explain the observed evolution of cold protoplanetary disks. External stresses can be naturally exerted on the disk via outflows (Blandford & Payne, 1982; Konigl, 1989), resulting in the loss of angular momentum from the system and mass accretion.

Another potential driver of the disk evolution could be the global density waves, ubiquitous in astrophysical disks. Previously, observations using the technique of Doppler tomography revealed presence of spiral waves excited by the gravity of the donor star in accretion disks of cataclysmic variables (Steeghs et al., 1997; Marsh & Horne, 1988). Global asymmetric emission pattern has also been found in a compact gaseous debris disk of a metal-rich white dwarf SDSS J122859.93+104032.9 (Manser et al., 2016). Recent observations of protoplanetary disks using direct imaging in optical and near-IR (Benisty et al., 2015), as well as sub-mm interferometry with ALMA (Christiaens et al., 2014; Pérez et al., 2016), show spiral arms to be a common phenomenon. They are believed to be excited by the gravity of embedded massive planets (Boccaletti et al., 2013; Collins et al., 2009; Garufi et al., 2013; Grady et al., 2013; Muto et al., 2012), stellar companions (Wagner et al., 2015; Dong et al., 2015a; Benisty et al., 2017), or by gravitational instability in a massive disk (Meru et al., 2017). Global spirals can also be driven by the non-uniform illumination of the disk by the central star (Montesinos et al., 2016), caused by the shadowing of the outer disk by an inclined inner disk (Casassus et al., 2015; Marino et al., 2015).

Spiral density waves induced by the gravity of a massive companion have also been routinely seen in numerical studies of disks (Sawada et al., 1986b, a), both in pure hydro (Savonije et al., 1994) and in the MHD simulations (Ju et al., 2016). Numerical simulations reveal that the high-amplitude spiral density waves rapidly evolve into spiral shocks, as the non-linear effects play important role in their propagation (Rafikov, 2002a; Dong et al., 2015a, b; Zhu et al., 2015). Recent work of Ju et al. (2016, 2017) on cataclysmic variables, subsequently extended to circumplanetary disks (Zhu et al., 2016), showed that spiral shocks can be important contributors to the angular momentum and mass transport even in well ionized disks with fully developed MRI. A similar conclusion was reached by Jiang et al. (2017) in their study of super-Eddington accretion around supermassive black holes.

On the theoretical side, properties of spiral shocks in accretion disks were first considered under the assumption of self-similarity, starting with the work of Spruit (1987), Spruit et al. (1987), and Larson (1990). More recently, self-similar shock solutions in disks were constructed by Hennebelle et al. (2016) with the goal of explaining spiral waves observed in simulations of disks with the envelope infall (Lesur et al., 2015). A limitation of the self-similar shock solutions is that they apply only to special disk models with certain scale-invariant profiles of the surface density. Real astrophysical disks need not be described by these idealized models.

Goodman & Rafikov (2001) explored the evolution of weakly nonlinear density waves excited in the protoplanetary disks by the low-mass planets. Their work, cast in the local, shearing sheet approximation, provided a formalism for describing the wave propagation, shock formation and its subsequent dissipation. Goodman & Rafikov (2001) have shown, in particular, that an ensemble of density waves launched by a population of planets can lead to evolution of a protoplanetary disk on Myr timescale. This work has been subsequently extended by Rafikov (2002a) to global density waves (not limited to the vicinity of the perturber), making wave evolution formalism applicable to disks with arbitrary distributions of temperature and density. This framework was then used in Rafikov (2002b) to consider gap opening by planets via the nonlinear dissipation of their density waves. Numerical work of Li et al. (2009), Yu et al. (2010), Dong et al. (2011), Duffell & MacFadyen (2012), Li & Li (2013), Fung & Chiang (2017) have largely confirmed these analytical results.

More recently, Rafikov (2016) analytically calculated the effect of a spiral shock of arbitrary strength (i.e. not only weakly nonlinear) on the disk, through which it propagates. He showed that shocks of even moderate strength can easily drive mass accretion through the disk at rates comparable to due to internal stresses (in agreement with the simulations of Ju et al. 2016). These analytical predictions for the shock-driven were confirmed by Ryan & MacFadyen (2017), who characterized shock strength by measuring the entropy jump across its front, as suggested in Rafikov (2016).

Observationally, the most accessible characteristic of the shock is not the entropy jump, but the density contrast across the shock. Even that quantity may be difficult to measure, especially in scattered light observations of protoplanetary disk, when the brightness contrast at the shock is more sensitive to the corrugations of the disk surface than to the variations of the surface density. Nevertheless, numerical simulations with radiation transfer post-processing can be used to calibrate the relation between the brightness contrast and the density jump, such as recently done by Dong & Fung (2017). Also, optically thin sub-mm observations using ALMA should be able to provide a direct probe of the density contrast across the shock.

This observational connection provides a motivation for our present work. Here we seek to numerically verify the analytical predictions of Rafikov (2016) formulated in terms of the (observable) density contrast across the shock. This provides a new way of understanding the shock-driven disk evolution, which is complementary to that used by Ryan & MacFadyen (2017).

Our work is structured as follows. In §2 we summarize the analytical background for our calculations, while in §3 we describe our numerical setup. The typical results of our simulations are discussed in §4, and in §5 we provide an exhaustive comparison between the numerical results and theoretical predictions of Rafikov (2016). The emergence of the secondary shocks in our simulations and their effect on disk evolution are described in §6. In §7 we provide an analytical prescription for computing the time-dependent angular momentum contribution emerging in our simulations, which allows us to formulate a general framework for following the shock-driven disk evolution in §8. We discuss and summarize our results in §9 and 10, respectively.

2. Basic equations.

We consider a two-dimensional (2D) fluid disk and characterize its properties in polar coordinates via the surface density and radial and azimuthal velocities and . We are interested in the effects of a spiral shock propagating through the disk on the angular momentum and mass transport. For simplicity, the disk is assumed to be unmagnetized and we neglect the gravitational effect of any external perturbers as well as the disk self-gravity. This does not mean that the density wave propagating through the disk cannot be launched via the gravitational coupling to an external perturber such as a planet or stellar companion. This assumption simply implies that we are studying the disk far from the wave excitation region, so that the torque induced on the density wave by the perturber’s gravity can be ignored.

Angular momentum conservation in a disk can be described quite generally via the equation (Balbus & Papaloizou, 1999; Ju et al., 2016)


where implies integration over , for any function . As the fluid motion does not depart too far from axisymmetry, it is convenient to define azimuthal velocity perturbation , where is the angular velocity in the absence of perturbation. Then one can easily show that the equation (1) reduces to


where is the mass accretion rate (defined to be positive for inflow), is the specific angular momentum, and is the Reynolds stress.

Note that the choice of is to some extent arbitrary — changing affects both the first term in the left hand side and the right hand side ( depends on ) of equation (2) in such a way that the relation (1) always holds true. It is customary to assume (Balbus & Papaloizou, 1999; Ju et al., 2016), and we will adopt this convention in our numerical work as well.

Disks evolving predominantly under the action of internal stress feature non-zero , which drives their evolution and mass transport (i.e. makes non-zero) as demonstrated by equation (2). Internal stress can be provided by turbulence driven by some instability in the disk, e.g. MRI.

Situation is different in disks evolving under the action of the global spiral shocks. In this case the angular momentum is directly injected into the disk at the shock fronts. Between the shock fronts the angular momentum is conserved in the absence of internal stresses. Even though such pattern of angular momentum injection is highly localized in azimuthal angle, the long-term effect of the shocks on the disk fluid can still be described via the azimuthally averaged quantities. In particular, we can assume that shock damping transfers the momentum carried by the density wave to the disk fluid at the rate per unit radial distance, where is the angular momentum flux carried by the wave.

In shock-mediated disks, the angular momentum deposition rate replaces the stress term in the right hand side of equation (2), so that the angular momentum evolution is governed by


As a result, the mass accretion rate in a shock-mediated disk can be expressed as


This expression generalizes a similar result of Rafikov (2016), who neglected the last time-dependent term in equation (4) finding (Lynden-Bell & Pringle, 1974)


This simplification is legitimate when the disk is in steady state, which is generally not true ( in particular, in our work, see §7). For that reason, here we will use a more complicated expression (4).

Rafikov (2016) demonstrated theoretically that the angular momentum deposition rate due to damping of a spiral shock with azimuthal wavenumber and pattern speed is given by


where is the pre-shock value of the isothermal sound speed at the radius , and is the pre-shock surface density. Also, an auxiliary function of the ratio of the post-shock () to pre-shock () pressure is given by


for an ideal gas with polytropic index . In practical applications one can often approximately set and in equation (6) to their mean (i.e. azimuthally averaged) values at a given radius.

Another convenient way to express irreversible heating is via the density contrast across the shock ( and being the post- and pre-shock values of gas density):


Note that in our 2D treatment of the disk should be replaced with .

In the case of an isothermal equation of state () one finds from expression (7), taking the limit , that


where is the relative density (or pressure) jump across the shock. This expression is valid for arbitrary amplitude of , and agrees with similar result found in Belyaev et al. (2013).

At small wave amplitudes () all these expressions reduce to the well-known cubic scaling (Spruit, 1987; Larson, 1990; Savonije et al., 1994)


However, at high wave amplitudes this approximation becomes inaccurate, as we show later (see §5.1).

The main goal of this work is to numerically verify the validity of equations (4), (6)-(9). We do this in §4-7.

3. Numerical setup

We run a set of two-dimensional numerical simulations in cylindrical polar coordinates using publicly available code Athena++4 (Stone et al. 2017, in preparation), the grid-based code, which solves equations of gas dynamics using high-order Godunov methods. This code is a new version of a popular astrophysical gas dynamics code Athena (Gardiner & Stone, 2005, 2008), and features improved performance and scalability as well as more flexible choice of coordinate systems. For all our simulations we use HLLE solver (Einfeldt, 1988). Athena++ was extensively used in simulations of protoplanetary disks (Dong et al., 2015b; Zhu et al., 2015; Arzamasskiy et al., 2017) and cataclysmic variables (Ju et al., 2016, 2017), and was proven to accurately reproduce the properties of density waves. To provide a direct comparison with the analytical predictions of Rafikov (2016), we consider a gaseous disk with no explicit viscosity (we describe the effects of numerical dissipation in §5.3).

The active domain of our simulations , i.e. it is limited to only half the disk in direction. We do that to increase resolution for a given computational cost. As we show in §5.3, high resolution is crucial in order to demonstrate agreement between theory and simulations. Our typical resolution is 447 cells per scale-height at the outer radius (or for brevity).

We typically run our simulations for , where is the orbital period at the outer radius. This is sufficient for the global wave pattern to fully develop inside the domain and reach a steady state. However, this not long enough for the disk material to get radially redistributed by shocks for const to be achieved, i.e. the disk itself does not relax to a global equilibrium (cf. Ryan & MacFadyen 2017). For that reason the first, time-dependent term in equations (2)-(3) plays important role in the angular momentum balance. This has important consequences discussed in §7.

Our strategy consists of running a number of simulations, in which only one parameter is varied as compared to a certain fiducial run. This allows us to get a clear picture of the effect of different physical inputs on the results. The numerical parameters used in our fiducial run are listed in Table 1.

Figure 1.— Two-dimensional map of the (logarithm of the) surface density perturbation emerging in one of our typical simulations at ( is the orbital period at the outer radius ) having an amplitude of at the outer boundary of the computational domain. This simulation features a fiducial set of parameters, namely a resolution of , aspect ratio at the outer radius, radially constant background surface density and temperature profiles, and the isothermal equation of state for the gas (see Table 1). The red dashed curve indicates the location of the spiral wake predicted by the linear (low-amplitude) theory (Rafikov, 2002a; Ogilvie & Lubow, 2002). The actual location of the wave crest (coincident with the shock position) is shown by the white dashed curve. It deviates from the linear theory prediction because of the non-linear effects (Rafikov, 2002a; Zhu et al., 2015).

3.1. Initial conditions

We initialize the disk with the equilibrium model described in Nelson et al. (2013), which has power-law profiles of the surface density and temperature:


where and correspond to the outer radius of the simulation domain. The hydrodynamic equilibrium implies that the angular velocity of the disk


where is the Keplerian frequency, and is the local isothermal sound speed. It is customary to describe the temperature of the disk using the aspect ratio:


For our fiducial run we use the uniform density profile with , as well as the globally isothermal equation of state , . The non-uniform density profiles are considered in §5.5, and we study other equations of state in §5.4. We adopt as our fiducial disk aspect ratio. We modify the aspect ratio and study its effect on the results in §5.2.

Parameters Values
Table 1Simulation set up

3.2. Boundary conditions and wave triggering

To launch the density wave we impose a boundary condition in the form of a density perturbation at the outer radius of the disk. The shape of the perturbation in sets the initial profile of the density wave, and is given by


where sets the rotation of the perturbation pattern at angular speed , and describes the initial width of the profile. We use for all of the simulations (this azimuthal range is resolved by about 90 grid points, see Table 1). The pattern speed is usually chosen to be for simplicity. We study the effects of nonzero in §5.6.

We set the velocity perturbation in and directions to zero at the outer boundary. At the inner disk boundary we reset all the quantities to their initial values.

In our simulations we vary the proportionality coefficient in equation (15) to have different perturbation amplitude . As the simulation evolves, the density perturbation propagates into the active domain in the form of a spiral density wave, which rapidly turns into a shock. The steep profile for assures that the shock develops very close to the outer boundary.

We do not expect our results to be affected by the details of the implementation of our boundary conditions, as long as the wave shocks close to the outer boundary. However, some fine features, like the appearance of a secondary shock (see §6), may be more sensitive to the boundary conditions. Careful examination of this issue is deferred to future study.

Since we simulate only one half of the disk, we use periodic boundary conditions to connect the values of fluid variables at and . As a result, the developing perturbation pattern features two identical density waves separated by in azimuth, see Figure 1. When computing angular momentum contributions using equation (2), we account for both shocks.

Figure 2.— Evolution of the profile of the density perturbation shown in Figure 1. Azimuthal cuts through the spiral wake at different radii are shown, referenced to the azimuthal angle corresponding to the linear prediction for the density wake location (Rafikov, 2002a; Ogilvie & Lubow, 2002). One can see that nonlinear evolution and damping cause continuous distortion of the wave profile, making it broader and reducing its amplitude.

3.3. Shock detection

To calculate the density jumps across the shocks we employ a simple algorithm to detect their location. At each radius, we calculate and identify the shock with the azimuthal region within which this derivative is greater than 1. Having done that, we calculate the number of shocks, the jump conditions across the shock, as well as the azimuthal width of the wake.

As shock has very steep density profile, the location of the shock is not very sensitive to the details of the detection algorithm5. We have verified by eye that our algorithm correctly reproduces both the location and the amplitude of the shock at all radial positions. The shock profile in coordinates resulting from this algorithm is illustrated in Figure 1.

4. Results

We now provide a brief description of certain characteristic features of our simulations, before giving a more in-depth, systematic analysis of their results in §5.

Figure 1 is a 2D color map of the surface density perturbation in one of our simulations with the fiducial set of parameters listed in Table 1. This particular run features surface density perturbation with the amplitude of at the outer edge of the simulation domain, see §3. One can see that the perturbation propagates through the radial extent of the disk in the form of a coherent spiral wave all the way to the inner boundary. In our runs we do not observe the development of the spiral wave instability, previously reported by Kim & Kim (2014), Bae et al. (2017), and Sormani et al. (2017). It is not clear whether this outcome is caused by the lack of an external gravitational perturber in our simulations, or their relatively limited radial extent.

Figure 3.— Radial profiles of the different angular momentum contributions appearing in equations (3) and (16), obtained from simulations with different initial wave amplitudes at the outer boundary: (a) and (b) . Red solid, blue dashed and green dotted lines show the advective term , stress term , and theoretical prediction for the stress contribution, correspondingly. The latter is computed via equations (6)-(9) using the density jump at the shock measured from simulations. The orange curve shows the time-dependent term , while the grey curve is the residual term , see equation (16).

The red dashed curve in Figure 1 shows the linear prediction for the shape of the spiral wake based on the calculation in Rafikov (2002a) [see equation (36) of that work]. One can see that the actual location of the wave crest deviates from the linear approximation in such a way as to make the spiral less tightly wound in the central part of the disk than the linear theory would predict. This deviation is caused by significant nonlinear distortion and broadening of the azimuthal profile of a high amplitude wave (Rafikov, 2002a). Radial profile of the surface density jump across the shock relative to the local background density , characterizing the evolution of the wave nonlinearity, is shown by the blue curve in Figure 5a, which is discussed further in §5.1.

Figure 4.— Radial profiles of the angular momentum contributions appearing in equations (3) and (16), obtained from simulations with different initial wave amplitude at the outer boundary (shown by different colors, see legend in the right panel). (a) Solid, dashed and dotted lines show the advective term , stress term , and theoretical prediction for the stress contribution, correspondingly (as in Figure 3 but on log vertical scale). (b) The absolute value of the time-dependent term (solid curves). Colored dots of corresponding color show theoretical prediction for the time-dependent term computed using equation (19), which is in excellent agreement with the numerical results. See §7 for more details. (c) The residual term defined by equation (16), which is caused by numerical dissipation in our simulations.

The evolution of the azimuthal profile of the wake is illustrated in Figure 2, which displays azimuthal cuts (at several radii) of the surface density (normalized by the density at the outer boundary) around the azimuthal angle corresponding to the linear prediction for the wake position. One can clearly see that, as a result of the nonlinear effects, the wake shape gets continuously distorted, resulting in the decay of its amplitude and profile broadening as the wave propagates towards the disk center. The importance of properly accounting for the wave non-linearity (as opposed to using purely linear results) for explaining the openness of the spiral arms driven by massive perturbers in the protoplanetary disks has been previously emphasized in Dong et al. (2015b) and Zhu et al. (2015).

As the main focus of our paper is on verifying the analytical predictions of Rafikov (2016) for the shock-driven mass and angular momentum transport, we extract from our simulations and analyze a number of relevant quantities. In particular, Figure 3 illustrates the radial profiles of the different angular momentum contributions appearing in the equation (2), directly measured in the simulation with the fiducial disk parameters and the perturbation amplitudes of 1 and 0.0625 at the outer radius. In the order that they appear in this equation, we call these contributions the time-dependent term (as it involves the time-derivative of disk variables) , the advective term , and the stress term . To ensure better statistics, all terms are averaged over after the wave pattern fully develops. It is clear that the radial profiles of the advective and stress terms differ in shape, which is caused by the substantial contribution coming from the time-dependent term. The difference in especially pronounced close to the outer disk radius, where the time-dependent term is most significant.

One can also see that the three terms do not sum up to zero, as they should according to the equation (2). The remaining non-zero residual (thin grey line) is explicitly defined as


see equation (3). It is roughly independent of radius and in simulation with has a magnitude comparable to other angular momentum contributions. The fact that is not equal to zero is caused by the numerical dissipation intrinsic for our finite-resolution runs; this is discussed in more details in §5.1 and 5.3.

In Figure 3 we also show the analytical prediction (with negative sign, so it can be directly compared with advective and stress terms) given by the equation (6) and calculated using the value of the density jump at the shock front measured in our simulations. Figure 3b reveals that has a radial profile similar to the stress term, but with a vertical offset approximately equal to the numerical residual term .

Figure 3a is identical to Figure 3b, except that it corresponds to a run with a higher amplitude of the perturbation at the outer edge, . Stronger perturbation clearly results in a significant increase in the amplitude of all physical angular momentum contributions, compared to Figure 3b. At the same time, the numerical residual does not increase nearly as much, rendering its contribution to the angular momentum balance (i.e. equation (2)) insignificant. This results in a considerably better agreement between the theoretical and the stress term at high wave amplitudes. Also, the time-dependent term is still important and leads to a significant difference between the profiles of the advective term and both the stress term and the theoretical prediction. We provide a more systematic discussion of the role of the amplitude on the wave propagation and associated transport properties in §5.1.

5. Comparison with theory

Figure 5.— (a) Radial profiles of the density jump across the shock measured in simulations with different values of the initial wave amplitude (indicated in panel). (b) Radial profiles of the auxiliary function (solid line), given by equation (9) for an isothermal equation of state and density jumps from the panel (a). Dashed lines show the weakly nonlinear approximation given by equation (10). This approximation clearly overpredicts the effect of the shock on the disk, especially at high wave amplitudes.

We now provide a more detailed comparison between the theory presented in Rafikov (2016) and our hydro simulations. The theory predicts that the effect of shocks on the disk depends on many physical parameters: strength of the shock, aspect ratio of the disk, radial profiles of the background density and temperature, equation of state, and so on. Our goal is to thoroughly test these dependencies. In addition, we will explore the sensitivity of our results to numerical parameters, such as the resolution of our simulations. This will help us formulate the conditions, under which one could expect simulations to reliably reproduce the actual effect of the global spiral shocks on disk evolution.

5.1. Variation of the wave amplitude.

We start by exploring the effect of the wave amplitude on the correspondence between the theory and simulations. In Figure 4 we plot different angular momentum contributions appearing in equation (2) as a function of radius, for different values of the wave amplitude at the outer radius.

In the left panel of this Figure we plot the stress term together with the theoretical prediction , as well as the advective term . One can see a very rapid increase in the magnitude of all angular momentum contributions with the amplitude of the initial surface density perturbation. The top and bottom sets of curves correspond to the values of different by a factor of 16, but the associated angular momentum contributions differ by more than three orders of magnitude near the outer edge.

Similar to the pattern described in §4, we see good agreement between the and . Not only does the theory correctly predict the variation of the shock-driven stress with , the two curves also closely follow each other in radius, with just a small vertical offset separating them. The offset becomes less pronounced for larger . Note that theoretical curve always underpredicts the numerical result, but only at the level of several tens of per cent or less at high wave amplitudes. The reasons for this discrepancy are discussed in §9.

In the linear regime () the agreement between and is much worse, as already mentioned in §4. To illustrate the reasons, in the right panel of Figure 4 we display the variation of the residual with the wave amplitude. One can see that at the lowest amplitudes is comparable to the stress and advective terms in magnitude. However, as increases, the residual term grows much slower than the physical angular momentum contributions. As a result, the quantitative agreement between and is much better for more nonlinear waves (higher ). We discuss the origin and behavior of further in §5.3.

Similar to Figure 3, the left panel of Figure 4 also shows a significant discrepancy between the advective and the stress terms. The qualitative difference is present at all wave amplitudes, and is especially pronounced in the outer disk. As in §4, this difference is naturally explained for all values of by the large time-dependent contribution . Its absolute value is shown in the middle panel of Figure 4. It is comparable to the stress and advective contributions in the outer disk, showing the importance of accounting for this angular momentum term in general. This implies that in our calculations we cannot accurately predict the radial behavior of just by using the steady state equation (5) and the shock dissipation prescription (6). Determination of the accretion rate throughout the disk in general requires the use of full equation (4), properly accounting for the time-dependent term, which is discussed further in §7.

In Figure 5a we display the radial profiles of the relative surface density perturbation for different . One can see a clear variation in the behavior as the wave nonlinearity changes. At the lowest amplitudes remains roughly constant as the wave travels inwards. However, at high values of density perturbation rapidly decays with : for the value of drops by about a factor of 3 as the wave propagates from to . This decay of the perturbation is caused by the much stronger nonlinear wave damping at higher values of . As a result, the angular momentum terms in Figure 4a also show much faster inward decay for higher .

Figure 6.— Same as Figure 4 but now illustrating the effect of changing disk scale height (values indicated in panel (c)) on the angular momentum flux contributions. Resolution per scale height is the same () for these two simulations. Angular momentum deposition rate by the shock shows faster decay with the distance travelled from the outer edge in a colder disk. Note the rise of the shock-induced stress in the inner part of the disk with , explained by the appearance of a secondary shock, see §6.

We also note a very sharp drop in the amplitude of the surface density perturbation at the outer boundary of the disk. This artifact of our boundary conditions arises because we do not initialize velocity perturbation consistently with density perturbation at the boundary. As a result, the actual amplitude of the wave close to the boundary is slightly different from . However, this does not affect the comparison between theory and simulations as is always computed using the actual local value of the surface density perturbation.

Figure 5b illustrates how the simple cubic approximation (10) describes the effect of the shock compared to the fully nonlinear prescription (7)-(9); we plot the full given in the isothermal case by the equation (9) as well as the small amplitude approximation (10). The two clearly agree with each other quite well at the low values of . However, as the wave nonlinearity increases, the cubic approximation considerably overpredicts the actual effect of the shock on the disk; the discrepancy is close to a factor of 2 for .

5.2. Sensitivity to .

Next we vary disk aspect ratio and explore the effect on agreement with theory. In Figure 6 we compare different angular momentum contributions for two simulations with , in which we vary from to . In doing so we kept the resolution per scale height the same, meaning that we had to decrease by a factor of 2 the number of grid points in every dimension in the case .

Comparing different sets of curves in Figure 6 corresponding to different , one can make several observations. First, irrespective of the value of , our simulations demonstrate good agreement between the theoretical () and numerical () shock-induced stress terms. The advective term clearly deviates from both because of the time-dependent contribution shown in Figure 6b.

Second, reduction of the disk scale height by a factor of 2 results in a dramatic reduction (by more than an order of magnitude at some radii) of all physical angular momentum contributions. Close to the outer edge of the disk this difference is consistent with our theoretical prediction (6) showing the explicit quadratic dependence of the shock-driven angular momentum flux on (or, equivalently, ).

Third, the rate at which the density wave transfers angular momentum to the disk fluid drops faster as the wave propagates in colder disk. This is caused by the faster evolution of the wave for lower and, subsequently, more rapid decay of the angular momentum carried by the wave (also implying faster decay of across the shock). This effect is similar to the faster decay of the waves with higher initial amplitude, obvious in Figures 4a and 5a. For that reason, at radii of the values of derived for two different differ more dramatically than near the outer edge of the disk.

Figure 7.— Same as Figure 4 but now illustrating the effect of changing resolution on the agreement between the theory and numerics.

Fourth, in colder disk the behavior of the stress term exhibits a qualitative change: its decay as the wave propagates inwards switches to growth at . This change is caused by the faster evolution of the wave in a colder disk, resulting in emergence of a secondary shock. We will discuss this phenomenon in more detail in §6, making only a couple of comments here. The new shock structure provides additional dissipation at its front, causing growth of as decreases (ultimately this new contribution to the stress term will also start to decay). In a hotter disk with we would find the same change in the behavior but at smaller radii (which are outside of our computational domain) because of the slower evolution of the wave. Nevertheless, regardless of the more complicated picture of the wave evolution in colder disks, our theory still reproduces numerical results very well (apart from a small vertical offset).

Figure 8.— Residual term plotted as a function of the wave non-linearity — wave amplitude at the outer edge of the computation domain. Different curves correspond to different sets of the disk () and simulation (resolution per scale height ) parameters, indicated on the panel. Dashed line shows scaling . This figure illustrates that the residual goes down with both the disk aspect ratio and the numeical resolution.

5.3. Effects of changing resolution.

We now discuss the role of numerical effects on our results. In Figure 7 we display individual contributions to the angular momentum budget obtained in simulations with three different resolutions.

Figure 7a shows that the shock-driven torque term is insensitive to resolution. Its theoretical counterpart (computed using numerically determined at the shock) varies with resolution, although relatively weakly, giving better agreement with simulations at higher resolution. Both the advective term and the time-dependent contribution shown in Figure 7b are affected by varying resolution only weakly.

What does change significantly with resolution is the residual term shown in Figure 7c (it is multiplied by a factor of 10 for better representation). One can see that increasing resolution consistently reduces the magnitude of . This naturally brings us to the conclusion that the residual term is caused by the numerical dissipation due to the finite resolution of our simulations.

The effect of the numerical residual is twofold. First, it gives rise to an unphysical contribution to the angular momentum balance in the disk, which corrupts the agreement between and . This effect is discussed further in §9. Second, higher numerical dissipation at lower resolution causes unphysical reduction of the density jump across the shock. Although small, this effect artificially decreases the value of at lower resolution, which is seen in Figure 7a.

Figure 9.— Same as Figure 4, but now showing the effect of varying the equation of state (we change the adiabatic index as shown in the right panel) on the agreement between the theory and simulations.

Figure 8 provides additional information on the behavior of the residual . It clearly shows that, everything else being equal (e.g. resolution, , etc.), the value of scales roughly linearly with the initial wave amplitude. This behavior is what one would expect from the usual viscous term in the fluid equations, reinforcing our interpretation of as being due to numerical dissipation. Note that all physical shock-induced angular momentum contributions exhibit a significantly steeper dependence on , see equations (7)-(9) and Figure 4a.

Figure 8 also clearly shows how the residual decreases as the resolution goes up. At a fixed resolution is also lower for smaller , simply because the angular momentum flux carried by the wave of fixed amplitude is lower in colder disks, see Figure 6.

5.4. Variation of the EOS.

Equation of state (EOS) of the disk fluid should have an important effect for shock-mediated transport. This can be seen simply in the explicit dependence of the theoretical prediction (7) on the adiabatic index in the non-isothermal case.

In Figure 9 we compare the outcomes of three simulations with the same starting amplitude but having EOS with , , and (isothermal). One can see that wave dissipation proceeds differently depending on the value of . Despite that, theoretical correctly reproduces the variation of as changes. Note that in non-isothermal runs we must use equation (7), which apparently performs as well as the equation (9) does in the isothermal case.

Figure 10.— Same as Figure 4 but now illustrating the effect of changing the slope of the background surface density , see equation (11), on the agreement between the theory and numerics.

On the other hand, in our non-isothermal runs the advective term deviates from the stress term much stronger than in the isothermal case. This is a result of a more pronounced role of the time-dependent contribution (see Figure 9b), which is likely caused by the continuous and spatially non-uniform injection of entropy into the disk by the shock in the non-isothermal simulations. This steady growth of entropy (absent in the isothermal case) results in a more pronounced role of the time-dependent term in these runs.

5.5. Variation of the background disk properties.

Most of the simulations shown in this paper feature radially uniform surface density profile (i.e. ). To demonstrate that this assumption does not affect our conclusions, in Figure 10 we show a series of runs with isothermal equation of state and outer amplitude , in which we vary the power-law slope of the background density profile (which requires accounting for the radial pressure support in equation (13)).

Figure 11.— Same as Figure 4 but now illustrating the effect of nonzero angular frequency of the imposed wave pattern on the agreement between the theory and numerics. Color scheme is illustrated in panel (c).

One can see that the variation of profile does surprisingly little to the numerical results, as both and the advective term exhibit only slight changes. The theoretical prediction varies even less. The agreement between and remains good.

We also run several non-isothermal simulations in which we vary the radial profile of , while keeping constant. Despite the nonuniform temperature profile, we again found good agreement between the theoretical prediction (6) and the numerically determined (we do not provide the figure illustrating this).

5.6. Rotation of the perturbation pattern.

The final parameter that we varied in our simulations is the angular frequency of the imposed wave pattern. All of the simulations so far used non-rotating pattern, , stationary in the inertial frame. In Figure 11 we show the effect of a non-zero on our results, expressing in units of  — Keplerian angular frequency at the outer edge of the disk. For these runs we impose perturbation at the outer boundary (the same as in the non-rotating case ), which uniformly rotates with and , allowing us to keep corotation outside the simulation domain.

One can see that variation of has a significant effect on the simulation results. The run with even exhibits the secondary shock (similar to that in Figure 6 for ) starting at around , resulting in the increase of inward of this radius. Nevertheless, despite these modifications, our analytical predictions (6)-(9) still work well, and the agreement between the stress term and for different values of in Figure 11 remains very good.

Figure 12.— Evolution of the wave profile, similar to Fig. 2, but now for a simulation with . Faster nonlinear evolution of the wave in a colder disk results in the emergence of a secondary shock at .

6. Secondary shocks.

One intriguing feature that emerges in some of our calculations is the secondary shock. It manifests itself in the angular momentum plots via the change in the behavior of with radius: the stress contribution, monotonically decaying with decreasing radius, suddenly starts to rise. This transition is clearly seen in our Figures 6a (for ) and 11a (for ). Similar behavior has also been observed in simulations of Ryan & MacFadyen (2017).

In Figure 12 we directly trace the appearance of this feature by showing the evolution of the wave profile in a simulation with . One can see that as the primary wave propagates inwards, a region of the negative surface density perturbation gradually develops ahead of the main shock (which itself features ). This behavior can be seen in the Figure 2 as well. In both cases the right side of the negative segment of the profile gradually steepens due to the nonlinear effects (Goodman & Rafikov, 2001; Rafikov, 2002a), until it breaks at in Figure 12.

The reason why this does not happen in Figure 2 can be related to a higher . Indeed, it has been shown both locally (Goodman & Rafikov, 2001) and globally (Rafikov, 2002a) that the rate at which nonlinear evolution of the weakly nonlinear density waves proceeds is set by the time-like coordinate


where is the aspect ratio at some fixed radius in the disk and is a function, monotonically increasing with the distance travelled by the wave, that absorbs the radial dependence of (its explicit form can be deduced using equations (32), (34) of Rafikov (2002a)). For a given starting amplitude and shape of the wave, a shock forms when a certain critical value of is reached. It is clear from equation (17) that this happens earlier (i.e. after traveling a shorter radial distance from the outer boundary) in colder disks with lower .

A very similar behavior has been seen in a number of recent studies of the disk-planet interaction (Zhu et al., 2015; Fung & Dong, 2015; Dong et al., 2017; Arzamasskiy et al., 2017), which generically demonstrate the emergence of a secondary shock at radii about halfway inward from the planetary orbit. In Bae et al. (2017) the secondary (or even tertiary) shocks were suggested to be the cause of the multiple ring-like structures in the protoplanetary disks, similar to those observed in HL Tau (ALMA Partnership et al., 2015) and TW Hya (Andrews et al., 2016). The evolutionary sequence leading to secondary shocks in these studies is identical to the picture outlined above: the development of the negative segment of the wave profile, which gradually evolves into a shock due to nonlinear effects. The second stage is a robust phenomenon, analogous to the nonlinear development of the N-wave in Goodman & Rafikov (2001).

What is less clear is the origin of the negative segment of the wave profile (giving rise to this whole evolution) in the first place. Fung & Dong (2015) and Lee (2016) suggested that the secondary shocks are excited by the planetary gravity at the ultraharmonic resonances. Our results cast doubt on this interpretation, or at least suggest that possible pathways to secondary shocks are not unique.

Indeed, in our case we see the secondary shocks to appear in simulations where external gravitational potential is completely absent as the wave is driven simply by the imposed pressure perturbation sharply localized at the outer boundary of the disk. Moreover, most of our simulations feature non-rotating wave patterns, so resonances with the pattern frequency are also excluded as a culprit behind the secondary shocks. Finally, the radial location at which the secondary shock emerges in our simulations is certainly non-unique and depends on , pattern speed of the perturbation and the initial wave amplitude. This is different from some previous studies, which typically find the secondary shock to appear at of the planetary semi-major axis (Zhu et al., 2015), although some dependence of the exact location on the disk temperature and planet mass was reported in Bae et al. (2017).

Unfortunately, at the moment we do not have a good theory explaining the emergence of the negative perturbations in our simulations. We observe this phenomenon in simulations with different initial amplitudes, including the lowest we tried (), which are essentially in the linear regime. This likely excludes the nonlinear effects, such as ultraharmonics (Fung & Dong, 2015; Lee, 2016), from being responsible for the appearance of the negative feature in our simulations. Our best guess at the moment is that it may be somehow related to the way, in which we are imposing our boundary condition for the perturbation at the outer edge of the simulation domain. Future work will show whether this explanation holds.

7. Time-dependent contribution.

Our simulations clearly demonstrate the important role played by the time-dependent term in relating the angular momentum and mass transport, see equation (3). In our case this term is particularly pronounced at large radii, close to the outer boundary of the simulation domain, where its contribution to the angular momentum balance is comparable to other terms, see Figures 3, 4b and others.

To better understand the nature of this term we note that the radial component of the Navier-Stokes equation can be written as , where we neglected the inertial terms altogether. Given that our convention assumes that , one then finds that (up to terms higher order in )


due to the radial pressure support in the disk.

Plugging this expression into the definition of the time-dependent contribution to the angular momentum balance and azimuthally averaging, one finds


where is the azimuthally averaged pressure.

We show how well this approximation works in Figures 4b, 6b, 7b, 9b-11b. There we display (with dots) the behavior of the time-dependent term given by the equation (19) using the radial profile of pressure measured in the simulations. Note that even though the majority of our simulations feature uniform initial profile of and radially-constant pressure, the time-dependent term is generally non-zero. This happens because even the weak modification of the by the shock-driven radial redistribution of mass turns out being sufficient for the right hand side of the equation (19) to provide a significant contribution.

Figure 4b and other similar plots make it clear that the analytical approximation for the time-dependent term works extremely well (the agreement is somewhat worse at lower resolution). For all wave amplitudes and at all radii the theoretical prediction falls right on top of the numerically determined values of the time-dependent term (the only mild exception being the lowest amplitude case for which the numerical dissipation might play some role). This is a very reassuring finding. It implies, for example, that the knowledge of the radial profile of pressure allows us to accurately predict the radial profile of using equations (4), (6), and (19) even in non-steady disks. Another important application of this finding is discussed next.

8. Shock-driven evolution of accretion disks.

We are now in position to formulate a fully time-dependent equation for surface density evolution in a shock-mediated disk. First, using our result (4), the continuity equation describing the evolution of the azimuthally averaged disk surface density can be written quite generally as


Second, equation (19) allows us to put this expression in a closed form:


Note that we can set in this equation without loss of accuracy. Indeed, one can show that the extra term inside the square brackets arising from the deviation of from due to the non-zero can be comparable to the last term (due to the time-dependent pressure support) only when the disk evolves slowly, on viscous timescale. However, in that case both the correction to and the last term are negligible compared to by a factor of at least anyway.

Together with the analytical prescription (6) for the equation (21) represents a closed form evolution equation for the (azimuthally averaged) surface density . Disk surface density enters this equation not only through the first, explicitly time dependent, term and , but also through the dependence of upon in the last term. Equation (21) thus represents a fully self-consistent framework for following the evolution of in the shock-mediated accretion disks, provided that both the thermodynamic state of the disk (i.e. the radial behavior of the sound speed ) and the radial profile of the shock strength (i.e. ) are fully specified.

Shock-mediated evolution of an accretion disk should eventually drive its inner regions towards a (quasi-)steady state (far from the corotation region corresponding to the pattern frequency ). In this case one can set in equation (21), resulting in . Equations (4) and (6) then imply that the surface density profile in this part of the disk should converge to


where we assumed Keplerian rotation. This expression illustrates that in steady state the radial behavior of the disk surface density is fully determined by the radial profiles of the background temperature (i.e. ) and of the pressure jump at the shock .

Note that equation (21) is not limited to shock-mediated disks, but applies to other types of disks as well. Recalling that in general one could apply this equation to disks evolving due to other mechanisms of the angular momentum transport, for which the stress tensor can be specified (e.g. some kind of effective viscosity, MRI, and so on).

Equation (21) is different from the classical form (without the time-dependent term) derived in Lynden-Bell & Pringle (1974), which is broadly used to study one-dimensional viscous disk evolution. This is because in their derivation Lynden-Bell & Pringle (1974) neglected pressure support altogether and set the specific angular momentum of the disk fluid to be equal to the Keplerian value. This is fully legitimate if the disk evolves slowly, e.g. on global viscous timescale. However, our simulations demonstrate that in more rapidly evolving disks (e.g. due to initial readjustment of the disk structure to the shock-imposed torque, as in our case) the time-dependent term may not be neglected and should be properly accounted for. It is likely that this term plays important role in various astrophysical disks experiencing rapid evolution, e.g. during outburst events.

9. Discussion

Results of simulations presented in §4 and 5 clearly support theoretical calculation of the angular momentum deposition by the spiral shocks of arbitrary strength presented in Rafikov (2016). In all cases when the shock amplitude is high and angular momentum injection by the wave (i.e. the stress term) dominates over the numerical effects (residual ) we find good agreement between and . The agreement is worse in simulations with small wave amplitudes, but even then the analytical prescription of Rafikov (2016) definitely reproduces at least the qualitative behavior and the magnitude of .

Our results strongly emphasize the importance of high resolution that the simulations trying to account for the disk shock phenomena must have. As demonstrated in Figure 3, at low wave amplitudes the physical contributions to the angular momentum balance become comparable with the numerical residual term (or even subdominant). At high wave amplitudes, in the fully nonlinear regime, situation improves because the magnitude of the physical angular momentum contributions grows faster than (Fig. 8). This likely explains why Ju et al. (2016) saw the numerical residual to be close to zero in their study of high-amplitude spiral shocks in disks of cataclysmic variables.

Figure 13.— (a) Solid, dashed and dot-dashed lines represent the stress term , and the theoretical predictions for computed using entropy jump (23) and density jump (8), correspondingly. (b) The difference between the density jump measured from simulation and theoretical prediction (24) computed from the measured pressure jump. The blue dashed line on this figure represents the run with two times lower resolution (as compared to the blue solid line). Both figures use the same color scheme.

It is also clear from Figure 4 that in most cases the residual term cannot fully explain the difference between the stress term and theoretical . We believe that this difference is also caused by numerical effects. In Figure 13a we again show the angular momentum terms for simulations with adiabatic equation of state () and different wave amplitudes. Solid line shows the stress term, while other lines represent two kinds of the theoretical prediction for the stress term: one computed from the density jump (8) and shown with dot-dashed line, and another computed via the entropy jump at the shock front (dashed curve), as suggested in Rafikov (2016) and previously implemented in Ryan & MacFadyen (2017). The latter can be written in terms of the density and pressure jumps at the shock as


We see surprising disagreement between the theoretical predictions computed using two methods. On the other hand, the angular momentum flux computed using entropy jump reproduce numerical stress term very well.

The only difference between equations (8) and (23) is the connection between the pressure and density jumps (Landau & Lifshitz, 1959):


Figure 13b shows that our simulations do not reproduce this relation exactly, i.e. the actual post-shock density is slightly different from . Although the error is very small, of the order of , it causes a big effect on . One can estimate the error in as


We verified that the difference between the two methods of analytical evaluation of the shock-induced stress can indeed be approximated by equation (25) with very good accuracy. The theoretical angular momentum flux6 computed using density jump and corrected using equation (25) coincides with computed using entropy jump almost perfectly.

Blue dashed line in Figure 13b shows the deviation in the post-shock density at the same shock amplitude as the blue solid line (), but measured in a simulation with two times lower resolution. One can see that lowering resolution by a factor of two almost doubles the error in density jump. This exercise demonstrates that the difference between the two methods of computing theoretical angular momentum contributions due to shocks is indeed numerical.

One possible explanation for this effect might come from the fact that shocks in our simulation are not true discontinuities. Shock front is in fact smeared over several grid points, making pinpointing exact boundaries of the shock region difficult. If our shock detection algorithm (see §3.3) makes a slight spatial error in the shock position, the fluid might expand adiabatically over this small distance. This will not affect the change of entropy (explaining why the calculation based on entropy jump is robust), but will cause the relationship between the density and pressure jumps to deviate from the theoretical expectation (24). If that were the reason, then we would expect the error in jump condition to go down with resolution, and to go up with as the number of grid points needed to capture the shock increases with its amplitude. Both trends are indeed observed in our simulations.

It is also worth reiterating that predicting the shock-driven mass accretion rate based on the shock strength is in general less trivial than predicting the angular momentum deposition rate. A direct connection between and established in Rafikov (2016) is valid only in steady state disks. Their relationship becomes more complicated in disks experiencing rapid evolution, as we demonstrated in §4 and 5. Nevertheless, by carefully examining individual contributions to the angular momentum balance we can understand how the time variability of the disk properties can be accounted for (see §7) and provide an explicit prescription for computing , see equations (4) and (19).

The general one-dimensional framework for the disk evolution represented by equation (21) should be useful for other studies of shock-mediated disks (as well as for accretion disks evolving due to internal stresses). It allows for fast and accurate exploration of many evolutionary models of such disks, avoiding computationally expensive multi-dimensional simulations (which could also be contaminated by the numerical artifacts).

The disk setup used in this work (both numerical and analytical) is explicitly two-dimensional. However, we do not expect our results to change in the case of a fully three-dimensional disk, as long as one is able to adequately describe the radial and vertical behaviors for the shock strength. In that case, like in Rafikov (2016), one could still use energy and angular momentum conservation to provide a connection between the entropy production at the shock and the angular momentum deposition rate. This connection is also independent of the cooling mechanism of the disk, as passage of gas through the shock is essentially instantaneous. Cooling behavior should affect only the background temperature of the disk, which can be accounted for by explicitly considering disk thermodynamics.

Our study also neglects the effect of magnetic fields, which modify the jump conditions at the shock and should affect the calculation of . However, as long as magnetic fields remain sub-thermal (which is almost always the case for the system we focus on), their effect on the shock jump conditions and entropy calculation should be subdominant. Angular momentum transport driven by magnetized turbulence through MRI (if it operates in the disk) can be easily accounted for in the framework of our approach by simply adding corresponding stress term in equations (20) and (21).

Note that throughout our study we did not specify the driving agent for the spiral wave. In many astrophysical disks it could be caused by a massive perturber — a planet or a binary companion (Dong et al., 2015b, a; Zhu et al., 2015; Arzamasskiy et al., 2017). In that case our results would apply to the coasting phase of the wave propagation, far from the perturber, when it is no longer significantly affected by its gravity.

Alternatively, recently Montesinos et al. (2016) proposed that spiral waves may be driven by the non-uniform illumination of the protoplanetary disks by the central star, caused e.g. by the misaligned inner disk. This possibility was invoked to explain the non-axisymmetric features seen in the disks of HD 142527 (Casassus et al., 2015; Marino et al., 2015) and HD 100546 (Benisty et al., 2017). Our results would apply equally well to spiral waves sourced by this or any other mechanism (e.g. gravitational instability, see Meru et al. (2017)).

9.1. Comparison with the previous work

Several numerical studies have recently explored the effect of global spiral shocks on disk evolution. Ryan & MacFadyen (2017) directly verified the analytical predictions of Rafikov (2016) for the spiral shock-driven angular momentum transport in mini-disks around the components of the supermassive black hole binaries. There are several differences between their work and ours. First, Ryan & MacFadyen (2017) run their simulation in full general relativity, while we use Newtonian gravity. Second, to characterize the angular momentum injection by the wave they used the entropy jump at the shock, as suggested in Rafikov (2016). In our present work we use the density jump at the shock front, which is easier to determine observationally, for that purpose. Third, Ryan & MacFadyen (2017) run their simulations until the disk reached a steady state, so that the time-dependent contributions were small in their case. Our simulations are not run for nearly as long, forcing us to explicitly account for the time-dependent terms, see §7.

Ju et al. (2016, 2017) studied the role of spiral shocks in driving accretion in disks of cataclysmic variables, both in hydro and MHD setups. Our analysis of the different torque contributions is very similar to that in Ju et al. (2016). Despite the improved numerical resolution, we see our results significantly affected by numerical dissipation, especially at low amplitudes (see Figure 3b), compared to the calculations in Ju et al. (2016). Presumably, the lack of noticeable residual is their simulations is caused by the higher amplitude of the spiral shocks driven by a massive companion.

Ju et al. (2017) also found that the shock-driven transport becomes inefficient when the disk aspect ratio is reduced. This is consistent with our results in §5.2, where we show that the angular momentum deposition by the spiral shock decays with the distance traveled by the wave considerably faster in disks with lower , see Figure 6. The accelerated decay is caused by the more rapid nonlinear evolution and associated damping of the waves propagating in colder disks, see the discussion around equation (17). As a result, we naturally expect the spiral waves reaching central regions of the disks in cataclysmic variables (which have ) to have very low values of the pressure jump at the shock, rendering them inefficient at transporting angular momentum far from the disk edge. Protoplanetary (Rafikov, 2016) and circumplanetary (Zhu et al., 2016) disks, which have higher , may provide better environments for the shock-driven transport.

10. Summary

We have carried out a numerical study aimed at understanding the effects of the global spiral shocks on astrophysical disks. The main goal of this work was verification of the analytical predictions of Rafikov (2016) for the angular momentum deposition in the disk by the global shocks. Our main results can be summarized as follows.

  • Our numerical results for the angular momentum deposition rate due to disk shocks show good agreement (at the level of tens of per cent) with the predictions of Rafikov (2016), evaluated using numerically determined shock strength. This agreement persists as we vary different relevant physical parameters — initial amplitude of the density wave, aspect ratio of the disk, equation of state of the disk fluid, radial profiles of the background surface density and temperature, angular frequency of the perturbation pattern (§5).

  • High resolution is essential for fully capturing the effect of global shocks on the disk structure, and we quantify the effect of numerical dissipation on our simulations (§5.3). Properly accounting for the detrimental effects of numerical dissipation is especially important for low amplitude shocks.

  • We generalize the relationship between the shock-driven angular momentum dissipation and the mass accretion rate , derived in Rafikov (2016) for steady disks, to cover the situations when the disk evolves rapidly. In this case the time-dependent terms provide important contribution to the angular momentum budget and must be explicitly accounted for.

  • We provide a simple analytical prescription for calculating the behavior of these time-dependent terms using the knowledge of the radial pressure profile in the disk (§7). This prescription is found to be in excellent agreement with the numerical results.

  • These findings allow us to formulate a closed form equation for the shock-driven evolution of accretion disks, valid even when the disk properties change on time scales shorter than the “viscous” timescale (§8). This framework can be used for accurate and efficient modeling of the shock-driven disk evolution in one (radial) dimension.

Successful confirmation of the analytical prescription for the shock-driven disk evolution obtained in Rafikov (2016) opens a way for semi-analytical studies of the long-term disk evolution driven by the spiral shocks, once the radial profile of the shock strength (i.e. density or pressure jump at the shock) is specified. Fully self-consistent and complete description of the evolution of the density wave as it propagates through the disk (and damps), which is needed for prescribing the shock amplitude at every radius, is a subject of future work.

We thank James Stone for helpful discussions and comments, as well as for making Athena++ publicly available. Financial support for this study has been provided by NSF via grants AST-1409524, AST-1515763, and NASA via grants 14-ATP14-0059, 15-XRP15-2-0139. All simulations were carried out using facilities supported by the Princeton Institute of Computational Science and Engineering.


  1. affiliation: Department of Astrophysical Sciences, Princeton University, Ivy Lane, Princeton, NJ 08540; leva@astro.princeton.edu
  2. affiliation: Centre for Mathematical Sciences, Department of Applied Mathematics and Theoretical Physics, University of Cambridge, Wilberforce Road, Cambridge CB3 0WA, UK
  3. affiliation: Institute for Advanced Study, Einstein Drive, Princeton, NJ 08540
  4. http://princetonuniversity.github.io/athena/
  5. Although some other aspects of our calculations may be sensitive to it, see §9.
  6. We do not show it in Figure 13a to avoid confusion.


  1. ALMA Partnership, Brogan, C. L., Pérez, L. M., et al. 2015, ApJ, 808, L3
  2. Andrews, S. M., Wilner, D. J., Zhu, Z., et al. 2016, ApJ, 820, L40
  3. Arzamasskiy, L., Zhu, Z., & Stone, J. M. 2017, ArXiv e-prints, arXiv:1710.11128
  4. Bae, J., Zhu, Z., & Hartmann, L. 2017, ArXiv e-prints, arXiv:1706.03066
  5. Balbus, S. A., & Hawley, J. F. 1991, ApJ, 376, 214
  6. Balbus, S. A., & Papaloizou, J. C. B. 1999, ApJ, 521, 650
  7. Belyaev, M. A., Rafikov, R. R., & Stone, J. M. 2013, ApJ, 770, 67
  8. Benisty, M., Juhasz, A., Boccaletti, A., et al. 2015, A&A, 578, L6
  9. Benisty, M., Stolker, T., Pohl, A., et al. 2017, A&A, 597, A42
  10. Blandford, R. D., & Payne, D. G. 1982, MNRAS, 199, 883
  11. Boccaletti, A., Pantin, E., Lagrange, A.-M., et al. 2013, A&A, 560, A20
  12. Casassus, S., Wright, C. M., Marino, S., et al. 2015, ApJ, 812, 126
  13. Chandrasekhar, S. 1960, Proceedings of the National Academy of Science, 46, 253
  14. Christiaens, V., Casassus, S., Perez, S., van der Plas, G., & Ménard, F. 2014, ApJ, 785, L12
  15. Collins, K. A., Grady, C. A., Hamaguchi, K., et al. 2009, ApJ, 697, 557
  16. Dong, R., & Fung, J. 2017, ApJ, 835, 38
  17. Dong, R., Li, S., Chiang, E., & Li, H. 2017, ArXiv e-prints, arXiv:1705.04687
  18. Dong, R., Rafikov, R. R., & Stone, J. M. 2011, ApJ, 741, 57
  19. Dong, R., Zhu, Z., Fung, J., et al. 2015a, ArXiv e-prints, arXiv:1512.04949
  20. Dong, R., Zhu, Z., Rafikov, R. R., & Stone, J. M. 2015b, ApJ, 809, L5
  21. Duffell, P. C., & MacFadyen, A. I. 2012, ApJ, 755, 7
  22. Einfeldt, B. 1988, SIAM Journal on Numerical Analysis, 25, 294
  23. Fung, J., & Chiang, E. 2017, ArXiv e-prints, arXiv:1701.08161
  24. Fung, J., & Dong, R. 2015, ApJ, 815, L21
  25. Gammie, C. F. 2001, ApJ, 553, 174
  26. Gardiner, T. A., & Stone, J. M. 2005, Journal of Computational Physics, 205, 509
  27. —. 2008, Journal of Computational Physics, 227, 4123
  28. Garufi, A., Quanz, S. P., Avenhaus, H., et al. 2013, A&A, 560, A105
  29. Goodman, J., & Rafikov, R. R. 2001, ApJ, 552, 793
  30. Grady, C. A., Muto, T., Hashimoto, J., et al. 2013, ApJ, 762, 48
  31. Hennebelle, P., Lesur, G., & Fromang, S. 2016, A&A, 590, A22
  32. Jiang, Y.-F., Stone, J., & Davis, S. W. 2017, ArXiv e-prints, arXiv:1709.02845
  33. Ju, W., Stone, J. M., & Zhu, Z. 2016, ApJ, 823, 81
  34. —. 2017, ApJ, 841, 29
  35. Kim, Y., & Kim, W.-T. 2014, MNRAS, 440, 208
  36. Konigl, A. 1989, ApJ, 342, 208
  37. Landau, L. D., & Lifshitz, E. M. 1959, Fluid mechanics
  38. Larson, R. B. 1990, MNRAS, 243, 588
  39. Lee, W.-K. 2016, ApJ, 832, 166
  40. Lesur, G., Hennebelle, P., & Fromang, S. 2015, A&A, 582, L9
  41. Li, H., & Li, S. 2013, in European Physical Journal Web of Conferences, Vol. 46, European Physical Journal Web of Conferences, 05003
  42. Li, H., Lubow, S. H., Li, S., & Lin, D. N. C. 2009, ApJ, 690, L52
  43. Lynden-Bell, D., & Pringle, J. E. 1974, MNRAS, 168, 603
  44. Manser, C. J., Gänsicke, B. T., Marsh, T. R., et al. 2016, MNRAS, 455, 4467
  45. Marino, S., Perez, S., & Casassus, S. 2015, ApJ, 798, L44
  46. Marsh, T. R., & Horne, K. 1988, MNRAS, 235, 269
  47. Meru, F., Juhász, A., Ilee, J. D., et al. 2017, ApJ, 839, L24
  48. Montesinos, M., Perez, S., Casassus, S., et al. 2016, ApJ, 823, L8
  49. Muto, T., Grady, C. A., Hashimoto, J., et al. 2012, ApJ, 748, L22
  50. Nelson, R. P., Gressel, O., & Umurhan, O. M. 2013, MNRAS, 435, 2610
  51. Ogilvie, G. I., & Lubow, S. H. 2002, MNRAS, 330, 950
  52. Pérez, L. M., Carpenter, J. M., Andrews, S. M., et al. 2016, Science, 353, 1519
  53. Rafikov, R. R. 2002a, ApJ, 569, 997
  54. —. 2002b, ApJ, 572, 566
  55. —. 2015, ApJ, 804, 62
  56. —. 2016, ApJ, 831, 122
  57. Ryan, G., & MacFadyen, A. 2017, ApJ, 835, 199
  58. Savonije, G. J., Papaloizou, J. C. B., & Lin, D. N. C. 1994, MNRAS, 268, 13
  59. Sawada, K., Matsuda, T., & Hachisu, I. 1986a, MNRAS, 221, 679
  60. —. 1986b, MNRAS, 219, 75
  61. Sormani, M. C., Sobacchi, E., Shore, S. N., Treß, R. G., & Klessen, R. S. 2017, MNRAS, 471, 2932
  62. Spruit, H. C. 1987, A&A, 184, 173
  63. Spruit, H. C., Matsuda, T., Inoue, M., & Sawada, K. 1987, MNRAS, 229, 517
  64. Steeghs, D., Harlaftis, E. T., & Horne, K. 1997, MNRAS, 290, L28
  65. Stoll, M. H. R., & Kley, W. 2014, A&A, 572, A77
  66. Turner, N. J., Fromang, S., Gammie, C., et al. 2014, Protostars and Planets VI, 411
  67. Urpin, V., & Brandenburg, A. 1998, MNRAS, 294, 399
  68. Velikhov, E. P. 1959, JETP, 36, 1398
  69. Wagner, K., Apai, D., Kasper, M., & Robberto, M. 2015, ApJ, 813, L2
  70. Yu, C., Li, H., Li, S., Lubow, S. H., & Lin, D. N. C. 2010, ApJ, 712, 198
  71. Zhu, Z., Dong, R., Stone, J. M., & Rafikov, R. R. 2015, ApJ, 813, 88
  72. Zhu, Z., Ju, W., & Stone, J. M. 2016, ApJ, 832, 193
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description