Pileups and Migration Rates for Planets in Low Mass Disks

TL;DR

This paper studies how low-mass planets interact with viscous disks, revealing a pileup effect behind the planet that influences migration rates and can be observed in protoplanetary disks.

Contribution

It provides the first detailed 2D hydrodynamical simulations of pileup effects in low-mass disks and introduces a new Type-II migration rate model.

Findings

Disk material piles up behind planets, acting as a leaky dam.

Simulation results match theoretical predictions for moderate gap depths.

Current theory is insufficient for very deep gaps, indicating areas for further research.

Abstract

We investigate how planets interact with viscous accretion disks, in the limit that the disk is sufficiently low mass that the planet migrates more slowly than the disk material. In that case, the disk's surface density profile is determined by the disk being in viscous steady state (VSS), while overflowing the planet's orbit. We compute the VSS profiles with 2D hydrodynamical simulations, and show that disk material piles up behind the planet, with the planet effectively acting as a leaky dam. Previous 2D hydrodynamical simulations missed the pileup effect because of incorrect boundary conditions, while previous 1D models greatly overpredicted the pileup due to the neglect of non-local deposition. Our simulations quantify the magnitude of the pileup for a variety of planet masses and disk viscosities. We also calculate theoretically the magnitude of the pileup for moderately deep gaps, showing good agreement with simulations. For very deep gaps, current theory is inadequate, and we show why and what must be understood better. The pileup is important for two reasons. First, it is observable in directly imaged protoplanetary disks, and hence can be used to diagnose the mass of a planet that causes it or the viscosity within the disk. And second, it determines the planet's migration rate. Our simulations determine a new Type-II migration rate (valid for low mass disks), and show how it connects continuously with the well-verified Type-I rate.

Tables

Table 1. Table 1: Overview of all simulations used in this study, grouped by K 𝐾 K . For each simulation, the Δ ​ T Δ 𝑇 \Delta T value in parenthesis corresponds to the result of the low resolution simulation. The quantities Σ Z , ± subscript Σ 𝑍 plus-or-minus \Sigma_{Z,\pm} refer to Σ Z subscript Σ 𝑍 \Sigma_{Z} at x ± subscript 𝑥 plus-or-minus x_{\pm} .

Name	$q$	$\alpha$	$K$	$\Delta \dot{M}$	$t_{\mathrm{avg}}$	$x_{-}$	$x_{+}$	${\Sigma }_{-}$	${\Sigma }_{+}$	$\Delta T_{\mathrm{LR}}$	$\Delta T_c$	$\Delta T$
Units				$\%$	$t_{\mathrm{orb}}$	$r_p$	$r_p$	${\Sigma }_{Z,-}$	${\Sigma }_{Z,+}$	$\dot{M}{\ell }_p$	$\dot{M}{\ell }_p$	$\dot{M}{\ell }_p$
q1x3a3x4	1e-03	3e-04	10667	9.8	1000	-0.207	0.275	5.35e-02	0.719	15.8	-4.84	11.0 (14.2)
q3x4a1x4	3e-04	1e-04	2880	8.4	3000	-0.166	0.198	7.17e-02	0.277	3.83	-0.672	3.15 (4.11)
q1x3a1x3	1e-03	1e-03	3200	1.2	3000	-0.177	0.228	7.23e-02	0.276	2.58	-0.581	2.00 (2.08)
q2x3a3x3	2e-03	3e-03	4267	0.15	100	-0.197	0.260	6.73e-02	0.239	1.99	-0.733	1.25 (1.54)
q3x4a3x4	3e-04	3e-04	960	0.80	3000	-0.145	0.169	0.124	0.260	1.39	-0.205	1.19 (1.34)
q1x3a3x3	1e-03	3e-03	1067	0.66	200	-0.150	0.191	0.108	0.218	0.913	-0.324	0.590 (0.662)
q2x3a1x2	2e-03	1e-02	1280	0.61	100	-0.171	0.213	0.116	0.191	0.621	-0.308	0.313 (0.385)
q1x4a1x4	1e-04	1e-04	320	1.1	3000	-0.118	0.126	0.148	0.272	0.963	-5.76e-02	0.905 (0.830)
q3x4a1x3	3e-04	1e-03	288	0.37	100	-0.112	0.133	0.173	0.258	0.596	-0.148	0.448 (0.459)
q1x3a1x2	1e-03	1e-02	320	0.14	100	-0.107	0.119	0.170	0.166	0.486	-0.307	0.179 (0.213)
q1x4a3x4	1e-04	3e-04	107	1.8	100	-6.70e-02	8.41e-02	0.205	0.275	0.457	-3.48e-02	0.422 (0.412)
q3x4a3x3	3e-04	3e-03	96	1.1	100	-7.28e-02	7.73e-02	0.281	0.278	0.340	-0.144	0.196 (0.203)
q3x5a1x4	3e-05	1e-04	29	6.7	10000	-6.70e-02	7.06e-02	0.458	0.516	0.228	-6.07e-03	0.222 (0.239)
q1x4a1x3	1e-04	1e-03	32	0.51	100	-6.70e-02	6.39e-02	0.461	0.473	0.219	-3.76e-02	0.182 (0.188)
q3x4a1x2	3e-04	1e-02	29	0.23	100	-6.70e-02	6.39e-02	0.570	0.546	0.166	-9.95e-02	6.62e-02 (6.74e-02)
q3x5a3x4	3e-05	3e-04	9.6	1.7	100	-6.70e-02	6.39e-02	0.733	0.748	0.108	-1.23e-02	9.53e-02 (9.37e-02)
q1x4a3x3	1e-04	3e-03	11	0.63	100	-6.11e-02	5.73e-02	0.718	0.706	0.101	-3.90e-02	6.17e-02 (6.21e-02)
q1x5a1x4	1e-05	1e-04	3.2	2.9	100	-6.70e-02	6.39e-02	0.893	0.898	4.56e-02	-5.96e-03	3.96e-02 (3.87e-02)
q3x5a1x3	3e-05	1e-03	2.9	0.67	100	-6.70e-02	6.39e-02	0.904	0.902	3.75e-02	-1.22e-02	2.54e-02 (2.57e-02)
q1x4a1x2	1e-04	1e-02	3.2	0.15	100	-6.11e-02	5.73e-02	0.892	0.879	3.52e-02	-2.40e-02	1.12e-02 (1.15e-02)
q1x5a3x4	1e-05	3e-04	1.1	0.86	100	-6.70e-02	6.39e-02	0.964	0.963	1.57e-02	-4.97e-03	1.07e-02 (1.08e-02)
q3x5a3x3	3e-05	3e-03	0.96	0.31	100	-6.70e-02	6.39e-02	0.966	0.964	1.28e-02	-5.74e-03	7.01e-03 (7.37e-03)
q1x5a1x3	1e-05	1e-03	0.32	0.27	100	-6.70e-02	6.39e-02	0.990	0.989	4.50e-03	-1.88e-03	2.62e-03 (3.20e-03)
q3x5a1x2	3e-05	1e-02	0.29	0.26	100	-6.70e-02	6.39e-02	0.990	0.990	3.55e-03	-1.31e-03	2.24e-03 (2.23e-03)
q1x5a3x3	1e-05	3e-03	0.11	0.17	100	-6.70e-02	6.39e-02	0.997	0.997	1.40e-03	-4.84e-04	9.17e-04 (1.08e-03)
q1x5a1x2	1e-05	1e-02	0.03	0.24	100	-6.70e-02	6.39e-02	1.000	1.000	3.80e-04	-8.20e-05	2.99e-04 (3.02e-04)

Equations

$$\frac{\partial }{\partial t}2\pi r⟨\Sigma \ell ⟩+\frac{\partial }{\partial r}\left(2\pi r⟨\Sigma v_r\ell ⟩+F_{\nu }\right)=t_{\mathrm{ex}},$$

$$t_{\mathrm{ex}}=-2\pi r⟨{\Sigma }^{\prime }\frac{\partial {\Phi }^{\prime }}{\partial \varphi }⟩,$$

$$F_{\nu }=-2\pi r^2⟨\nu \Sigma \left(r\frac{\partial \Omega }{\partial r}+\frac{1}{r}\frac{\partial v_r}{\partial \varphi }\right)⟩.$$

$$F_{\nu }\approx -2\pi r^3\nu ⟨\Sigma ⟩{\partial }_r⟨\Omega ⟩,$$

$$F_{\nu }\approx 3\pi \nu ⟨\Sigma ⟩⟨\ell ⟩.$$

$$\Delta T={\int }_0^{\infty }{t}_{\mathrm{ex}}dr.$$

$$\frac{\partial }{\partial t}2\pi r⟨{\Sigma }^{\prime }{\ell }^{\prime }⟩+\frac{\partial }{\partial r}{F}_{\mathrm{wave}}=t_{\mathrm{ex}}-t_{\mathrm{dep}},$$

$$F_{\mathrm{wave}}$$

$$=2\pi r^2\left(⟨\Sigma ⟩⟨v_r^{\prime }{v}_{\varphi }^{\prime }⟩+⟨{\Sigma }^{\prime }{v}_r^{\prime }{v}_{\varphi }^{\prime }⟩+⟨v_r⟩⟨{\Sigma }^{\prime }{v}_{\varphi }^{\prime }⟩\right).$$

$$\frac{d}{dr}{F}_{\mathrm{wave}}\approx t_{\mathrm{ex}}-t_{\mathrm{dep}}(\mathrm{WSS}).$$

$$\frac{\partial }{\partial t}2\pi r⟨\Sigma ⟩⟨\ell ⟩+\frac{\partial }{\partial r}\left(-\dot{M}⟨\ell ⟩+F_{\nu }\right)=t_{\mathrm{dep}},$$

$$\dot{M}=-2\pi r⟨\Sigma v_r⟩,$$

$$\frac{\partial }{\partial t}2\pi r⟨\Sigma ⟩=\frac{\partial }{\partial r}\dot{M},$$

$$\dot{M}$$

$$\frac{d}{dr}\left(-\dot{M}⟨\ell ⟩+F_{\nu }\right)$$

$$F_{\nu }(r)=\dot{M}⟨\ell ⟩+F_{\ast }+{\int }_{r_i}^r{t}_{dep}(r^{\prime })d{r}^{\prime }.$$

$$F_{\nu }(r)=\dot{M}⟨\ell ⟩+{\int }_{r_i}^r{t}_{dep}(r^{\prime })d{r}^{\prime }(\mathrm{VSS}).$$

