Disk Accretion Driven by Spiral Shocks
Abstract
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 shockdriven 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 timedependent 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 onedimensional (in radius) evolution of the shockmediated accretion disks, which can be applied to a variety of astrophysical systems.
Subject headings:
accretion, accretion disks — protoplanetary disks — planetdisk interactions — hydrodynamics — shock waves1. 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 magnetorotational 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 metalrich white dwarf SDSS J122859.93+104032.9 (Manser et al., 2016). Recent observations of protoplanetary disks using direct imaging in optical and nearIR (Benisty et al., 2015), as well as submm 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 nonuniform 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 highamplitude spiral density waves rapidly evolve into spiral shocks, as the nonlinear 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 superEddington accretion around supermassive black holes.
On the theoretical side, properties of spiral shocks in accretion disks were first considered under the assumption of selfsimilarity, starting with the work of Spruit (1987), Spruit et al. (1987), and Larson (1990). More recently, selfsimilar 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 selfsimilar shock solutions is that they apply only to special disk models with certain scaleinvariant 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 lowmass 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 shockdriven 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 postprocessing 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 submm 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 shockdriven 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 timedependent angular momentum contribution emerging in our simulations, which allows us to formulate a general framework for following the shockdriven disk evolution in §8. We discuss and summarize our results in §9 and 10, respectively.
2. Basic equations.
We consider a twodimensional (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 selfgravity. 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)
(1) 
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
(2) 
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 nonzero , which drives their evolution and mass transport (i.e. makes nonzero) 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 longterm 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 shockmediated 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
(3) 
As a result, the mass accretion rate in a shockmediated disk can be expressed as
(4) 
This expression generalizes a similar result of Rafikov (2016), who neglected the last timedependent term in equation (4) finding (LyndenBell & Pringle, 1974)
(5) 
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
(6) 
where is the preshock value of the isothermal sound speed at the radius , and is the preshock surface density. Also, an auxiliary function of the ratio of the postshock () to preshock () pressure is given by
(7) 
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 preshock values of gas density):
(8) 
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
(9) 
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).
3. Numerical setup
We run a set of twodimensional numerical simulations in cylindrical polar coordinates using publicly available code Athena++
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 scaleheight 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, timedependent 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.
3.1. Initial conditions
We initialize the disk with the equilibrium model described in Nelson et al. (2013), which has powerlaw profiles of the surface density and temperature:
(11)  
(12) 
where and correspond to the outer radius of the simulation domain. The hydrodynamic equilibrium implies that the angular velocity of the disk
(13) 
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:
(14) 
For our fiducial run we use the uniform density profile with , as well as the globally isothermal equation of state , . The nonuniform 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 

