A SIMPLE ANALYTICAL MODEL FOR GAPS IN PROTOPLANETARY DISKS

Disk–satellite interactions have been studied for decades. The general problem of an orbiting point mass interacting gravitationally with a surrounding disk of either gas or solids is a fundamental physics problem with a multitude of applications in astrophysics. Saturn's rings, for example, are an exquisite testbed for disk–satellite interactions in the case where the disk is composed entirely of solids (Goldreich & Tremaine 1978). For the hydrodynamic case, how a planet interacts with the gaseous disk which spawned it is of great importance to the question of how the planets formed. It also provides a vital keystone in our interpretation of observational data, both of extrasolar planetary systems and of debris disks around newly formed stars. Additionally, disk–satellite interactions are of importance to the orbital evolution of a black hole binary embedded in a circumbinary disk.

The linear theory of hydrodynamical disk–planet interactions has been very successful so far (Goldreich & Tremaine 1978, 1980; Ward 1986). Predictions of the disk's linear response to the planet and subsequent predictions of the disk perturbation's gravitational influence on the planet have been accurately reproduced in numerical studies (e.g., Tanaka et al. 2002; Cresswell & Nelson 2006; Bitsch & Kley 2010; D'Angelo & Lubow 2010).

The details of nonlinear theory, however, have not all been reproduced in numerical studies. In linear theory, a spiral wave is excited by the planet, and passes politely through the disk. After each wave front passes, it leaves no trace of its presence. A nonlinear wave, however, can shock. When it does so, it imparts its angular momentum to the disk, causing fluid elements to move away from the planet's position. This process carves out a low-density annulus in the vicinity of the planet's orbit. This low-density region is known as a "gap," and its detailed properties have eluded a robust and numerically reproducible semi-analytical description for some time.

Many attempts have been made to describe this three-dimensional system with a one-dimensional model to predict the gap structure, i.e., the surface density Σ as a function of radius r. Varnière et al. (2004) took the disk evolution equations of Lynden-Bell & Pringle (1974) and added a planetary torque term, to predict the width and shape of gaps. The calculated gaps turned out to be much wider than those seen in numerical studies. Crida et al. (2006) attempted to fix this by adding a component due to "pressure torque" in the disk. However, both of these studies predicted gaps which were much too deep. Additionally, these studies considered the role of planetary torque as a barrier holding back the viscous accretion of gas. Such considerations predict strict criteria for gap opening (Lin & Papaloizou 1979, 1993; Bryden et al. 1999; Crida et al. 2006), but semi-analytical and numerical studies have demonstrated that planets can violate these criteria and still open gaps on long timescales (Goodman & Rafikov 2001; Dong et al. 2011a, 2011b; Duffell & MacFadyen 2012).

Similarly problematic assumptions have been made for 1D models of circumbinary disks. Disk models of Syer & Clarke (1995) and later Kocsis et al. (2012) implicitly assume that if the planet does not migrate with the gas, there is a "pile-up" of gas at the outer edge of the gap. In contrast, the disk model of Lubow & D'Angelo (2006) correctly accounted for gas flow across the gap.

Not only are many of these assumptions incorrect (see, for example, Duffell et al. 2014), but they also result in gap profiles that are inconsistent with numerical experiments. The gap profiles in these studies generally predict exponentially deep gaps; that is, the surface density near the planet is depleted by a factor $\sim \mathrm{exp}(-{(q/{q}_{\mathrm{gap}})}^{2})$ relative to the unperturbed state, where q is the planet-to-star mass ratio, and ${q}_{\mathrm{gap}}$ is some critical threshold for gap opening. On the other hand, several recent numerical parameter surveys have shown that gap depth scales as a power-law with q, the disk aspect ratio $h/r$ and the dimensionless viscosity α (Duffell & MacFadyen 2013; Fung et al. 2014), so that gaps are not as depleted as in the 1D models.

Fung et al. (2014) provided an argument for the power-law scalings found in Duffell & MacFadyen (2013), pointing out that the torque produced by the planet would be proportional to the surface density in the gap near the planet's position ${{\rm{\Sigma }}}_{p}$, assuming that the torque was dominated by resonances excited inside the gap. Very recently, Kanagawa et al. (2015) produced a 1D model for gaps which repeated the arguments of Fung et al. (2014), recovering the power-law scaling of Duffell & MacFadyen (2013).

However, Kanagawa et al. (2015) used an incomplete description of the planetary torque. For most of their study, the torque was specified as the excitation torque in the disk, as calculated by linear theory, whereas the torque that controls gap opening is the deposition torque, i.e., the rate of transfer of angular momentum from the planetary wake to the disk by shocks. They attempted to fold this into their study by including a parameter x_d that measured the distance the wave travels from the planet before shocking. The resultant gap profile was then some function of the quantity x_d, which was treated as a free parameter.

Fortunately, the process of excitation and deposition of torque in the disk has been studied by Goodman & Rafikov (2001, hereafter GR01). Therefore, there is no need to speculate about this process, at least in the weakly nonlinear regime. In the present study, the planetary angular momentum flux derived by GR01 is applied to a disk model similar to Kanagawa et al. (2015).

The model is described in Section 2. In Section 3, the resultant gap profile ${\rm{\Sigma }}(r)$ is compared with gap profiles from a recent parameter survey (Duffell 2015) and agreement is found, particularly in the "partial gap" regime. This provides confidence that a complete 1D disk model may be possible after taking into account higher-order effects. A summary and discussion are given in Section 4.

These results may be testable observationally in the not-too-distant future, as the ALMA project may be capable of measuring surface density profiles in transition disks.

The initial approach here is essentially equivalent to that of Varnière et al. (2004), though the underlying equations will be derived in a slightly different manner. Rather than using the equation of mass flux to derive a differential equation for ${\rm{\Sigma }}(r)$, a planetary torque term is inserted into the equation of angular momentum flux, and the fluxes at infinity are assumed to be unaffected by the presence of the planet. The system is then solved algebraically.

The equations of mass and angular momentum conservation in a steady-state viscous disk are as follows (ignoring the planet, at first):

$$\begin{eqnarray}&&\dot{M}=-2\pi r{\rm{\Sigma }}{v}_{r}=\mathrm{const}\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&\dot{J}=-\dot{M}{r}^{2}{\rm{\Omega }}-2\pi {r}^{3}{\rm{\Sigma }}\nu \displaystyle \frac{d{\rm{\Omega }}}{{dr}}=\mathrm{const},\end{eqnarray} \tag{ 2 }$$

where $\dot{M}$ is the mass accretion rate, $\dot{J}$ is the angular momentum flux, v_r is the radial drift velocity, Ω is the orbital angular velocity of the gas, and ν is the kinematic viscosity. All quantities are vertically and azimuthally averaged. The accretion rate is defined so that positive $\dot{M}$ implies inward accretion. The first term in (2) is the flux due to advection, and the second term is due to viscous torques.

Assuming Keplerian Ω (as will generally be done for simplicity in this study), Equation (2) becomes

$$\begin{eqnarray}&&\dot{J}=-\dot{M}{r}^{2}{\rm{\Omega }}+3\pi {\rm{\Sigma }}\nu {r}^{2}{\rm{\Omega }}.\end{eqnarray} \tag{ 3 }$$

There are two possible solutions such that $\dot{M}$ and $\dot{J}$ are both spatially uniform. The first is $\dot{M}=0$, $\dot{J}=3\pi {\rm{\Sigma }}\nu {r}^{2}{\rm{\Omega }}\;=$ constant, and the second is $\dot{M}=3\pi {\rm{\Sigma }}\nu \;=$ constant, $\dot{J}=0$, so that viscous angular momentum transport outward exactly cancels advective transport inward. The latter solution will be concentrated on in this study, as there is zero accretion in the former solution. However, the $\dot{M}=0$ solution does not pose any additional challenges, and what follows would be nearly identical (a linear combination of the two is also possible, but again this is ignored in the current study for simplicity).

Thus, for this (unperturbed) solution, the surface density Σ is inversely proportional to the viscosity, ν. If ν is spatially uniform, then Σ will also be, in steady state.

Once the constants $\dot{M}$ and $\dot{J}$ are known, one can add a planetary torque as an additional source of angular momentum flux:

$$\begin{eqnarray}&&\dot{J}=-\dot{M}{r}^{2}{\rm{\Omega }}+3\pi {\rm{\Sigma }}\nu {r}^{2}{\rm{\Omega }}+{{\rm{\Phi }}}_{p}(r)=0.\end{eqnarray} \tag{ 4 }$$

${{\rm{\Phi }}}_{p}(r)$ is the angular momentum flux generated by the planet; the torque density deposited by shocks is given by the negative derivative of this function. Mass accretion onto the planet is neglected.

Note that it is assumed that the planet does not affect the values of $\dot{M}$ and $\dot{J}$ at infinity, so that in particular the gap does not act as a "barrier" preventing mass flux through the system. The surface density ${\rm{\Sigma }}(r)$ adjusts itself to maintain these values of $\dot{M}$ and $\dot{J}$ in the presence of the planet. The value of $\dot{M}$ is the unperturbed value, the same as it would be in absence of the planet:

$$\begin{eqnarray}&&\dot{M}=3\pi \nu {{\rm{\Sigma }}}_{0}(r)\end{eqnarray} \tag{ 5 }$$

(the viscosity ν can also depend on radius, but this notation was omitted for simplicity). The planet's influence is entirely housed in the planetary angular momentum flux, ${{\rm{\Phi }}}_{p}(r)$. ${{\rm{\Phi }}}_{p}(r)$ is determined by how far the wave propagates before shocking, and how quickly its angular momentum is transferred to the disk after shocking. This process has been studied by Goodman & Rafikov (2001), and their calculations provide a clear prescription for ${{\rm{\Phi }}}_{p}(r)$, which applies in the weakly nonlinear regime. The specific details of this function will be discussed further in Section 2.1, but the overall scaling is given by the one-sided torque emanating from the planet:

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{p}(r)={{\rm{\Sigma }}}_{p}{a}^{4}{{\rm{\Omega }}}_{p}^{2}{q}^{2}{{\mathcal{M}}}^{3}f(r),\end{eqnarray} \tag{ 6 }$$

where a is the radial position of the planet, ${{\rm{\Sigma }}}_{p}$ is the surface density at the planet's position ${\rm{\Sigma }}(a)$, q is the planet-to-star mass ratio, and ${\mathcal{M}}$ is the Mach number (the inverse of the disk aspect ratio $h/r$). f(r) is a dimensionless function of r, discussed later. Near the planet, f evaluates to $f(a)={f}_{0}$, where f₀ is some dimensionless order-unity constant. Note that there is an implicit assumption that the waves are excited within the gap so that the angular momentum flux is proportional to the density in the gap, ${{\rm{\Sigma }}}_{p}$.

Evaluating (4) at the planetary position (also incorporating (5) for the mass flux) gives

$$\begin{eqnarray}&&-3\pi \nu {{\rm{\Sigma }}}_{0}{a}^{2}{{\rm{\Omega }}}_{p}+3\pi {{\rm{\Sigma }}}_{p}\nu {a}^{2}{{\rm{\Omega }}}_{p}+{{\rm{\Sigma }}}_{p}{a}^{4}{{\rm{\Omega }}}^{2}{q}^{2}{{\mathcal{M}}}^{3}{f}_{0}=0.\end{eqnarray} \tag{ 7 }$$

This can be inverted to find a formula for gap depth:

$$\begin{eqnarray}&&\displaystyle \frac{{{\rm{\Sigma }}}_{p}}{{{\rm{\Sigma }}}_{0}}={(1+{f}_{0}\displaystyle \frac{{q}^{2}{{\mathcal{M}}}^{3}{a}^{2}{{\rm{\Omega }}}_{p}}{3\pi \nu })}^{-1}.\end{eqnarray} \tag{ 8 }$$

This recovers the gap depth scaling reported in Duffell & MacFadyen (2013). This exercise was also effectively carried out by Fung et al. (2014) and by Kanagawa et al. (2015). More specifically,

$$\begin{eqnarray}&&\displaystyle \frac{{{\rm{\Sigma }}}_{p}}{{{\rm{\Sigma }}}_{0}}=\displaystyle \frac{1}{1+{f}_{0}K(q)/(3\pi )},\end{eqnarray} \tag{ 9 }$$

where

$$\begin{eqnarray}&&K(q)\equiv {q}^{2}{{\mathcal{M}}}^{5}/\alpha \end{eqnarray} \tag{ 10 }$$

and α is the dimensionless Shakura–Sunyaev viscosity parameter, $\alpha =\nu {{\mathcal{M}}}^{2}/({a}^{2}{\rm{\Omega }})$. Figure 1 compares gap depths from Equation (9) with steady-state 2D solutions from a numerical parameter survey over q, α, and ${\mathcal{M}}$ (Duffell 2015), demonstrating the accuracy of this formula over a wide range of parameter space. For these 84 disk models, the only significant deviation from Equation (9) is in the $K\sim 10-1000$ range, where higher-order effects become important. For larger values of K, Fung et al. (2014) found that this scaling eventually changes, suggesting that higher-order effects become important for very large $K\gt {10}^{4}$. This figure was used to constrain the one free parameter ${f}_{0}=0.45$. This value of f₀ is also consistent with the one-sided torque calculated numerically in these disk models (${f}_{0}\approx 0.451\pm 0.534/{\mathcal{M}}$, where the ± sign indicates the torque excited on either side of the planet, with the greater torque excited in the outer disk). After choosing ${f}_{0}=0.45$, the model is completely specified.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

2.1. Density as a Function of Radius

Now with this value of ${{\rm{\Sigma }}}_{p}$, one can return to the original $\dot{J}$ equation:

$$\begin{eqnarray}&&-3\pi \nu {{\rm{\Sigma }}}_{0}(r){r}^{2}{\rm{\Omega }}+3\pi {\rm{\Sigma }}(r)\nu {r}^{2}{\rm{\Omega }}+{{\rm{\Sigma }}}_{p}{a}^{4}{{\rm{\Omega }}}_{p}^{2}{q}^{2}{{\mathcal{M}}}^{3}f(r)=0.\end{eqnarray} \tag{ 11 }$$

Given ${{\rm{\Sigma }}}_{p}/{{\rm{\Sigma }}}_{0}$, which was calculated above (9), one can solve (11) for ${\rm{\Sigma }}(r)$:

$$\begin{eqnarray}&&{\rm{\Sigma }}(r)={{\rm{\Sigma }}}_{0}(r)\left(1-\displaystyle \frac{f(r)K(q)/(3\pi )}{1+{f}_{0}K(q)/(3\pi )}\sqrt{a/r}\right).\end{eqnarray} \tag{ 12 }$$

The function f(r) is a scaled-out version of the angular momentum flux due to the shocking of the planetary wake, which has been calculated semi-analytically by Goodman & Rafikov (2001). The wave conserves its angular momentum until it shocks and deposits its angular momentum into the disk. Therefore, the function f(r) is simply a constant, until the distance from the planet that the wave shocks. The explicit shape of this function is shown in Figure 3 of GR01, and its extension to a global disk is shown in Figure 3 of Rafikov (2002). The function f(r) can be described in terms of the parameter $\tau (r)$:

$$\begin{eqnarray}f(r)=\left\{\begin{array}{lc}{f}_{0} & \tau (r)\lt {\tau }_{\mathrm{sh}}\\ {f}_{0}\sqrt{{\tau }_{\mathrm{sh}}/\tau (r)} & \tau (r)\gt {\tau }_{\mathrm{sh}}\end{array}\right.,\end{eqnarray} \tag{ 13 }$$

where the shock position ${\tau }_{\mathrm{sh}}$ was calculated by Goodman & Rafikov (2001), and is given by

$$\begin{eqnarray}&&{\tau }_{\mathrm{sh}}=1.89+0.53/(q{{\mathcal{M}}}^{3}).\end{eqnarray} \tag{ 14 }$$

The parameter $\tau (r)$ represents an appropriately scaled distance from the planet. In shearing-box coordinates, this is given by

$$\begin{eqnarray}&&\tau (r)=\displaystyle \frac{{2}^{3/4}}{5}| \displaystyle \frac{3}{2}{\mathcal{M}}(r/a-1){| }^{5/2}.\end{eqnarray} \tag{ 15 }$$

To apply this to a global disk, in which surface density and sound speed vary as ${\rm{\Sigma }}\propto {r}^{-{p}_{1}}$ and $c\propto {r}^{-{p}_{2}}$, one must perform the following integral:

$$\begin{eqnarray}&&\tau (r)=\displaystyle \frac{3}{{2}^{5/4}}{{\mathcal{M}}}^{5/2}\left|{\displaystyle \int }_{1}^{r/a}| {s}^{3/2}-1{| }^{3/2}{s}^{(5{p}_{2}+{p}_{1})/2-11/4}{ds}\right|,\end{eqnarray} \tag{ 16 }$$

which can be integrated analytically and expressed in terms of hypergeometric functions. Numerical calculations presented in the next section will assume a uniform isothermal disk with ${p}_{1}={p}_{2}=0$.

Equations (12)–(14) and (16) completely specify the value of the surface density at every point in the disk, given the parameters q, α, and $h/r$, and the background surface density profile ${{\rm{\Sigma }}}_{0}(r)$.

The accuracy and range of applicability of this simple analytical model are tested by comparing directly with numerical calculations from the extensive parameter survey of Duffell (2015). In that survey, three parameters were varied: mass ratio q, Mach number ${\mathcal{M}}$ (or equivalently $h/r=1/{\mathcal{M}}$), and viscosity ν (or equivalently, $\alpha =\nu {{\mathcal{M}}}^{2}/({a}^{2}{{\rm{\Omega }}}_{p})$. This study used a uniform viscosity, $\nu \;=$ constant.

Figure 2 shows variation with mass ratio for a disk with ${\mathcal{M}}=20$ and $\alpha =0.01$. Mass ratios from $q=3.2\times {10}^{-5}$ to $q={10}^{-3}$ are studied showing a wide range of applicability of the 1D model. Near $q={10}^{-3}$ (Jupiter-mass), the model begins to break down as the gap is significantly wider than predicted, but the depth of the gap is still reasonably accurate. The non-smooth gap edges are due to the assumption that shock formation occurs at a single radius, i.e., that the angular momentum flux is described by a broken power-law (Equation (13)), with a sharp transition at the shock position.

The reason for the deviation around Jupiter-mass is likely that all of the excitation is assumed to come from within the gap, whereas once the gap becomes deep enough, the dominant resonances live just outside the gap edges. It is likely that a great deal of improvement would come from modifying the function ${{\rm{\Phi }}}_{p}(r)$ to account for excitation of waves in the presence of the detailed shape of the gap.

Figure 3 shows variation of the gap profile with viscosity, for a $q=6.4\times {10}^{-5}$ ($\sim 20{M}_{\oplus }$) planet in a disk with $h/r=0.05$ (${\mathcal{M}}=20$). The gap depth and shape is accurately captured for this range of viscosities.

Zoom In Zoom Out Reset image size

Finally, Figure 4 shows dependence on disk temperature (equivalently ${\mathcal{M}}$ or $h/r$) for $q=6.4\times {10}^{-5}$, $\alpha =0.01$ (α is fixed, but $\nu =\alpha {a}^{2}{{\rm{\Omega }}}_{p}/{{\mathcal{M}}}^{2}$ varies with ${\mathcal{M}}$). The shape of the gap is reasonably well-captured, but again once the gap becomes deep enough, there is a similar deviation from the model (this same deviation occurs at around $K\sim 100$ in Figure 1; note that the last panel in this figure has K = 42).

Zoom In Zoom Out Reset image size

This analytic model should be improved for the case of deep gaps (this is planned for a future study), but it accurately captures the "partial gap" regime, and the basic fundamentals are correct in that the scalings of Duffell & MacFadyen (2013) are reproduced, and low-mass planets can also open gaps in sufficiently low-viscosity disks, which has also been seen numerically (Dong et al. 2011a, 2011b; Duffell & MacFadyen 2012).

A simple analytical model has been presented for gap profiles in protoplanetary disks. In this model, effects due to planetary torque come directly from the detailed calculations of shock propagation by GR01 and Rafikov (2002).

The model is solved algebraically, and compared with long-duration numerical calculations. The model's performance in the "partial gap" regime is unprecedented in its agreement with numerical calculations and its reproduction of empirically found scalings. For deep gaps, the model does not correctly predict the gap shape. For example, Equation (8) does not extend to the inviscid limit, and it does not take into account effective viscosities generated by planet-induced instabilities. The failure of the model for deep gaps is most likely due to the fact that resonances located several scale heights from the planet become dominant once the gap suppresses the resonances nearest the planet. Additionally, effects due to non-Keplerian orbital motion should become important for very deep gaps. In practice the error that this assumption introduces is typically sub-dominant to the error in the assumed form of the planetary angular momentum flux.

For multiple gap-opening planets, the situation becomes more complicated, but still may be solvable. If the planets are sufficiently far apart, numerical calculations (Duffell & Dong 2015) have found that the gap profile is given by the product of the gap profiles from each planet by itself. However, if the planets are close enough that the gaps merge, the combined gap can be much shallower than the predicted depth, because each planet provides a time-dependent stirring to its neighbor's gap (Duffell & Dong 2015).

One could envision an iterative procedure for folding these effects into the analytic formula. For example, given the profile ${\rm{\Sigma }}(r)$, one could calculate the contribution of the torque from all resonances in the disk. This could give a new form for ${{\rm{\Phi }}}_{p}(r)$, which would provide a corrected surface density via Equation (11). This process could be repeated until a converged gap profile is found. Such a procedure will be attempted in a future study.

Resources supporting this work were provided by the NASA High-End Computing (HEC) Program through the NASA Advanced Supercomputing (NAS) Division at Ames Research Center. I am grateful to Eugene Chiang, Roman Rafikov, Andrew MacFadyen, Steve Stahler, Eliot Quataert, Kazuhiro Kanagawa, Alessandro Morbidelli, and Aurelien Crida for helpful comments. I would also like to thank the anonymous referee for a helpful review.