$$F_{\nu }=\dot{M}\ell \times \{\begin{array}{ll}1& r<r_i\\ 1+\frac{\Delta T}{\dot{M}\ell }& r>r_o\end{array}$$

$$\Sigma$$

$${\Sigma }_Z\equiv \frac{\dot{M}}{3\pi \nu },$$

$$T_{\pm }$$

$$C_{+}$$

$$C_{-}$$

$${\left(F_{\nu }-\dot{M}\ell \right)}_{r_p}$$

$$=$$

$$\frac{{\Sigma }_p}{{\Sigma }_{Z,p}}\approx \frac{1}{1+|C_{-}|K/(3\pi )}\approx \frac{1}{1+0.04K},$$

$$K\equiv \frac{q^2}{\alpha h^5},$$

$$\Delta T$$

$$\approx$$

$$t_{\mathrm{ex}}=\pm C\Sigma \frac{q^2}{x^4}\left(1+2.26x\right),$$

$$T_{\pm }\approx \pm C{\Sigma }_{\pm }\frac{q^2}{|x_{\pm }{|}^3}\left(1+2.26x_{\pm }\right),$$

$$\Delta T=4.3_{2.8}^{6.6}{q}^{1.05\pm 0.06}{\alpha }^{-0.91\pm 0.04}\dot{M}{\ell }_p,$$

$$\Sigma \approx {\Sigma }_Z\left(1+\frac{\Delta T}{\dot{M}\ell }\right).$$

Full text

planet-disk interactions, protoplanetary disks, accretion disks

1 Introduction

Protoplanetary disks are being observed in ever-increasing detail, e.g., via imaging and spectral studies (Williams & Cieza, 2011; Espaillat et al., 2014; ALMA Partnership et al., 2015; Andrews et al., 2016, 2018). The inferred properties of these disks can be used to test theories for protoplanetary disk evolution and planet formation. For example, many imaged disks exhibit bands and gaps that may be sculpted by planetary mass companions (e.g., ALMA Partnership et al., 2015; Kanagawa et al., 2015a; de Boer et al., 2016; Isella et al., 2016; Fedele et al., 2017; Dong et al., 2018; Zhang et al., 2018, and others). In addition, SEDs (and images) of so-called transitional disks reveal that they have inner holes, which might be emptied out by planets (Zhu et al., 2011; Andrews et al., 2011; Dodson-Robinson & Salyk, 2011; Zhu et al., 2012; Espaillat et al., 2014; van der Marel et al., 2016; Haffert et al., 2019).

The planet-disk interaction problem has been studied extensively (for reviews see e.g., Lin & Papaloizou, 1993; Kley & Nelson, 2012; Baruteau et al., 2014). Planets torque material in the disk by launching spiral waves, which then damp and deposit their angular momentum in the disk. If the planet is sufficiently massive, the torques open up a gap in the disk (Lin & Papaloizou, 1986a). An inevitable corollary of the planet torquing the disk is the disk torquing the planet, and the resulting migration of the planet. Historically, the migration timescales were thought to fall into two main categories: Type I for low mass planets that are unable to open gaps in their disks; and Type II for high mass planets that open nearly infinitely deep gaps, locking them into the disk’s viscously driven accretion (Ward, 1997).

Simplified 1D models have been constructed for the mutual evolution of the planet and disk when the gap is very deep (e.g., Syer & Clarke, 1995; Ivanov et al., 1999; Ward, 1997; Liu & Shapiro, 2010; Kocsis et al., 2012a, b). But, as we shall show in this paper, such models have a serious difficulty due to the fact that they usually assume (either explicitly or implicitly) “local deposition”: i.e., that waves damp immediately after being launched. In reality, however, deposition is non-local, as waves transport angular momentum from where they are excited to where they are damped. Under local deposition, gaps become extremely deep—with the depth depending exponentially on the planet mass (e.g., Tanigawa & Ikoma, 2007; Liu & Shapiro, 2010; Fung et al., 2014; Kanagawa et al., 2015b). In early work (Syer & Clarke, 1995; Ivanov et al., 1999; Ward, 1997), it was assumed that massive enough planets opened infinitely deep gaps, and thus no material could flow across the planet’s orbit. In Kocsis et al. (2012a, b), exponentially deep gaps were considered. They made the same low-disk-mass approximation that we make in this paper: that the disk remains in viscous steady state as the planet migrates. However, their work is based on the local approximation. We discuss their work further in §5.5. A key result of many of the aforementioned 1D models (e.g., Syer & Clarke, 1995; Liu & Shapiro, 2010; Kocsis et al., 2012a, b) is that they produce an enhancement of gas exterior to the planet’s orbit (termed a pileup), and this pileup follows the planet as it migrates inwards.

More recently, there has been a significant amount of progress in understanding deep gaps around large planets. Based initially on hydrodynamical simulations, several authors (e.g., Crida et al., 2006; Fung et al., 2014) found that even very massive planets do not open exponentially clean gaps. Instead, their results follow a (non-exponential) scaling relationship which can be derived analytically if one assumes most of the angular momentum injected by the planet comes from nearby the planet (Kanagawa et al., 2015b; Duffell, 2015). Dong & Fung (2017), Kanagawa et al. (2017), and Duffell (2019) expanded the parameter space covered by these simulations and have found similar scaling relations – even showing that gaps in 3D are similar to gaps in 2D (Fung & Chiang, 2016). However, none of the hydrodynamical studies have reproduced the pileup effect seen in the local deposition 1D models. That seems somewhat puzzling as one would expect that a very massive planet should act as a barrier to accreting material, and even if the barrier is partial, it should slow down the accretion of the gas. This might be important observationally, as the pileup may be responsible for the inner hole in transitional disks (or, more precisely, the pileup is the observed part of the disk with the inner hole).

A second outstanding question for planets that open deep gaps—in addition to the existence of a pileup—is what is their migration rate? This has been addressed recently by several authors. Duffell et al. (2014) and Dürmann & Kley (2015) found that gap opening planets are not locked into the disk’s viscous evolution, but instead migrate at a range of rates set primarily by the disk-to-planet mass ratio, with larger disk masses resulting in faster migration. Kanagawa et al. (2018) found similar results and provided an empirical formula for the migration rate which smoothly connects the non-gap opening regime to the deep gap regime. In slight tension with these results, however, the simulations of Robert et al. (2018) showed that the deep gap migration rates, while not being exactly the Type II rate, were still proportional to the disk’s viscosity.

In the present work, we address both of these questions, focusing on a particularly simple case: when the disk is sufficiently low-mass that the planet migration rate is slower than the disk material. As we shall see, this results in a particularly clean setup, since the disk’s viscous steady state structure can be studied while ignoring planet migration.

The outline of the paper is as follows. In §2 we set up the planet-disk interaction problem in low mass disks, and study it analytically to the extent that we can. Our main result is that there is a single quantity that needs to be determined: the total amount of angular momentum put into the disk by the planet ($\Delta T$) when the disk is in viscous steady state. This quantity controls both the pileup and the planet’s migration rate. In §3–4 we turn to hydrodynamical simulations with the primary goal of determining $\Delta T$: in §3 we outline our numerical method, focusing on our new boundary conditions which allow the disk to settle into the steady-state solution described in §2; and in §4 we present the results of a suite hydrodynamical simulations. In §5, we consider some implications of our simulations, including the planet migration rate. Finally, we summarize and list some open questions in §6 and §7.

2 Planets in Low-mass disks

Our basic assumption throughout this paper is that the accretion disk is sufficiently low mass that the planet migrates more slowly than the disk material. As we show in §5.1, for the parameter-range that we consider, a disk qualifies as low-mass if it is slightly less massive than the planet; in fact, in some cases the disk can even be more massive than the planet and still qualify as low-mass. Such disks are relevant for Jupiter-mass planets, and may also be relevant for terrestrial planets during the late stages of planet formation. We shall also assume that the planet is circular with radius $r_{p}$ and mean motion $\Omega_{p}$.

Our primary goal is to calculate the planet-disk torque ($\Delta T)$111We denote it $\Delta T$ because it is a sum of positive torque on exterior material and negative torque on interior material. once the disk has reached viscous steady state (VSS). This torque is important for two reasons. First, it affects the surface density profile of the disk in VSS, leading to a pileup of material outside of the planet’s orbit. Second, it forces the planet to migrate by removing the planet’s angular momentum. As the planet migrates, the disk passes through successive VSS solutions. Therefore in the aforementioned low-disk-mass limit, we may obtain the planet’s migration rate by considering the dynamics on timescales long enough for the disk to reach VSS, while neglecting dynamics on the longer migration time.

Three dimensionless parameters affect $\Delta T$ in a non-trivial way: the planet-star mass ratio ($q$), the strength of viscosity (e.g., as parameterized by the Shakura & Sunyaev $\alpha$), and the aspect ratio of the disk ($h\equiv H/R$). A fourth potential parameter is the accretion rate of the disk in VSS ($\dot{M}$), or equivalently the overall amplitude of the surface density. But with the fairly standard assumptions that we shall make, $\Delta T\propto\dot{M}$ 222In particular, the proportionality $\Delta T\propto\dot{M}$ relies on the assumption that the disk is locally isothermal. We suspect that using a more realistic equation of state will not change our results significantly, but leave the verification to future work. , i.e., the dependence on $\dot{M}$ is trivial. As a result, we shall calculate the dimensionless torque $\Delta T/(\dot{M}\ell_{p})$ where $\ell_{p}=r_{p}^{2}\Omega_{p}$ is the specific angular momentum of the planet, and this will be a function of three dimensionless parameters ($q,\alpha,h$).

2.1 Excitation and Deposition of Angular Momentum

With the planet’s orbit fixed and circular, there are two timescales on which the disk’s properties evolve—the wave and viscous timescales. On the faster wave timescale, the planet excites waves in the disk, and these propagate away from the planet where they damp by viscosity or shocks. On this timescale wave steady state (WSS) is reached, meaning that the wave pattern becomes stationary in the rotating reference frame of the planet. On the viscous timescale, the azimuthally-averaged (“mean”) surface density reacts to the damping of the waves, and VSS is reached. The wave timescale is $\sim{\rm(orbital\ time)}/h$, because the group velocity of pressure waves is of order the sound speed (Ogilvie & Lubow, 2002). The viscous timescale is $\sim{\rm(orbital\ time)}/(\alpha h^{2})$, and therefore significantly longer.

Angular momentum is transferred from the planet to the disk in a two-stage process: (i) Excitation: the planet excites waves at Lindblad resonances, where it transfers angular momentum to the waves; and (ii) Deposition: after the waves propagate, they deposit their angular momentum in the disk.333 For clarity, we ignore here the complication that angular momentum is also transferred to co-orbital material, which does not launch propagating waves. We show below that co-orbital torques are typically sub-dominant.

The distinction between excitation and deposition is sometimes ignored in the literature (e.g., Ward, 1997; Liu & Shapiro, 2010; Kocsis et al., 2012a, b, although the former reference discusses some of the effects of this distinction). Nonetheless, it is of crucial importance (e.g., Lunine & Stevenson, 1982; Greenberg, 1983; Goodman & Rafikov, 2001; Rafikov, 2002a, b; Muto et al., 2010; Duffell, 2015; Kanagawa et al., 2015b, 2017; Ginzburg & Sari, 2018). In particular, whereas excitation is straightforward to calculate from linear theory (Goldreich & Tremaine, 1980), it is deposition that controls the surface density profile of the disk.

We make the above discussion quantitative via equations that track angular momentum transfer. The reader uninterested in technical details may skip the remainder of this subsection without great loss. We consider a 2D disk, in which the dynamical variables are $\{\Sigma,v_{r},v_{\phi}\}$ in standard notation. Where convenient, we shall also employ $\ell\equiv rv_{\phi}$ (specific angular momentum) and $\Omega\equiv v_{\phi}/r$ as surrogates for $v_{\phi}$. Variables are decomposed into “mean” and “wave” components, denoted by brackets and primes respectively, e.g., $\Sigma=\langle\Sigma\rangle+\Sigma^{\prime}$, where $\langle\Sigma\rangle\equiv\oint\Sigma d\phi/(2\pi)$.

In Appendix A, we start from the general 2D equations of motion for a disk with shear viscosity (Eqs. (A1)-(A2)) to derive exact equations for three angular momentum densities: the total ($\propto\left\langle\Sigma\ell\right\rangle$), wave ($\propto\left\langle\Sigma^{\prime}\ell^{\prime}\right\rangle$), and mean flow ($\propto\left\langle\Sigma\right\rangle\left\langle\ell\right\rangle$), with the following results.

Total:

[TABLE]

where the excitation torque density444 We adopt the convention of representing torque densities (i.e., torque per unit radius) by lower case $t$, and torques (or angular momentum fluxes) by capitalized $F$ or $T$. is

[TABLE]

in which $\Phi^{\prime}$ is the wave component of the planet’s potential, and the viscous torque is

[TABLE]

An approximate form for $t_{\rm ex}$ is derived by Goldreich & Tremaine (1980), resulting in the well-known “standard torque formula” ($t_{\rm ex}\propto\Sigma q^{2}/|r-r_{p}|^{4}$). We shall show from our simulations that, with minor modifications, the standard torque formula is of satisfactory accuracy, even when the waves are nonlinear.

For $F_{\nu}$, we may typically neglect its non-wave contribution to approximate

[TABLE]

as we shall verify in our hydrodynamical simulations (see Appendix A.1 and also Kanagawa et al. (2017)). A further approximation is to assume that $\langle\Omega\rangle$ is nearly Keplerian, as is typically true everywhere except near the bottom of a deep gap, resulting in the familiar form

[TABLE]

Equation (1) shows that the disk locally conserves angular momentum, aside from that input by the planet ($t_{\rm ex}$). The total torque that the planet applies to the disk is the quantity that we ultimately desire:

[TABLE] 2. 2.

Wave:

[TABLE]

where the flux of angular momentum carried by the waves is

[TABLE]

In the latter expression, the first term $\propto\langle v_{r}^{\prime}v_{\phi}^{\prime}\rangle$ is the usual wave flux that is conserved in the linearized adiabatic problem without viscosity (e.g., Goldreich & Tremaine, 1979); the second is a triple correlation that can become of comparable importance when the waves are nearly nonlinear; and the third is generally negligible outside of the co-orbital zone because $\left\langle v_{r}\right\rangle$ is ${\cal O}(\alpha h^{2})$ (for an example see Appendix A.1).

The quantity $t_{\rm dep}$ in Eq. (7) is the deposition torque density; $t_{\rm dep}$ is displayed explicitly in Eq. (A9), but for present purposes it suffices to note that it vanishes wherever there are no waves.555Our expression for $t_{\rm dep}$ does not vanish when the viscosity vanishes. That is a consequence of using a locally isothermal equation of state. If the more realistic locally adiabatic equation of state were used, our expression for $t_{\rm dep}$ would vanish at zero viscosity (Miranda & Rafikov, 2019). Of course, the case of zero viscosity is not physical: at small viscosity, dissipation occurs at shocks (Goodman & Rafikov, 2001).

The wave equation may typically be simplified: since the wave timescale is shorter than the viscous one, when considering viscous evolution one may drop the $\partial/\partial t$ in that equation, yielding the wave steady state (WSS) equation

[TABLE]

To appreciate the implication, let us focus for definiteness on orbital radii $r>r_{\rm p}$, in which case the outer torque excited by the planet is $T_{+}=\int_{r_{p}}^{\infty}t_{\rm ex}dr=\int_{r_{p}}^{\infty}t_{\rm dep}dr$, where the latter relation follows from the fact that the wave flux vanishes at infinity and at the planet (ignoring co-orbital torques). In other words, the total exterior torque excited by the planet is equal to that deposited into the mean flow, with $F_{\rm wave}$ the intermediary that transports angular momentum from where it is excited (Lindblad resonances) to where it is deposited. 3. 3.

Mean flow:

[TABLE]

where

[TABLE]

is the mass accretion rate666Our sign convention is such that $\dot{M}>0$ corresponds to inwards mass accretion, i.e. $v_{r}<0$..

Equations (1)–(12) follow from the general 2D equations without approximation (aside from those denoted explicitly with $\approx$ that we have made for simplicity)—in particular, they do not assume that wave quantities are smaller than mean ones, or a specific form for the shear viscosity or equation of state. Note that Eq. (1) is equal to the sum of Eqs. (7) and (11) because $\langle\Sigma\ell\rangle=\langle\Sigma^{\prime}\ell^{\prime}\rangle+\langle\Sigma\rangle\langle\ell\rangle$.

Kanagawa et al. (2017) perform a similar decomposition to that presented above, but for the steady-state equations. As a result, they do not distinguish between WSS and VSS. Equation (11) is well-known, e.g., Ward (1997) and Rafikov (2002b). But the decomposition into waves and mean flow allows for an explicit general expression for $t_{\rm dep}$, given in Eq. (A9).

2.2 Viscous Evolution

The mean flow equation (Eq. 11) governs the evolution of the mean surface density on the viscous timescale. It is equivalent to the standard angular momentum equation for a planet-less viscous accretion disk (e.g., Lynden-Bell & Pringle, 1974), aside from the term $t_{\rm dep}$, which is caused by the transfer of angular momentum from waves to the mean flow as they damp. The fluxes in this equation are the viscous $F_{\nu}$, as in Eq. (1), and $-\dot{M}\left\langle\ell\right\rangle$, which represents the inward advective transport of the mean flow’s $\left\langle\ell\right\rangle$777 Note that $\dot{M}\propto\left\langle\Sigma v_{r}\right\rangle=\left\langle\Sigma\right\rangle\left\langle v_{r}\right\rangle+\left\langle\Sigma^{\prime}v_{r}^{\prime}\right\rangle$, and therefore the waves participate in the advection. If one wished, the $\left\langle\Sigma^{\prime}v_{r}^{\prime}\right\rangle$ term could be transferred into the definition of $t_{\rm dep}$, while at the same time adding it into $F_{\rm wave}$. We choose not to do so because we wish the viscous evolution equations (Eqs. (11) and (13)) to form a closed set for the quantities $\dot{M}$ and $\left\langle\Sigma\right\rangle$, once $t_{\rm dep}$ is known..

As the disk evolves towards VSS, its viscous evolution is determined by Eq. (11), together with mass conservation:

[TABLE]

where $\dot{M}$ is defined in Eq. (12). Equations (11) and (13) form a closed set of equations for $\left\langle\Sigma\right\rangle$ and $\dot{M}$ because (i) $\left\langle\ell\right\rangle$ is nearly Keplerian, aside from a (typically small) correction which may be determined from $\left\langle\Sigma\right\rangle$ by radial pressure balance; (ii) $F_{\nu}$ may be approximated by Eq. 4; and (iii) the $t_{\rm dep}$ profile can be calculated from the $\langle\Sigma\rangle$ profile, because the waves may be considered to be in WSS. However, determining the $t_{\rm dep}$ profile theoretically is difficult, even under the WSS assumption, and we shall resort to numerical simulations to determine it (§3-4)888Goodman & Rafikov (2001), and Rafikov (2002a) calculate $t_{\rm dep}$ for sub-thermal-mass planets ($q\lesssim h^{3}$) and in the absence of viscosity; see also Duffell (2015) and Ginzburg & Sari (2018) who use a somewhat crude approximation to $t_{\rm dep}$ in order to determine self-consistent gap profiles. But those results are not directly applicable to the higher-mass planets and viscous disks that we consider in this paper..

2.3 Viscous Steady State (VSS)

In VSS, Eqs. (11) and (13) imply

[TABLE]

For a given $t_{\rm dep}$ profile, the solution of the second equation is trivially

[TABLE]

Here, $F_{*}$ is an integration constant and $r_{i}$ is arbitrary, but we shall choose it to be sufficiently inwards of the planet’s orbit that the waves have damped by then, and so $t_{dep}$ vanishes at $r<r_{i}$999We ignore here the possibility that waves can sometimes reach the inner edge of the disk if the viscosity is small enough (Rafikov, 2002a).. The constant $F_{*}$ represents the angular momentum flux injected at the inner edge of the disk, e.g., by the star. It is often chosen to yield $F_{\nu}=\Sigma=0$ near the star’s surface (Shakura & Sunyaev, 1973; Lynden-Bell & Pringle, 1974). But because the first term in Eq. (16) increases as $\left\langle\ell\right\rangle\propto r^{1/2}$, we may discard the constant $F_{*}$ more than a few stellar radii away from the star. The VSS solution may therefore be written as

[TABLE]

provided $r$ is far enough from the star (as we assume to be true henceforth). We shall make use of this VSS solution extensively in our analysis below. Note that the profile of $F_{\nu}(r)$ immediately determines the $\left\langle\Sigma\right\rangle$ profile after inserting an approximate form for $F_{\nu}$ (Eq. (4) or (5)). In other words, it is $t_{\rm dep}$ (rather than $t_{\rm ex}$) that directly controls the surface density profile of the disk.

2.3.1 VSS solution far from the planet: connecting $\Delta T$ to the pileup and the migration rate

We may apply the above solution to determine the density profile far from the planet:

[TABLE]

where $r_{o}$ is the distance beyond which $t_{dep}\approx 0$ and $\Delta T$ is the total deposited torque, which must be equal to the excited torque (Eq. 10). We have dropped angled brackets, because the waves are damped in this domain. Using $F_{\nu}\approx 3\pi\nu\Sigma\ell$ (Eq. 5), the VSS surface density is

[TABLE]

The quantity

[TABLE]

that appears in Eq. (19) is the well-known solution for a planet-less accretion disk far from the star (Lynden-Bell & Pringle, 1974)—or equivalently one with $F_{*}=0$. For ease of reference below, we call it the “zero-angular-momentum-flux” (or ZAM) solution, whence the subscript Z.

Steady-state solutions given by Eqs. (18) and (19) with $\Delta T\neq 0$ have previously been studied in the context of the disk inner boundary where $\Delta T$ is the torque of the central star on the disk (e.g., Shakura & Sunyaev, 1973; Lynden-Bell & Pringle, 1974), and in the context of circumbinary disks where $\Delta T$ is the total torque of the binary on the disk (e.g., Syer & Clarke, 1995; Kocsis et al., 2012a, b; Rafikov, 2013, 2016; Miranda et al., 2017; Tang et al., 2017; Muñoz et al., 2019).

To better understand the VSS solution described by Eqs. (18) and (19) we plot an illustrative example in Figure 1. The curves are taken from one of our hydrodynamical simulations described in §4. The top panel shows that far inside of the planet, $\Sigma=\Sigma_{Z}$, while far outside there is a pileup relative to $\Sigma_{Z}$. The bottom panel shows that far outside the planet $F_{\nu}-\dot{M}\ell=\Delta T$, which is spatially constant; the value of $\Delta T$ determines the height of the pileup in $\Sigma$.

The torque deposited into the disk, $\Delta T$, comes at the expense of the planet’s orbit, implying that the planet’s instantaneous migration rate is $\dot{r}_{p}=-2r_{p}\Delta T/(M_{p}\ell_{p})$ where $M_{p}$ is the planet mass and where we assume that the planet maintains a circular orbit and does not accrete any material. A consequence of the above is that whenever $\Delta T>0$ the planet will migrate inwards, and will be accompanied by a pileup outside of its orbit.

2.3.2 Calculation of $\Delta T$ in moderately deep gaps, and the difficulty

with very deep gaps

We apply here the VSS equation to determine $\Delta T$ for gaps that are “moderately deep” (to be defined shortly). The results will be shown to match those from simulations for gaps that are $\gtrsim 25\%$ of the background density. Although moderately deep gaps are only of moderate interest—particularly if one is interested in large pileups—we present the theory here because it helps clarify the results of the simulations to be presented shortly, and it also shows why the theory is much more difficult for deeper gaps. The theory for moderately deep gaps was developed by Duffell (2015) and Kanagawa et al. (2015b). We mostly follow their approach, but extend it to calculate the two-sided torque ($\Delta T$).

The standard torque formula is $t_{\rm ex}\propto\Sigma q^{2}/|r-r_{p}|^{4}$, with a cutoff at $|r-r_{p}|\lesssim h$ (Goldreich & Tremaine, 1980)101010We ignore co-orbital torques in this section, but consider their impact in §4.5 in the context of our simulations. . Therefore, provided that the gap is not too deep, most of the excited torque comes from a distance $\sim h$ from the planet. For such “moderately deep” gaps, one may set the inner excited torque to $T_{-}\equiv\int_{0}^{r_{p}}t_{\rm ex}dr={\rm const}\times q^{2}\Sigma_{p}/h^{3}$, where $\Sigma_{p}$ is the value at the planet.111111We implicitly assume that $\Sigma$ does not vary significantly between $r_{p}$ and $r_{p}-h$, which is expected to be true even when there is a gap, because the lengthscale of the gap is set by $t_{\rm dep}$ rather than $t_{\rm ex}$. Our simulation results support this expectation (see Figure 5 below). A similar argument applies to the outer excited torque ($T_{+}$), but with a different constant. In order to obtain these constants, one must account for the detailed shape of $t_{\rm ex}$ near the torque cutoff, which we do by numerically solving the linear equations of motion; details are in Appendix B.121212To isolate the one-sided Lindblad torques we make use of the property that the linear wave flux far from the planet carries all of the Lindblad torque. We solve the linear equations on top of a background $\Sigma=\Sigma_{Z}$, with a small enough viscosity such that we can measure the wave angular momentum flux far from the planet (Korycansky & Pollack, 1993). Since the wave flux is not conserved in the locally isothermal equations (see e.g., Lee, 2016; Miranda & Rafikov, 2019), we also assume $c_{s}=h={\rm const}$ for the linear calculation. The effect of the $c_{s}$ gradient should be sub-dominant as the important Lindblad resonances are located $\sim h$ away from the planet. We find

[TABLE]

where these $C_{\pm}$ are applicable for $h=0.05$ and $\Sigma_{Z}\propto r^{-1/2}$, which are the values we shall use in our simulations. More general expressions can be found in, e.g., Tanaka et al. (2002).

We may now obtain $\Sigma_{p}$ from the VSS equation (Eq. 17) at $r_{p}$:

[TABLE]

where the second equality follows from the excited torque being equal to the deposited torque (Eq. 10). Setting $F_{\nu}=3\pi\nu\Sigma\ell$, $\dot{M}=3\pi\nu\Sigma_{Z}$, and $\nu=\alpha h^{2}\ell$ we obtain for the depth of the gap

[TABLE]

where

[TABLE]

is a commonly used parameter that measures the relative strength of a planet’s gravitational torque ($\propto q^{2}/h^{3}$ at distance $h$) to the disk’s viscous torque ($\propto\alpha h^{2}$ at distance $h$) (Ward, 1997; Duffell, 2015; Kanagawa et al., 2015b, 2017). Inserting $\Sigma_{p}$ into Eq. (21) yields the one-sided torques $T_{+}$ and $T_{-}$. The two-sided torque $\Delta T=T_{+}+T_{-}$ is then

[TABLE]

We see that the reason for the existence of a non-vanishing $\Delta T$ is that the exterior torques exceed the interior torques ($|C_{+}|>|C_{-}|$) by ${\cal O}(h)$, as is well-known from studies of Type I migration (Goldreich & Tremaine, 1980; Ward, 1997). Duffell (2015) and Kanagawa et al. (2015b) previously derived Eq. (26), while Kanagawa et al. (2018) derived Eq. (29) with the numerical coefficient extracted from their simulations.

We show below that this first-principles prediction for gap depth and $\Delta T$ agrees well with simulation results for moderately deep gaps.131313We have assumed in our derivation that the torque cutoff occurs at $h$, which is true for sub-thermal-mass planets. We have not explored the case of super-thermal-mass planets because in our simulations the planets that produce moderately deep gaps are mostly sub-thermal. At first glance, it might appear surprising that one may predict the gap depth and torques without any knowledge of the $t_{\rm dep}$ profile. The reason is that once one knows where the torques are excited, one may calculate the ratio of excited torque to surface density at that location. Since the total excited torque is equal to the total deposited torque, and since the deposited torque determines the surface density, one then has a closed system of equations. We may illustrate this by equating the timescale for a planet to open a gap with that required for viscosity to close it, as must be true in VSS. The latter time is $t_{\rm close}\sim x_{\rm gap}^{2}/\nu$, where $x_{\rm gap}$ is the width of the gap, which is set by the width of the $t_{\rm dep}$ profile. The former is $t_{\rm open}\sim L/T_{-}$, where $L$ is the angular momentum required to vacate material from the gap, $L\sim(\Sigma_{Z,p}-\Sigma_{p})x_{\rm gap}^{2}$ (setting $r_{p}=\Omega_{p}=1$ here for simplicity). Equating the two timescales, we see that $x_{\rm gap}^{2}$ cancels, and with $T_{-}\sim\Sigma_{p}q^{2}/h^{3}$ and $\nu\sim\alpha h^{2}$, we find $\Sigma_{p}/\Sigma_{Z,p}\sim 1/(1+K)$, in agreement with the form of Eq. (26); a slightly more careful calculation can nearly reproduce the order-unity coefficient multiplying $K$.

The above discussion suggests that constructing a theory for very deep gaps will be difficult. Once a gap is sufficiently deep, most of the torque will be excited beyond $h$ from the planet (Petrovich & Rafikov, 2012; Ginzburg & Sari, 2018). In order to determine that distance, one needs to know the amount of torque deposited inside of that distance, since that will affect the surface density there, which in turn controls where the torque is excited. But, as previously emphasized, understanding deposition is difficult. An additional, though related, difficulty is that in our derivation for moderately deep gaps we assumed that $\Sigma$ was nearly constant between the inner and outer excitation locations ($r_{p}-h$ and $r_{p}+h$). But for very deep gaps that is no longer true, and the jump in $\Sigma$ from its inner to its outer excitation location is also determined by the deposition of torque within that zone. Further discussion of very deep gaps, in light of our simulation results, will be presented in §5.4.

3 Numerical Method

Our main goal in the remainder of the paper is to calculate $\Delta T$ with hydrodynamical simulations, and to understand and explain the results theoretically. Henceforth, we shall set the planet’s orbital radius and mean motion to $r_{p}=1$ and $\Omega_{p}=1$, which sets the length and time units. Note that one could also set $\dot{M}$ to unity, because the viscous evolution equations are linear in $\dot{M}$ for the locally isothermal equation of state that we adopt. However, we prefer to keep dependences on $\dot{M}$ explicit to avoid potential confusion (where helpful, we also keep some dependences on $r_{p}$ and $\Omega_{p}$ explicit).

Our numerical setup is mostly standard, with one main exception: the inner and outer boundary conditions are based on the VSS solution far from the planet (Eq. 19), which allows us to find a pileup where others have not.

We use the GPU-accelerated FARGO3D code (Benítez-Llambay & Masset, 2016) to evolve the 2D hydrodynamical equations of motion, Eqs. (A1)-(A2). We take the equation of state to be locally isothermal, $P=c_{s}^{2}\Sigma=h^{2}(r\Omega_{K})^{2}\Sigma$, where $\Omega_{K}$ is the Keplerian orbital frequency and $h$ is the aspect ratio, which we fix at $h=0.05$. Viscosity is modelled explicitly, with kinematic shear viscosity $\nu=\alpha c_{s}^{2}/\Omega_{K}=\alpha h^{2}r^{2}\Omega_{K}$ and constant $\alpha$.

The planet is modeled as a softened gravitational potential with softening length $\epsilon=0.6h=0.03$. This value approximates the vertically-averaged 3D potential to within $\sim 10\%$ for distances $|r-1|>h$ (Müller et al., 2012). We present further details of our numerical setup in Appendix C.

3.1 Boundary Conditions

Boundary conditions are implemented in FARGO3D with layers of ghost cells interior to our inner boundary, located at $r_{i}$, and exterior to our outer boundary, located at $r_{o}$. At the outer boundary, we wish to supply the system with a steady mass accretion rate $\dot{M}$ without injecting any angular momentum at the inner boundary, i.e., we want $F_{*}=0$ in Eq. (17). At the inner boundary the solution should match onto the ZAM solution (Eq. 20), where the $\dot{M}$ that appears in the solution should be $\dot{M}|_{r_{i}}$ (rather than the injected $\dot{M}|_{r_{o}}$). To account for this, we match onto the ZAM solution by ensuring that $\nu\Sigma$ is constant across the boundary, i.e., we set the value of $\Sigma$ in ghost cells at $r<r_{i}$ such that $\nu\Sigma$ in those ghost cells is constant and equal to the value in the cell at $r=r_{i}$. For $v_{r}$ in the inner ghost cells, we set it to the ZAM expression $v_{r}=3\nu/(2r)$, i.e. the value consistent with $\Sigma$ and $\dot{M}|_{r_{i}}$. Finally, for $v_{\phi}$, we set it to its (pressure-supported) Keplerian value by extrapolating from the cell at $r=r_{i}$. We also ensure that the flow is axisymmetric at $r_{i}$ by adopting a wave-killing zone between $r_{i}<r<r_{i,wkz}$, where we damp $v_{r}$ to its azimuthal average (see Eq. C1). By damping only $v_{r}$ and not $v_{\phi}$ or $\Sigma$, we ensure that the wave-killing procedure conserves angular momentum (and mass), and therefore all of the angular momentum excited by the planet—and carried by density waves—is deposited into the disk within the computational domain.

At the outer boundary, we seek to match onto the exterior VSS solution $\Sigma=\Sigma_{Z}(1+\Delta T/(\dot{M}\ell))$ (Eq. 19), where $\dot{M}$ is the mass to be injected at $r_{o}$, and $\Delta T$ is unknown beyond the fact that it should be a constant number in VSS. One way to avoid the difficulty of not knowing $\Delta T$ beforehand is to extend $r_{o}$ to a sufficiently large value that $\Delta T\ll\dot{M}\ell$, in which case the exterior VSS solution is $\Sigma\approx\Sigma_{Z}$ (as was done in e.g., Miranda et al., 2017; Muñoz et al., 2019). But rather than making the computational domain so large, we note that in the exterior VSS solution $dF_{\nu}/dr=\dot{M}d\ell/dr$ (Eq. 15 with $t_{\rm dep}=0$), which provides a condition on the gradient of $\Sigma$ at $r_{o}$—given an input value for $\dot{M}$. To apply this condition, we set in the ghost cells (at $r>r_{o}$) $F_{\nu}=F_{\nu}|_{r_{o}}+(\ell-\ell|_{r_{o}})\dot{M}$, i.e., we set $\Sigma$ in the ghost cells according to the relation $3\pi\nu\Sigma\ell=(3\pi\nu\Sigma\ell)|_{r_{o}}+(\ell-\ell|_{r_{o}})\dot{M}$. To set the value of $v_{r}$ in the ghost cells, we then use $v_{r}=-\dot{M}/(2\pi r\Sigma)$, and for $v_{\phi}$ we set it to the pressure-corrected and extrapolated Keplerian value, as for the inner disk. As was the case for the inner boundary, we enforce a wave-killing-zone near the outer boundary between $r_{o,wkz}<r<r_{o}$.

We note that Miranda et al. (2017) used, for a subset of their circumbinary disk simulations, an outer boundary condition similar in spirit to ours: they fixed $\Sigma$ at the outer boundary to the value given by Eq. (19), where the $\Delta T$ was measured from a previous iteration of the simulation.

3.2 Iterative approach to VSS

In order to reach VSS we must integrate the equations of motion for several viscous times at the outer boundary, which is prohibitively long. For $\alpha\sim 10^{-4}$ an outer viscous time is $\sim 5$ million planet orbits. To run that long on, for example, one K80 GPU takes a wallclock time of over a year, at our typical timestep of 0.003 planet orbits. Therefore, for our small $\alpha$ simulations we adopt an iterative approach. At each iterative step, we start from a profile for $\langle\Sigma(r)\rangle$. We then run a FARGO3D simulation to WSS, which is much shorter than VSS (§2.1). The result of that simulation determines $t_{\rm dep}$, which in turn determines $\langle\Sigma(r)\rangle$ from the VSS equation (Eq. 17), and hence can be used to initiate the next iterative step.

We consider a simulation to have reached VSS when the time-averaged $\dot{M}$ throughout the domain is within 10% of the forced $\dot{M}$ at the outer boundary. In our simulations with $q\gtrsim 2\times 10^{-3}$ the disk becomes eccentric (Goodchild & Ogilvie, 2006; Kley & Dirksen, 2006; Kley et al., 2008; Fung et al., 2014; Teyssandier & Ogilvie, 2017). The disk eccentricity makes defining a steady-state difficult as there is a long precession timescale in the disk and it is uncertain if the disk eccentricity should persist in steady-state. For these reasons we simply omit these simulations from our analysis, leaving a more detailed study of such disks to future work.

3.3 Resolution

In all of our simulations, the computational domain extends from $(r_{i},r_{o})=(0.3,3.68)$ with uniform spacing in azimuth and $\ln(r)$. The edges of the wave-killing zones are at $(r_{i,wkz},r_{o,wkz})=(0.46,3.0)$. The number of grid cells in each dimension is $N_{\phi}\times N_{r}=1005\times 401$, which provides near square cells (i.e. $\Delta r\sim r\Delta\phi$) and corresponds to eight cells per scale-height. To check convergence, we have also run each simulation at the cruder resolution $N_{\phi}\times N_{r}=502\times 200$. We find that the low resolution total torque agrees with the high resolution result to within $\sim 10\%$, on average, and $\sim 30\%$ in the worst case (see Table 1).

4 Numerical Results

We run a suite of FARGO3D simulations to viscous steady state. Figure 2 shows the $q$ and $\alpha$ values that we cover, with each filled circle representing a converged VSS simulation. Table 1 in Appendix D summarizes the simulation results. Note that we name simulations according to their $q$ and $\alpha$ values in an obvious notation; e.g., simulation “q1x3a3x4” has $q=10^{-3}$ and $\alpha=3\times 10^{-4}$.

For some purposes below, it will prove convenient to group simulations according to their value of $K$, as simulations with similar $K$ values have similar gap depths and one-sided torques (§2.3.2). Figure 2 shows that our chosen parameters group into clusters with nearly (though not identically) the same values of $K$.

The ultimate result from these simulations is $\Delta T$, the values of which are listed in Table 1. They are also plotted versus $K$ below (Figure 10). But we refrain from a discussion of $\Delta T$ until after we have described the simulation results in more detail.

4.1 Standard simulation overview

We focus first on a single “standard” simulation, q1x3a1x3 (i.e., $q=\alpha=10^{-3}$ implying $K=3200$). Its pileup factor is $\Delta T/(\dot{M}\ell_{p})=2.0$. In Figure 3, we show the 2D VSS surface density for this simulation with several gas streamlines overplotted. One may observe the deep gap ($<0.1\%$) surrounding the planet, with trailing spiral arms visible in the inner and outer disks.

In Figure 4 we show the principle torque balances from §2 for this simulation. All of the quantities shown have been averaged over $3,000$ orbits of the planet. Details of the torque calculations and averaging procedure are given in Appendix C. In each panel, the vertical solid lines mark the extent of our computational domain, and the dashed vertical lines mark the start of the wave-killing regions.

We now walk through each of the panels. Panel (a) shows the azimuthally averaged $\Sigma$ profile. Of particular note is the gas pileup where $\Sigma$ is roughly a factor of two larger than $\Sigma_{Z}$. Panels (b) and (c) show the torques from the WSS equation (Eq. 10). In panel (b) we show the differential torques $t_{\rm ex}$, $t_{\rm dep}$, and $dF_{\rm wave}/dr$ while in panel (c) we show the integrated torques (or fluxes), $T_{\rm ex}\equiv\int_{0}^{r}dr^{\prime}\,t_{\rm ex}(r^{\prime})$ , $T_{\rm dep}\equiv\int_{0}^{r}dr^{\prime}\,t_{\rm dep}(r^{\prime})$, and $F_{\rm wave}$. In all cases, we compute the deposited torque profile using Eq. (10) rather than Eq. (A9). In computing the torque profiles we neglect the contribution from waves with $m=1$.141414More specifically, we omit the contribution from $m=1$ to the sums defined in Appendix C (Eq. C2). We do this simply for aesthetic reasons: the $m=1$ contribution to $t_{\rm ex}$ is highly oscillatory, but hardly affects $t_{\rm dep}$ —as shown by the dashed black line in panel (b), which plots the $m=1$ contribution to $t_{\rm dep}$.

Most of the excited torque comes from near the peak of the $t_{\rm ex}$ profile. More precisely, the one-sided torques are dominated by the values at the inner and outer peaks of the torque per unit logarithmic distance, $xt_{\rm ex}$, where $x=r-1$ is the distance to the planet. For this simulation the inner and outer peaks occur at $x_{-}\approx-0.18$ and $x_{+}\approx 0.23$, shown as the dotted vertical lines, respectively. At larger distances, $t_{\rm ex}$ becomes oscillatory due to the dominance of isolated low-$m$ Lindblad resonances. However, the torque from these oscillatory regions mostly cancels, as may be seen in the plot of $T_{\rm ex}$.

Panel (c) illustrates the distinction between torque excitation and deposition (§2.1): the $T_{\rm dep}$ profile is broader than $T_{\rm ex}$, because waves transport angular momentum away from the planet.

Panels (d) and (e) show the torques from the VSS equation (Eq. 15), with $t_{\rm dep}$ repeated from earlier panels. That the three torques nearly sum to zero in panel (e) illustrates that our simulation has reached VSS. In fact, the deviation from zero in that panel is mainly due to the neglect of the $m=1$ mode. The detailed shape of these profiles near the planet play an important role in determining $\Delta T$ when the gap is very deep—a point we return to in §5.4. Panel (f) shows that $\dot{M}$ is nearly constant throughout the disk, implying that mass transport has reached steady state—in addition to angular momentum transport.

We convert the $T_{\rm dep}$ profile in panels (c) and (e) into a $\Sigma$ profile via the VSS equation (Eq. 16), with $F_{\nu}=3\pi\nu\Sigma l$, i.e., ignoring non-Keplerian contributions, and plot the result in panel (a) as a blue-dashed line. The agreement with the true $\Sigma$ profile is excellent, except for a small disagreement near the bottom of the gap where the non-Keplerian effects are evidently important.

For this simulation, the outer wave killing zone has little effect, because the waves have already damped before reaching $r_{o,wkz}$, as evidenced from the fact that both $t_{\rm dep}$ and $F_{\rm wave}$ are nearly zero by then. Conversely, the inner wave killing zone has a dramatic effect on the $\Sigma$ profile: it forces it to rise to $\Sigma_{Z}$ across an artificially short distance. But one may see that this artificiality has negligible effect on the value of $\Delta T$, or on quantities such as the depth of the gap. In a realistic disk with no wave killing zone and $r_{i}\rightarrow 0$, the waves would deposit their angular momentum at smaller $r$, resulting in a more gradual rise of $\Sigma$ inwards. But the same amount of angular momentum would still be deposited, because our artificial wave-damping prescription conserves angular momentum. In other words, the (non-wave-killing) computational domain need only capture most of the wave excitation rather than the wave deposition, in order to correctly determine the torques, and hence $\Delta T$. To illustrate this point further, the black circles in panel (a) show the surface density calculated from Eqs. (10) and (17) using the values of $F_{\rm wave}$ at the wave-killing boundaries. These agree with the true surface density profile. Nonetheless, we emphasize that our $\Sigma$ profile is incorrect at $r<r_{i,wkz}$, and the resulting error will be seen to be more dramatic in some of our other high-$K$ simulations.

4.2 Radial profiles at different $q$ and $\alpha$

In Figure 5 we show the $\Sigma/\Sigma_{Z}$, $t_{\rm ex}$, and $t_{\rm dep}$ profiles for our VSS simulations with $K\gtrsim 30$. We group simulations by their $K$, even though the value of $K$ varies slightly within each group. Each simulation is colored by its $\alpha$ value – a scheme which we adopt for the remainder of the paper. In plotting the $t_{\rm ex}$ and $t_{\rm dep}$ profiles we have removed the contribution from within the planet’s Hill sphere (which has a radius of $(q/3)^{1/3}$) for $q=2\times 10^{-3}$, and replaced the missing bit with dashed lines. We do so because the profiles have large spikes near the planet that hide the other profiles. But these tend to cancel and so are likely not of great importance. One deduces the following from these figures:

• •

$K$ is an excellent ordering parameter: simulations with similar $K$ have similar gap depths and torque density profiles. This is to be expected from the theory for moderately deep gaps (§2.3.2), but it continues to hold true for very deep gaps ($K\gtrsim 100$). Furthermore, as expected, simulations with larger $K$ tend to have both deeper gaps and larger pileups. The largest pileup we find is $\sim 10$, at $K\sim 10^{4}$.

• •

As before, vertical lines in the $\Sigma$ plots show the locations where inner and outer excited torques predominantly come from ($x_{\pm}$), defined as where $xt_{\rm ex}$ reaches its inner and outer extrema. We argued in §2.3.2 that these locations are of key importance. We see that higher $K$ systems excite their torque farther from the planet. Figure 6 shows these locations for all of our simulations. Evidently, simulations with $K\lesssim 100$ have their torque excited at the torque cutoff ($h$), and hence qualify as moderately deep gaps (§2.3.2). As $K$ increases, the excitation site is pushed further out, because the gap at $h$ becomes so deep that there is negligible wave excitation there.

• •

Simulations at a given $K$ have different $\Delta T$, as inferred from the relative heights of their pileups; i.e., their two-sided torque differs even though their one-sided torque is quite similar. This is only superficially paradoxical, because the differences in one-sided torques across different simulations become amplified in forming the two-sided torque. The sense of variation is that, at fixed $K$, simulations with lower $\alpha$ (and hence lower $q$) have larger pileups (and hence higher $\Delta T$). As we shall show below, this trend is systematic, and is not caused by the variation of $K$ within each group.

4.3 Torque excitation

In order to calculate $\Delta T$ from first principles, one may proceed in an iterative way. First, given a background surface density profile, one determines the excited torques ($T_{\pm}$). Second, one calculates where that excited torque is deposited, after it has been carried further away from the planet by waves (i.e., $t_{\rm dep}$); $t_{\rm dep}$ then determines the $\Sigma$ profile via the VSS equation (Eq. 17). Finally, one uses that new $\Sigma$ profile to calculate the new $T_{\pm}$, and iterates until convergence. In this subsection, we focus on the first step. In particular, we show that given $\Sigma$ one may predict $T_{\pm}$ quite simply—without running a full hydrodynamical simulation. Somewhat surprisingly, the case of a very deep gap is even simpler than that of a moderate gap.

Following, e.g., Goldreich & Tremaine (1980) (but see also Artymowicz, 1993; Korycansky & Pollack, 1993; Ward, 1997; Tanaka et al., 2002; Rafikov & Petrovich, 2012; Petrovich & Rafikov, 2012), we calculate $T_{\pm}$, given $\Sigma$, under the assumption that the waves are linear. We therefore linearize the equations of motion, and solve them numerically, subject to outgoing boundary conditions. See Appendix B for details. This is similar to what we have done in §2.3.2, except here we use the $\langle\Sigma\rangle$ and $\left\langle\Omega\right\rangle$ profiles from the hydrodynamical simulation in VSS as the background. The top left panel of Figure 7 compares the profile of $t_{\rm ex}/\langle\Sigma\rangle$ from the linear solution (dashed line) with that from one of our hydrodynamical simulations with $K\sim 10$. The agreement is almost perfect, because the waves launched in the simulation are indeed linear at this modest value of $K$. This demonstrates that one need not solve a full hydrodynamical simulation to obtain $T_{\pm}$ for this value of $K$—only the much simpler linear solution is needed (even though it is still numerical). The lower left panel repeats the exercise, but for a simulation with $K\sim 3,000$. Now, the linear and hydro solutions disagree close to the planet, demonstrating that the waves are very nonlinear there. But near where the torques are excited—i.e., in the vicinity of $x_{\pm}$—the linear and hydro solutions agree quite well. Hence for this simulation, too, the simple linear solution suffices to predict $T_{\pm}$, once $\left\langle\Sigma\right\rangle$ is specified.

One might wish for an even simpler— and purely analytic—prediction for $T_{\pm}$. In Appendix E, we derive a simple extension to the “standard torque formula” of Goldreich & Tremaine (1980) that accounts for the leading asymmetry between inner and outer torque at large distances from the planet, starting from the more general expression derived by Artymowicz (1993) and Ward (1997). Our result is

[TABLE]

where $C\approx 2.5$; the above expression is independent of $h$. This is shown as dotted lines in the left panels of Figure 7. As seen in the figure, this formula fails near the torque cutoff, and hence is inadequate to explain $T_{\pm}$ in the low-$K$ simulation. But it matches the high-$K$ simulation well at $x_{\pm}$. To see the behavior more clearly, in the right panels of the figure we re-plot on a log-scale, and also add our other simulations at the two $K$ values. The values of $t_{\rm ex}/\left\langle\Sigma\right\rangle$ at $x_{\pm}$ (indicated by the circles) are close to the values from Eq. (30) near $x_{\pm}$. Therefore, for high-$K$ one may predict $T_{\pm}$ by multiplying Eq. (30) at $x_{\pm}$ by the surface density at those locations ($\Sigma_{\pm}$):

[TABLE]

after dropping an order-unity coefficient. In §5.4 below, we compare this prediction for $T_{\pm}$ with the actual values in all of our $K\gtrsim 100$ simulations.

To summarize, we have shown that the linear calculation suffices to determine $T_{\pm}$ in all of our simulations, and the much simpler standard torque formula (with the added asymmetry) suffices for the high-$K$ simulations. The latter result might appear surprising in light of studies showing that the standard torque formula can be quite inaccurate—it can even give the wrong sign for $t_{\rm ex}$ at certain distances from the planet (e.g., Dong et al., 2011; Rafikov & Petrovich, 2012; Petrovich & Rafikov, 2012, see also the right panels of Figure 7 at $|x|\gg|x_{\pm}|$). Nonetheless, those inaccuracies evidently do not have a large effect on $T_{\pm}$—at least for the range of parameters spanned by our simulations.

4.4 Separating Lindblad from Co-orbital Torques

We shall separate Lindblad from co-orbital torques in the simulations, in order to show that (i) Lindblad torques are well-understood for moderate gaps, and (ii) co-orbital torques are usually sub-dominant, across all simulations. Before doing so, we describe here how we separate out the two torques.

Previous treatments have separated the torques by calling torque excited inside the horseshoe zone the co-orbital torque, and that excited outside the Lindblad torque (e.g., Paardekooper & Papaloizou, 2009). However, by examining 2D plots of $t_{\rm ex}$ (not shown), such a distinction appears ambiguous: there is no clear boundary separating one type of torque from another. Instead, we have found that the distinction becomes much clearer when examining $t_{\rm dep}$. Figure 8 shows 2D maps of $t_{\rm dep}$ for two simulations with $K\sim 1,000$. The black dashed lines show the separatrices. These mark the transition from librating to circulating streamlines in the planet’s co-rotating frame. The zoomed-in insets of Figure 8 show that the distinction between Lindblad and co-orbital torques is quite apparent: Lindblad torques show up as outwardly projecting arms, and co-orbital torques as the nearly elliptical structure near the planet (caused by the U-turn of fluid that follows nearly horseshoe orbits). One may understand why $t_{\rm dep}$ is more useful for separating out the two torques as follows: before contributing to $t_{\rm dep}$, the Lindblad torque propagates away from its point of excitation along the spiral arms, away from the co-orbital zone. Figure 8 also shows that the separatrix is only an approximate dividing line between the spiral-type and elliptical-type pattern. As such, we separate the two contributions by eye for each simulation.

4.5 Total torques

We present here the torques from our suite of VSS simulations. Figure 9 shows the measured torques as a function of $K$. The left panel shows the one-sided inner (leftwards pointing triangles) and outer Lindblad torques (rightwards pointing triangles). The right panel shows the final $\Delta T$ (solid points), which is what is relevant for the pileup and the planet’s migration; it also shows the Lindblad component of $\Delta T$ (open points). For the most part, the co-orbital torques constitute a tens of percent correction to the total $\Delta T$. However, they are significantly more important in the high $\alpha$ simulations ($\alpha\gtrsim 3\times 10^{-3}$).

In both panels, we plot as solid red lines our “moderate gap” predictions for the Lindblad torques from §2.3.2. The agreement is excellent at $K\lesssim 100$, both for the one-sided torques, and for the Lindblad component of $\Delta T$. We also show in the right panel the canonical Type I scaling (which corresponds to the no gap limit, $K\ll 100$) as the dashed red lines, taken from Tanaka et al. (2002). The upper one is for Lindblad only, and agrees with the prediction (and simulations). The lower one includes co-orbital torques, and appears to provide a better match to the low $K$ simulations151515 Our $K\ll 100$ Lindblad prediction for $\Delta T$ is the same as that of Tanaka et al. (2002) because the $C_{\pm}$ we calculated from linear theory (Eqs. 22–23) differ negligibly from theirs. That same linear calculation also produces the corotation torque, and we have verified that we get the same result as Tanaka et al. (2002) for that as well. .

Proceeding to $K\gtrsim 100$, we see that the one-sided inner torques asymptote to values of $T_{-,LR}\approx-\dot{M}\ell_{p}$, while $T_{+,LR}>\dot{M}\ell_{p}$. In VSS, we know that the inner torque cannot exceed $\dot{M}\ell_{p}$ and deviates from $\dot{M}\ell_{p}$ by the gap depth (Eq. 17 at $r=1$), which for our highest $K$ simulations is less than $0.1\%$ (Figure 5). Unlike the inner torque, the outer torque has no such restriction. Clearly, the moderate-gap $\Delta T$ prediction is no longer valid once $K\gtrsim 100$, as $\Delta T$ continues to increase with $K$ sub-linearly. By contrast, Kanagawa et al. (2018) find that $\Delta T$ is nearly constant (for fixed $h$) at those $K$ values. As we show below in §5.3, this discrepancy is due to the Kanagawa et al. (2018) results not being in VSS.

From the right panel, we see that the $\Delta T$’s from simulations with the same $\alpha$ (i.e., with the same color) trace out distinct lines, with the height of the line dropping with increasing $\alpha$. This implies that there is variation in $\Delta T$ at fixed $K$, with higher $\alpha$ simulations having lower $\Delta T$, as already suggested by the pileups seen in Figure 5. Motivated by this, we fit $\Delta T$ to a power law in $q$ and $\alpha$ above $q=10^{-4}$. The result is,

[TABLE]

where the errors are statistical. We have chosen to omit the $K\sim 10^{4}$ point from our fit as it has an unusually large $\Delta T/(\dot{M}\ell_{p})$. Given that the dependence is nearly $\propto\dot{M}\ell_{p}\alpha^{-1}$, we show in Figure 10, $\Delta T\times\alpha$ as a function of $q$ (not $K$). The lines in the figure show the best fit $\Delta T$ values for each of the $\alpha$ values.

Given $\Delta T$, we can calculate the expected pileup magnitude and $\Sigma$ profile outside of the location where $T_{\rm dep}\approx\Delta T$ using Eq. (19),

[TABLE]

In the surface density profiles of Figure 5 we show the value of the pileup at r=3.5 calculated from our $\Delta T$ scaling (Eq. 32) with horizontal arrows. These are in good agreement with the true $\Sigma$ values from the simulations.

5 Additional Features of the VSS Solutions

5.1 Planet Migration Rate and Validity of the VSS Assumption

The two-sided torque $\Delta T$ must come at the expense of the planet’s angular momentum. Hence for positive $\Delta T$, the planet migrates inwards at a rate

[TABLE]

From the results of §4, there are two different regimes for $\Delta T$: a moderate gap regime, and a deep gap regime.

For moderate gaps, $\Delta T$ is given by Eq.  (28)161616We again ignore the contribution from co-orbital torques as they are a minor correction (Figure 9)., which we rewrite as

[TABLE]

where the latter expression is specialized to $h=0.05$. The migration rate is therefore

[TABLE]

where $M_{d}=4\pi r_{p}^{2}\Sigma_{Z,p}$ is a measure of the local disk mass and $\tau_{\rm orb}=2\pi/\Omega_{p}$ is the orbital period of the planet. This is very similar to the standard Type I rate, aside from the extra gap reduction factor in the denominator. More precisely, the standard Type I rate has a coefficient of $2.3$ at $h=0.05$ if one ignores co-orbital torques; co-orbital torques change it to 1.6. (e.g., Tanaka et al., 2002; Kley & Nelson, 2012). Our new Type I migration rate includes the reduction effect of the gap; from Figure 9, it is valid at $K\lesssim 100$.

In the deep gap regime the situation is quite different, as the torque scaling switches to Eq. (32) for the $K\gtrsim 100$ simulations that we have run (up to $K\lesssim 10^{4}$). Using that fit to $\Delta T$, the migration rate is

[TABLE]

again for $h=0.05$. Remarkably, we find that in VSS planets migrate at a rate which is roughly independent of their mass and the disk’s viscosity and is only dependent on the disk-to-star mass ratio.

It is instructive to compare the VSS migration rate above with prior Type II results. These typically predict that the planet migrates at the same rate as the disk’s viscous accretion rate, although sometimes with an additional mass reduction factor when the disk is less massive than the planet (Syer & Clarke, 1995; Ward, 1997; Ivanov et al., 1999; Edgar, 2007; Armitage, 2010). The general migration rate expression (Eq. 34) may be rewritten as

[TABLE]

where the disk’s viscous time is $\tau_{\rm visc}=r_{p}^{2}/\nu_{p}$ and the bracketed factor is the dimensionless pileup factor determined from our simulations (Figure 9). Therefore the planet’s migration rate differs from the disk’s viscous accretion rate by two factors: $M_{d}/M_{p}$ and the pileup factor.

We conclude this subsection by examining the criterion for VSS to be valid. The basic assumption for VSS is that the planet migrates more slowly than the disk material (see also Kocsis et al., 2012a, b). For pileup factors ($\Delta T/(\dot{M}\ell_{p})$) that are of order a few or less, one therefore requires $|\dot{r}_{p}/r_{p}|\lesssim 1/\tau_{\rm visc}$, which implies from Eq. (39) that $M_{d}\lesssim M_{p}\left(\Delta T/\dot{M}\ell_{p}\right)^{-1}$. In other words, for order-unity pileups the VSS assumption is valid when the disk is less massive than the planet. For much larger pileups there is a more stringent constraint, because material at the peak of the pileup moves more slowly that $r_{p}/\tau_{\rm visc}$. But since the biggest pileups that we have found are $\sim 10$, we shall not consider very large pileups here.

5.2 Gap depth and width

In this subsection we compare our gap depths and widths to previously published results, and provide a new gap depth scaling that matches our VSS simulations.

The top panel of Figure 11 shows the gap depths for all of our simulations. To calculate the gap depth we first remove the circumplanetary disk region, which we crudely define as the region where both $|x|$ and $|\phi|$ are $<{\rm max}(h,(q/3)^{1/3})$, and then take the azimuthally-averaged surface density at the planet. The dotted line shows the prediction for moderately deep gaps, $\Sigma_{p}=1/(1+0.04K)$ from Eq. (26), which agrees with previous studies (e.g., Kanagawa et al., 2015b; Duffell, 2015; Kanagawa et al., 2017). Above $K\sim 300$ our gap depths are significantly deeper than Eq. (26) due to the separation of $|x_{-}|$ and $h$ (cf. Figure 6 and §5.4). A similar “two-step" effect has also been seen for gaps around low-mass ($q\lesssim 10^{-5}$), low-viscosity ($\alpha\lesssim 10^{-3}$) planets (Ginzburg & Sari, 2018). We fit our gap depths to a corrected scaling relation, $\Sigma_{p}=1/(1+0.04K+(K/K_{c})^{2})$, where $K_{c}=180$ is a fit parameter.

The bottom panel of Figure 11 shows gap widths, defined as the width of the region where $\Sigma_{p}<0.5\Sigma_{Z,p}$. Kanagawa et al. (2016) (see also Kanagawa et al., 2017) empirically determined that the gap width follows $\Delta_{g}=0.41(h^{2}K)^{1/4}$, which we show as the dotted line. For low viscosity disks, we find more radially extended gaps than predicted from that relation. From Figure 5, we see that such “extra wide” gaps are very asymmetric with respect to the planet. In fact, for many of our high-$K$ simulations we find that the inner boundary of the gap extends past our inner wave-killing boundary, and hence in reality could be much wider than found in our simulation (as explained at the end of §4.1).

An alternative empirical gap depth and width relation that separates the $q$ and $\alpha$ dependence has recently been developed by Duffell (2019). We find that his relation matches our depths in the low and intermediate $K$ regime, and, in particular, reproduces the variation seen at fixed $K$. Only at our highest $K$ and lowest $\alpha$ values ($\alpha=10^{-4}$ at $K\sim 3000$ and $\alpha=3\times 10^{-4}$ at $K\sim 10^{4}$), does his relation overpredict the depth of the gap.

5.3 Comparison with $\Delta T$’s from previous work

In §2, we argued that simulations must have correct boundary conditions in order to produce the correct pileup in surface density. But we have also argued (§4 and §5.4) that the value of $\Delta T$ should be set largely by what happens where the torque is excited ($x_{\pm}$), which occurs quite close to the planet. Therefore, prior simulations that did not adopt correct boundary conditions might produce values of $\Delta T$ that are comparable to ours.

Figure 12 compares our values of $\Delta T$, with those from simulations by Dürmann & Kley (2015) and Kanagawa et al. (2018). In contrast to our VSS boundary conditions, Kanagawa et al. (2018) set all fluid quantities equal to their initial conditions, which corresponds to the ZAM solution, at the outer boundary, and use an open boundary condition at the inner boundary. They also use wave-killing regions near both boundaries where they damp all quantities to their initial conditions. Dürmann & Kley (2015) also fix all quantities to their initial conditions at the outer boundary, but fix only $v_{r}$ and $v_{\phi}$ to their initial conditions at the inner boundary. For their wave-killing prescription, they damp only $v_{r}$ and $v_{\phi}$ to their initial conditions and not $\Sigma$. An additional difference is that Kanagawa et al. (2018) allow their planets to migrate, as do Dürmann & Kley (2015) for a subset of their simulations.

Focusing on $K$ values above $100$, we see that the non-migrating Dürmann & Kley (2015) torques are quite close to our VSS torques with the exception of the $K\sim 1,000$ planets, suggesting that the boundary conditions are not essential to achieving the correct $\Delta T$, for these $K$ values. By contrast, the Kanagawa et al. (2018) torques tend to be quite different than our VSS torques, as do the migrating planets of Dürmann & Kley (2015), suggesting that those migrating planets are not in VSS—likely due to the disk-to-planet mass ratio being too large (§5.1).

5.4 Towards a Theory of Very Deep Gaps

The theory for moderately deep gaps (§2.3.2) assumes that torque excitation happens at the torque cutoff, i.e., that $|x_{\pm}|\sim h$. But once the gap becomes sufficiently deep, torque is mostly excited further away from the planet, i.e., in the gap wall, where the surface density is higher. As mentioned previously, to predict where the gap wall occurs requires knowing where torque is deposited, because that determines $\Sigma$ via the VSS equation (Eq. 17). But torque deposition is difficult to calculate from first principles. Therefore we cannot yet present a complete theory for very deep gaps. Nonetheless, we are able to take a few steps towards such a theory.

First, as suggested in §4.3, one may obtain the excited torques quite simply given the $\Sigma$ profile. In fact, one only needs four key numbers characterizing the $\Sigma$ profile: the aforementioned $x_{-}$ and $x_{+}$, as well as the $\Sigma$ values at $x_{\pm}$, i.e., $\Sigma_{-}$ and $\Sigma_{+}$. We illustrate these important quantities in Figure 13. Assuming that $t_{\rm ex}/\Sigma$ follows Eq. (30) , $x_{\pm}$ correspond to the locations where $d\ln\Sigma/d\ln|x|\approx 3$, and furthermore, the excited torques, $T_{\pm}$, follow from Eq. (31) once the four numbers are known. Figure 14 compares this prediction for $T_{\pm}$ with what is found in the simulations, where the “prediction” makes use of the values of $x_{\pm}$ and $\Sigma_{\pm}$ extracted from the simulations. Above $K\sim 100$, the agreement is quite good, particularly at $K\gtrsim 3,000$, thus confirming our claim that given $\Sigma$, one can calculate the one-sided torques.

Second, as we now show, $\Sigma_{-}$ may be determined by integrating the VSS equation from $r=0$ to $r=r_{p}+x_{-}$. The calculation is nearly the same as for the moderate gap case (cf. the discussion surrounding Eq. 26). Starting from Eq. (24), but taking the upper limit of the integral to be $r_{p}+x_{-}$ yields

[TABLE]

where $\Sigma_{Z,-}$ is $\Sigma_{Z}$ at $x_{-}$ and where171717Note that the asymmetric coefficient here is $1.26$ as opposed to $2.26$ because there is a factor of $r$ when converting between $F_{\nu}$ and $\Sigma$ for Keplerian disks.

[TABLE]

Equation (40) is the extension of Eq. (26) to very deep gaps, but now it only provides a consistency relation between $\Sigma_{-}$ and $x_{-}$. In Figure 15, we show that the measured values of $\Sigma_{-}$ agree well with the values provided by Eq. (40) given $x_{-}$, for simulations with $K>100$.

We have now reduced the problem to three unknowns, $x_{-}$, $x_{+}$, and $\Sigma_{+}$. One may additionally integrate the VSS equation between $x_{-}$ and $x_{+}$ to relate the jump in surface density ($\Sigma_{+}-\Sigma_{-}$) to the torque deposited between $x_{-}$ and $x_{+}$. In other words, if one can determine the torque deposited between $x_{-}$ and $x_{+}$ (which is likely highly non-local and non-linear), as well as the values of $x_{-}$ and $x_{+}$, one will have a complete theory.

5.5 Models based on local deposition are inadequate, particularly at high $K$

A number of previous papers (Liu & Shapiro, 2010; Kocsis et al., 2012a, b) have constructed 1D models for what we call VSS. However, these are based on the assumption that $t_{\rm ex}=t_{\rm dep}$, i.e., they ignore the fact that waves transport angular momentum from where they are excited (at $\sim x_{\pm}$) to where they are deposited. While that might seem a minor point, it leads to extremely erroneous results, as we demonstrate briefly here.

To show this, we set $t_{\rm dep}=t_{\rm ex}$ in the VSS equation (Eq. 15), and use Eq. (30) for $t_{\rm ex}$. For simplicity, we also set $\dot{M}=0$ and $\nu\ell={\rm const}$ such that $F_{\nu}=3\pi\alpha h^{2}\Sigma$. These approximations are purely for demonstration purposes since the full Eq. (15) yields very similar results. It is straightforward to show that the solution of the VSS equation in the inner disk yields $\Sigma=e^{-f_{-}K(h/|x|)^{3}}$, where we have set $\Sigma$ at the inner boundary to one and $f_{-}\approx 0.08$ from Eq. (30) evaluated at $x=-h$. Similarly, in the outer disk $\Sigma=\Sigma_{p}e^{f_{+}K(h/x)^{3}}$, where $f_{+}\approx 0.1$. To connect the inner and outer disk we assume that $F_{\nu}={\rm const}$ between $x=-h$ and $x=+h$ (as was done in e.g., Kocsis et al., 2012a, b). The total torque for the local model is then $\Delta T\approx e^{(f_{+}-f_{-})K}-1$. We see that the gap becomes exponentially deep and the torque exponentially large for $K\gtrsim 50$. We find this same divergence of $\Delta T$ in the full Eq. (15) with $\dot{M}\neq 0$ and retaining all $r$ dependencies (a similar result is found in Liu & Shapiro, 2010). Such a divergence is incorrect. For example, at $K\sim 10^{4}$, the local model would predict $\Delta T/\dot{M}\ell_{p}\sim e^{200}$ for the pileup factor, whereas we find a value $\Delta T/\dot{M}\ell_{p}\sim 10$, an enormous discrepancy. We may conclude that local deposition is grossly inadequate, particularly at large $K$181818We note that Kocsis et al. (2012a, b) do not find such large $\Delta T$ values because they both reduce the coefficient of $t_{\rm ex}$, as well as include radially dependent $\alpha$ and $h$ profiles. .

6 Summary

We examined the planet-disk interaction problem in disks of low enough mass that the planet’s migration time is slower than the disk’s viscous accretion time. Our main results are as follows:

• •

One may study such disks by treating the planet’s orbit as fixed, and examining the disk’s properties in viscous steady state (VSS). This is a particularly clean setup to study the planet-disk interaction problem. The key question becomes what is the total torque injected by the planet ($\Delta T$) in VSS, for a particular set of problem parameters (principally, $q,\alpha,h$)? The value of $\Delta T$ determines both the pileup of disk material exterior to the planet’s orbit, and the migration rate of the planet.

• •

We predicted $\Delta T$ for moderately deep gaps (§2.3.2). We then ran a series of hydrodynamical simulations that reached VSS for a variety of parameters. The results of the simulations agreed with the theory for moderately deep gaps. But for very deep gaps, the theory is inadequate. Empirically, for very deep gaps our simulations yielded the approximate relation $\Delta T/\dot{M}\ell_{p}\approx 4(q/\alpha)$ when $q\gtrsim 10^{-4}$.

• •

We calculated the resulting planet migration rate, showing how the well-understood Type I rate smoothly transitions into a new Type II rate as the gap formed by the planet becomes increasingly deep.

7 Open Questions

We have left a considerable number of open questions to future investigations. Some of these are as follows.

• •

What is the VSS result for parameter values not examined in this paper, and is it possible to achieve a pileup factor larger than the largest we found in our simulations (i.e., $\Delta T/\dot{M}\ell_{p}\sim 10$ for $K\sim 10^{4}$)? Our simulations only explored a limited range of parameters: we set $h=0.05$ and $q$ and $\alpha$ along the grid of filled circles in Figure 2. We expect that both higher $q$ and lower $\alpha$ might lead to a higher pileup factor. At higher $q$, we found that the disk transitioned to an eccentric state (see also Goodchild & Ogilvie 2006; Kley & Dirksen 2006; Kley et al. 2008; Fung et al. 2014; Teyssandier & Ogilvie 2017). How realistic is that result, and if it is realistic, what is the VSS for an eccentric disk? A potential difficulty is that eccentric disks behave quite differently in 2D and 3D (e.g., Ogilvie, 2008; Lee et al., 2019). Regarding lower $\alpha$, we have not been able to reach VSS for $\alpha<10^{-4}$ with our simulations because of their computational cost. Of course, if $\alpha$ is too small, the time to reach VSS might be longer than the age of the disk.

• •

Do 3D effects significantly affect the pileup? How important is accretion onto the planet? What is the effect of using a more realistic equation of state? Miranda & Rafikov (2019) show that an adiabatic (rather than locally isothermal) equation of state leads to a different $t_{\rm dep}$ profile. We suspect the change will be minor because $\Delta T$ is most sensitive to what happens very close to the planet, where the effect of equation of state is likely minor.

• •

Are the surface density profiles for disks in VSS consistent with those inferred from observations of protoplanetary disks? For example, could inferred inner holes be the result of the deficit of material within a planet’s orbit relative to a pileup outside of it? And could some of the gaps and rings imaged at large radii—that are often attributed to planets (e.g., Zhang et al., 2018)—be the result of a pileup outside of the planet? In this paper, we have only addressed gas dynamics, and to make detailed comparison with observations one must also understand how the dust behaves. Hence we leave comparison with observations to future work.

• •

How wide are the inner gaps? We find that for many of our highest $K$ simulations, the gap extends into our inner wave-killing zone. Future work should extend the inner boundary to smaller radii to determine more realistic wave deposition profile. This may prove useful for diagnosing whether an observed gap is due to a planet or by some other process.

We thank the referee, Roman Rafikov, for a thorough reading of the manuscript and many helpful comments, and Diego Muñoz for many insightful discussions. This research was supported in part through the computational resources and staff contributions provided for the Quest high performance computing facility at Northwestern University which is jointly supported by the Office of the Provost, the Office for Research, and Northwestern University Information Technology. Y.L. acknowledges NASA grant NNX14AD21G and NSF grant AST-1352369.

Appendix A Steady-state derivation

We derive the equations for the three angular momentum densities (total, wave, and mean flow) that are needed for §2.1. The 2D equations of motion for a fluid with surface density $\Sigma$, velocity $\mathbf{v}$, and pressure $P$ are,

[TABLE]

where $\Phi$ is the external gravitational field and $\mathbf{S}=\mathbf{\nabla}\mathbf{v}+\mathbf{\nabla}\mathbf{v}^{T}-2/3(\mathbf{\nabla}\cdot\mathbf{v})$ is the stress tensor. Specializing to cylindrical coordinates, $(r,\phi)$, $\Sigma$ and the specific angular momentum, $\ell=rv_{\phi}$ evolve according to,

[TABLE]

where $S_{r\phi}=r\partial_{r}\Omega+r^{-1}\partial_{\phi}v_{r}$ and, for convenience, we have combined the pressure and viscous stress into $\mathnormal{f}=-P+\nu\Sigma S_{\phi\phi}$. Together, Eqs. (A3) and (A4) describe the evolution of the total angular momentum density,

[TABLE]

and have the azimuthal averages,

[TABLE]

where $\dot{M}$ is defined in Eq. (12); $F_{\nu}\equiv-2\pi r^{2}\left\langle\nu\Sigma S_{r\phi}\right\rangle$, as displayed in Eq. (3); and the excited torque density $t_{\rm ex}$ is defined in Eq. (2). To obtain the evolution equations for the wave angular momentum $\left\langle\Sigma^{\prime}\ell^{\prime}\right\rangle$, we add $\ell^{\prime}\partial_{t}\Sigma$ to $\Sigma^{\prime}\partial_{t}\ell$ and take the azimuthal average,

[TABLE]

The term inside the radial derivative is the wave flux of angular momentum, defined in Eq. (8) and $t_{\rm dep}$ is the deposition rate of angular momentum by the waves,

[TABLE]

Note that since $t_{\rm dep}$ depends on wave quantities, it should reach a steady-state value on the timescale for the waves to reach steady-state. Finally, to obtain the evolution of the axisymmetric angular momentum we subtract Eq. (A8) from Eq. (A7), which results in Eq. (11).

A.1 Approximations

The wave flux in Eq. (8) is made of three terms defined in Eq. (9). In the left panel of Figure 16, we show the total $F_{\rm wave}$ from an example simulation (described in §4; solid line) and each of its terms. The $\propto\left\langle v_{r}^{\prime}v_{\phi}^{\prime}\right\rangle$ term (dashed line) is dominant nearly everywhere in the disk except for the region closest to the planet where the triple correlation term (dotted line) becomes dominant. As expected, the $\propto\left\langle v_{r}\right\rangle$ term (dashed-dotted line) is nearly zero throughout the disk.

The viscous flux in Eq. (3) is made of two terms,

[TABLE]

where the second term depends only on wave quantities. We show in the right panel of Figure 16 that this wave term is negligible everywhere in the disk except for very close to the planet.

Appendix B Linear solution

We describe how we solve the linear response of a disk to a planet, which is needed in §2.3.2 and §4.3. We decompose variables as follows:

[TABLE]

where we have neglected the azimuthal average of the radial velocity (see Figure 16). The gravitational potential is similarly transformed to $\phi_{m}$. To obtain the linear equations of motion we expand Eqs. (A1) and (A2) with $P=c_{s}^{2}(r)\Sigma$ to first order in $(u_{m},v_{m},\sigma_{m})$. The result is (Goldreich & Tremaine, 1979; Korycansky & Pollack, 1993; Tanaka et al., 2002)

[TABLE]

where, $\kappa^{2}=4\left\langle\Omega\right\rangle^{2}+rd\left\langle\Omega\right\rangle^{2}/dr$. The viscous accelerations, $f_{r,\phi}$, are determined by linearly expanding $(\mathbf{\nabla}\cdot(\nu\Sigma\mathbf{S}))/\Sigma$ in Eq. (A2). Once $\left\langle\Sigma\right\rangle$ and $\left\langle\Omega\right\rangle$ are specified, we solve Eqs. (B4)-(B6) as a boundary value problem for each $m$. Since the equations are linear, we use a simple matrix method where we discretize the equations onto a radial grid, specify outgoing wave boundary conditions at both boundaries, and solve the resulting tri-diagonal system of equations. For the outgoing wave boundary conditions, we assume that far from the planet the waves are in the WKB limit and follow,

[TABLE]

where $k_{m}$ is the positive root of the WKB dispersion relation, $c_{s}^{2}k_{m}^{2}=\kappa^{2}-m^{2}(\left\langle\Omega\right\rangle-\Omega_{p})^{2}$. This matrix inversion method is different than the shooting method typically used to solve the linear planet-disk equations of motion (e.g., Korycansky & Pollack, 1993; Tanaka et al., 2002; Rafikov & Petrovich, 2012; Petrovich & Rafikov, 2012).

Appendix C Numerical Appendix

Here we describe the numerical setup of our hydrodynamical simulations. FARGO3D solves Eqs. (A1) & (A2) on a staggered mesh where the density lies at the center of the cell and the velocities lie at the edges of the cell in their respective directions. For simulating accretion disks, FARGO3D uses the fast advection algorithm of its predecessor to significantly increase the CFL limited timestep by removing the dominant Keplerian azimuthal velocities (Masset, 2000). Typically, the timestep is constrained by the radial sound crossing time in the inner cells.

Our boundary conditions described in §3.1 apply only for the azimuthally averaged density and velocities. To ensure that there are no waves at the boundaries we adopt wave-killing zones (de Val-Borro et al., 2006). In these regions, we artificially enhance wave damping, such that the waves vanish at the computational boundaries. In particular, at the end of each timestep we additionally evolve the radial velocity according to

[TABLE]

where the local damping timescale $\tau=1/(30\Omega_{K}(r))$ and $R(r)$ is a quadratic function which is zero in the bulk of the domain, and rises to unity near the boundaries (de Val-Borro et al., 2006). Our choice to damp only the radial velocity ensures that we conserve both mass and angular momentum in the wave-killing zones. This is in contrast to most of the other gap-opening studies which utilize wave-killing zones that additionally damp $\ell$ and $\Sigma$ (e.g., Duffell & MacFadyen, 2013; Dürmann & Kley, 2015; Kanagawa et al., 2017).

C.1 Flux and torque calculation

Here we outline our numerical calculation of $\dot{M}$, $F_{\nu}$, $F_{\rm wave}$, $t_{\rm ex}$, presented in the main text. In short, these quantities are taken from their respective steps in the FARGO3D algorithm and a running time average is computed at each timestep as FARGO3D evolves the equations of motion over the averaging periods given in Table 1. As an example of this process, we focus here on calculating $F_{\rm wave}$ (Eq. 8) as it is the most involved. At each timestep, FARGO3D updates the angular momentum of a cell from the angular momentum fluxes in the $r$ direction. These fluxes are computed by reconstructing the cell-centered angular momenta to the radial faces of each cell (for details of this process see Benítez-Llambay & Masset, 2016). During these updates we store the reconstructed values of $\Sigma$ and $\ell$ on the cell faces, as well as the $v_{r}$ values. Using these we compute,

[TABLE]

where $(v_{r}\Sigma^{*})_{m}$ and $\ell^{*}$ are the $m$-th components of the Fourier transforms of $v_{r}\Sigma^{*}$ and $\ell^{*}$, and where the stars indicate that these values are reconstructed. Note that in the sum we only retain the $m>0$ terms. This procedure is done every timestep during the averaging period, with each timestep contributing a new value to the running average for each $m$ contribution to $F_{\rm wave}$. An analogous procedure is done for $\dot{M}$ and $F_{\nu}$ in the update functions for $\Sigma$ and $v_{\phi}$, respectively. We choose to calculate $F_{\rm wave}$, $F_{\nu}$, and $\dot{M}$ in this way so that (i) the time-averages exactly correspond to the changes in total angular momentum and mass in a given cell over the averaging period, and (ii) so that we may separate the contributions from different $m$ values.

Appendix D Simulation Table

Appendix E Ward Torque

We derive Eq. (30), which is the leading order asymmetric correction to the standard torque formula (Goldreich & Tremaine, 1980), starting from the more general torque profiles of Artymowicz (1993) and Ward (1997). Ward (1997) gives the excited torque density as (his Eq. 14, but see also Eq. 51 and 54 of Artymowicz (1993)),

[TABLE]

where $\xi=mc_{s}/(r\kappa)$, and the upper (lower) sign corresponds to the outer (inner) disk. The potential $\psi$ is,

[TABLE]

Here, $K_{0}$ and $K_{1}$ are modified Bessel functions of the second kind, $\Lambda=m|x/r_{p}|/\sqrt{(r/r_{p})}$ and $f=|\Omega-\Omega_{p}|/\Omega_{K}$. This form of $\psi$ was derived in Ward (1997) (as opposed to Artymowicz (1993)), and so from now on we shall associate this particular form of the torque density to Ward. Specializing to a sound speed profile of $c_{s}=hr\Omega_{K}$ we can rewrite $\xi$ as $\xi=mh(\Omega_{K}/\kappa)\approx mh$ for a nearly Keplerian disk. Equations (E1) and (E2) are evaluated at effective Lindblad resonances defined by $D_{\star}=\kappa^{2}-m^{2}(\Omega-\Omega_{p})^{2}+(mc_{s}/r)^{2}=0$. In terms of $m$ and $f$ this resonance condition is ,

[TABLE]

i.e. for a given distance to the planet there is a corresponding value of $m$. Note that $m$ diverges to infinity as $f\rightarrow h$ as $|x|\rightarrow 2/3h$, but the torque does not diverge due to the exponential decay of the Bessel functions with $\Lambda\propto m\rightarrow\infty$.

The Goldreich & Tremaine (1980) approximation to the excited torque follows from setting $r=r_{p}$ in Eqs. (E1) and (E2) unless it appears as $x=r-r_{p}$ in which case $|x|\approx 2r_{p}/(3m)$. With these approximations $\psi=2K_{0}(2/3)+K_{1}(2/3)$ and the torque density becomes,

[TABLE]

where the numerical constant $\mathcal{C}=(32/81)(2K_{0}(2/3)+K_{1}(2/3))^{2}\approx 2.5$. This is the $t_{\rm ex}$ profile of Goldreich & Tremaine (1980). Note that this is symmetric with respect to the sign of $x$. Ward (1997) showed that the leading order correction to the Goldreich & Tremaine (1980) torque follows from the $m^{4}\psi^{2}$ term in Eq. (E1),

[TABLE]

where $\psi_{0}=2K_{0}(2/3)+K_{1}(2/3)$. To convert $m$ to $x$ we expand the resonance condition (Eq. E3) to $(r_{p}/|x|)^{3}$ order,

[TABLE]

The final torque density with the leading asymmetry is then,

[TABLE]

In Figure 17, we compare the full Ward (1997) torque given by Eq. (E1) as a function of $x$ for a constant surface density disk against the Goldreich & Tremaine (1980) approximation given by Eq. (E4). There are a few important points to highlight here. First, there is no torque asymmetry due to $t_{\rm GT}$. Second, far from the torque cutoff, Eq. (E8) is a good approximation to $t_{\rm WW}$.