Domain  
0.2  
Resolution 
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
(15)  
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.
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 algorithm
4. Results
We now provide a brief description of certain characteristic features of our simulations, before giving a more indepth, 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.
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.
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 nonlinearity (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 shockdriven 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 timedependent term (as it involves the timederivative 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 timedependent term. The difference in especially pronounced close to the outer disk radius, where the timedependent 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 nonzero residual (thin grey line) is explicitly defined as
(16) 
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 finiteresolution 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 timedependent 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
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 shockdriven 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 timedependent 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 timedependent 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 .
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 () shockinduced stress terms. The advective term clearly deviates from both because of the timedependent 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 shockdriven 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.
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).
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 shockdriven 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 timedependent 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 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 shockinduced angular momentum contributions exhibit a significantly steeper dependence on , see equations (7)(9) and Figure 4a.
5.4. Variation of the EOS.
Equation of state (EOS) of the disk fluid should have an important effect for shockmediated transport. This can be seen simply in the explicit dependence of the theoretical prediction (7) on the adiabatic index in the nonisothermal 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 nonisothermal runs we must use equation (7), which apparently performs as well as the equation (9) does in the isothermal case.
On the other hand, in our nonisothermal 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 timedependent contribution (see Figure 9b), which is likely caused by the continuous and spatially nonuniform injection of entropy into the disk by the shock in the nonisothermal simulations. This steady growth of entropy (absent in the isothermal case) results in a more pronounced role of the timedependent 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 powerlaw slope of the background density profile (which requires accounting for the radial pressure support in equation (13)).
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 nonisothermal 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 nonrotating pattern, , stationary in the inertial frame. In Figure 11 we show the effect of a nonzero 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 nonrotating 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.
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 timelike coordinate
(17) 
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 diskplanet 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 ringlike 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 Nwave 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 nonrotating 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 nonunique 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 semimajor 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. Timedependent contribution.
Our simulations clearly demonstrate the important role played by the timedependent 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 NavierStokes 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 )
(18) 
due to the radial pressure support in the disk.
Plugging this expression into the definition of the timedependent contribution to the angular momentum balance and azimuthally averaging, one finds
(19) 
where is the azimuthally averaged pressure.
We show how well this approximation works in Figures 4b, 6b, 7b, 9b11b. There we display (with dots) the behavior of the timedependent 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 radiallyconstant pressure, the timedependent term is generally nonzero. This happens because even the weak modification of the by the shockdriven 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 timedependent 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 timedependent 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 nonsteady disks. Another important application of this finding is discussed next.
8. Shockdriven evolution of accretion disks.
We are now in position to formulate a fully timedependent equation for surface density evolution in a shockmediated 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
(20)  
Second, equation (19) allows us to put this expression in a closed form:
(21)  
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 nonzero can be comparable to the last term (due to the timedependent 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 selfconsistent framework for following the evolution of in the shockmediated 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.
Shockmediated 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
(22) 
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 shockmediated 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 timedependent term) derived in LyndenBell & Pringle (1974), which is broadly used to study onedimensional viscous disk evolution. This is because in their derivation LyndenBell & 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 shockimposed torque, as in our case) the timedependent 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 highamplitude spiral shocks in disks of cataclysmic variables.
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 dotdashed 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
(23) 
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):
(24) 
Figure 13b shows that our simulations do not reproduce this relation exactly, i.e. the actual postshock 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
(25) 
We verified that the difference between the two methods of analytical evaluation of the shockinduced stress can indeed be approximated by equation (25) with very good accuracy. The theoretical angular momentum flux
Blue dashed line in Figure 13b shows the deviation in the postshock 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 shockdriven 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 onedimensional framework for the disk evolution represented by equation (21) should be useful for other studies of shockmediated 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 multidimensional simulations (which could also be contaminated by the numerical artifacts).
The disk setup used in this work (both numerical and analytical) is explicitly twodimensional. However, we do not expect our results to change in the case of a fully threedimensional 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 subthermal (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 nonuniform illumination of the protoplanetary disks by the central star, caused e.g. by the misaligned inner disk. This possibility was invoked to explain the nonaxisymmetric 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 shockdriven angular momentum transport in minidisks 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 timedependent contributions were small in their case. Our simulations are not run for nearly as long, forcing us to explicitly account for the timedependent 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 shockdriven 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 shockdriven 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 shockdriven 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 timedependent 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 timedependent 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 shockdriven 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 shockdriven disk evolution in one (radial) dimension.
Successful confirmation of the analytical prescription for the shockdriven disk evolution obtained in Rafikov (2016) opens a way for semianalytical studies of the longterm 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 selfconsistent 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.
Footnotes
 affiliation: Department of Astrophysical Sciences, Princeton University, Ivy Lane, Princeton, NJ 08540; leva@astro.princeton.edu
 affiliation: Centre for Mathematical Sciences, Department of Applied Mathematics and Theoretical Physics, University of Cambridge, Wilberforce Road, Cambridge CB3 0WA, UK
 affiliation: Institute for Advanced Study, Einstein Drive, Princeton, NJ 08540
 http://princetonuniversity.github.io/athena/
 Although some other aspects of our calculations may be sensitive to it, see §9.
 We do not show it in Figure 13a to avoid confusion.
References
 ALMA Partnership, Brogan, C. L., Pérez, L. M., et al. 2015, ApJ, 808, L3
 Andrews, S. M., Wilner, D. J., Zhu, Z., et al. 2016, ApJ, 820, L40
 Arzamasskiy, L., Zhu, Z., & Stone, J. M. 2017, ArXiv eprints, arXiv:1710.11128
 Bae, J., Zhu, Z., & Hartmann, L. 2017, ArXiv eprints, arXiv:1706.03066
 Balbus, S. A., & Hawley, J. F. 1991, ApJ, 376, 214
 Balbus, S. A., & Papaloizou, J. C. B. 1999, ApJ, 521, 650
 Belyaev, M. A., Rafikov, R. R., & Stone, J. M. 2013, ApJ, 770, 67
 Benisty, M., Juhasz, A., Boccaletti, A., et al. 2015, A&A, 578, L6
 Benisty, M., Stolker, T., Pohl, A., et al. 2017, A&A, 597, A42
 Blandford, R. D., & Payne, D. G. 1982, MNRAS, 199, 883
 Boccaletti, A., Pantin, E., Lagrange, A.M., et al. 2013, A&A, 560, A20
 Casassus, S., Wright, C. M., Marino, S., et al. 2015, ApJ, 812, 126
 Chandrasekhar, S. 1960, Proceedings of the National Academy of Science, 46, 253
 Christiaens, V., Casassus, S., Perez, S., van der Plas, G., & Ménard, F. 2014, ApJ, 785, L12
 Collins, K. A., Grady, C. A., Hamaguchi, K., et al. 2009, ApJ, 697, 557
 Dong, R., & Fung, J. 2017, ApJ, 835, 38
 Dong, R., Li, S., Chiang, E., & Li, H. 2017, ArXiv eprints, arXiv:1705.04687
 Dong, R., Rafikov, R. R., & Stone, J. M. 2011, ApJ, 741, 57
 Dong, R., Zhu, Z., Fung, J., et al. 2015a, ArXiv eprints, arXiv:1512.04949
 Dong, R., Zhu, Z., Rafikov, R. R., & Stone, J. M. 2015b, ApJ, 809, L5
 Duffell, P. C., & MacFadyen, A. I. 2012, ApJ, 755, 7
 Einfeldt, B. 1988, SIAM Journal on Numerical Analysis, 25, 294
 Fung, J., & Chiang, E. 2017, ArXiv eprints, arXiv:1701.08161
 Fung, J., & Dong, R. 2015, ApJ, 815, L21
 Gammie, C. F. 2001, ApJ, 553, 174
 Gardiner, T. A., & Stone, J. M. 2005, Journal of Computational Physics, 205, 509
 —. 2008, Journal of Computational Physics, 227, 4123
 Garufi, A., Quanz, S. P., Avenhaus, H., et al. 2013, A&A, 560, A105
 Goodman, J., & Rafikov, R. R. 2001, ApJ, 552, 793
 Grady, C. A., Muto, T., Hashimoto, J., et al. 2013, ApJ, 762, 48
 Hennebelle, P., Lesur, G., & Fromang, S. 2016, A&A, 590, A22
 Jiang, Y.F., Stone, J., & Davis, S. W. 2017, ArXiv eprints, arXiv:1709.02845
 Ju, W., Stone, J. M., & Zhu, Z. 2016, ApJ, 823, 81
 —. 2017, ApJ, 841, 29
 Kim, Y., & Kim, W.T. 2014, MNRAS, 440, 208
 Konigl, A. 1989, ApJ, 342, 208
 Landau, L. D., & Lifshitz, E. M. 1959, Fluid mechanics
 Larson, R. B. 1990, MNRAS, 243, 588
 Lee, W.K. 2016, ApJ, 832, 166
 Lesur, G., Hennebelle, P., & Fromang, S. 2015, A&A, 582, L9
 Li, H., & Li, S. 2013, in European Physical Journal Web of Conferences, Vol. 46, European Physical Journal Web of Conferences, 05003
 Li, H., Lubow, S. H., Li, S., & Lin, D. N. C. 2009, ApJ, 690, L52
 LyndenBell, D., & Pringle, J. E. 1974, MNRAS, 168, 603
 Manser, C. J., Gänsicke, B. T., Marsh, T. R., et al. 2016, MNRAS, 455, 4467
 Marino, S., Perez, S., & Casassus, S. 2015, ApJ, 798, L44
 Marsh, T. R., & Horne, K. 1988, MNRAS, 235, 269
 Meru, F., Juhász, A., Ilee, J. D., et al. 2017, ApJ, 839, L24
 Montesinos, M., Perez, S., Casassus, S., et al. 2016, ApJ, 823, L8
 Muto, T., Grady, C. A., Hashimoto, J., et al. 2012, ApJ, 748, L22
 Nelson, R. P., Gressel, O., & Umurhan, O. M. 2013, MNRAS, 435, 2610
 Ogilvie, G. I., & Lubow, S. H. 2002, MNRAS, 330, 950
 Pérez, L. M., Carpenter, J. M., Andrews, S. M., et al. 2016, Science, 353, 1519
 Rafikov, R. R. 2002a, ApJ, 569, 997
 —. 2002b, ApJ, 572, 566
 —. 2015, ApJ, 804, 62
 —. 2016, ApJ, 831, 122
 Ryan, G., & MacFadyen, A. 2017, ApJ, 835, 199
 Savonije, G. J., Papaloizou, J. C. B., & Lin, D. N. C. 1994, MNRAS, 268, 13
 Sawada, K., Matsuda, T., & Hachisu, I. 1986a, MNRAS, 221, 679
 —. 1986b, MNRAS, 219, 75
 Sormani, M. C., Sobacchi, E., Shore, S. N., Treß, R. G., & Klessen, R. S. 2017, MNRAS, 471, 2932
 Spruit, H. C. 1987, A&A, 184, 173
 Spruit, H. C., Matsuda, T., Inoue, M., & Sawada, K. 1987, MNRAS, 229, 517
 Steeghs, D., Harlaftis, E. T., & Horne, K. 1997, MNRAS, 290, L28
 Stoll, M. H. R., & Kley, W. 2014, A&A, 572, A77
 Turner, N. J., Fromang, S., Gammie, C., et al. 2014, Protostars and Planets VI, 411
 Urpin, V., & Brandenburg, A. 1998, MNRAS, 294, 399
 Velikhov, E. P. 1959, JETP, 36, 1398
 Wagner, K., Apai, D., Kasper, M., & Robberto, M. 2015, ApJ, 813, L2
 Yu, C., Li, H., Li, S., Lubow, S. H., & Lin, D. N. C. 2010, ApJ, 712, 198
 Zhu, Z., Dong, R., Stone, J. M., & Rafikov, R. R. 2015, ApJ, 813, 88
 Zhu, Z., Ju, W., & Stone, J. M. 2016, ApJ, 832, 193