An empirically-derived formula for the shape of planet-induced gaps in protoplanetary disks

TL;DR

This paper derives an empirical formula for the shape of planet-induced gaps in protoplanetary disks using hydrodynamics simulations, enabling quick generation of synthetic disk profiles for various planetary and disk parameters.

Contribution

It introduces a novel empirical formula for the surface density profile of gaps in protoplanetary disks, based on extensive numerical simulations across key parameters.

Findings

Excellent agreement between the formula and numerical gaps across parameter space.

The formula accurately predicts gap depth and width for a wide range of planet masses.

A user-friendly code is provided for synthetic disk profile generation.

Abstract

This study uses numerical hydrodynamics calculations and a novel method for densely sampling parameter space to measure the precise shape of a gap opened by a planet in a gaseous disk, as a function of planet-to-star mass ratio, disk Mach number, and disk viscosity. Formulas for gap depth and width are determined, which are combined to form a complete formula for surface density as a function of radius in the disk. This new analytical formula is compared with numerically-derived gaps opened by planets ranging from very low masses up to a few times Jupiter's mass, and excellent agreement is found over a wide range of parameter space relevant to planet-disk interactions. A simple-to-use code is presented to rapidly generate synthetic disk profiles.

Equations

$${\partial }_t\Sigma +\nabla \cdot (\Sigma \vec{v})=0$$

$${\partial }_t(\Sigma v_j)+\nabla \cdot (\Sigma \vec{v}{v}_j+P{\hat{x}}_j-\nu \Sigma \vec{\nabla }{v}_j)=-\Sigma \vec{\nabla }\varphi$$

$$P=c^2\Sigma$$

$$q(t)=q_0{e}^{t/\tau },$$

$$\Sigma (R_{+}^d)=\Sigma (R_{-}^d)=d{\Sigma }_0.$$

$${\Sigma }^{\mathrm{DM}}/{\Sigma }_0=\frac{1}{1+\frac{0.45}{3\pi }\frac{q^2{\mathcal{M}}^5}{\alpha }}.$$

$$\Sigma /{\Sigma }_0=\frac{1}{1+\frac{0.45}{3\pi }\frac{q^2{\mathcal{M}}^5}{\alpha }\delta (q)}.$$

$$\delta (q)\equiv \frac{{\Sigma }_0/\Sigma -1}{0.45q^2{\mathcal{M}}^5/(3\pi \alpha )}$$

$$\delta(q)=\left\{\begin{array}[]{l@{\quad\quad}l}1&q<q_{\rm NL}\\ (q/q_{\rm NL})^{-1/2}+(q/q_{w})^{3}&q>q_{\rm NL}\\ \end{array}\right.$$

$$q_{\mathrm{NL}}=1.04{\mathcal{M}}^{-3},$$

$$q_w\sim {\alpha }^{1/5}/{\mathcal{M}}^{14/5},$$

$$q_w\sim q_{\mathrm{NL}}(\alpha \mathcal{M}{)}^{1/5}.$$

$$q_w=34q_{\mathrm{NL}}(\alpha \mathcal{M}{)}^{1/2}.$$

$${\Sigma }_{\mathrm{gap}}=0.015{\Sigma }_0{e}^{-600q},$$

$$((r/a{)}^{1/6}-1{)}^6=\frac{(q/q_{\Sigma }{)}^3-1}{D^3},$$

$$D(\mathcal{M},\alpha )=7{\mathcal{M}}^{3/2}/{\alpha }^{1/4}.$$

$$(\Delta /a{)}^6\sim {\left(\frac{q}{q_{\Sigma }D}\right)}^3.$$

$$q_{\Sigma }^2\sim \alpha /{\mathcal{M}}^5$$

$$\Delta /a\sim \frac{\sqrt{q}({\mathcal{M}}^5/\alpha {)}^{1/4}}{D^{1/2}}.$$

$$\Delta /a\sim \sqrt{q\mathcal{M}}/{\alpha }^{1/8}.$$

$$\Delta /a\sim \sqrt{q}{\mathcal{M}}^{3/4}/{\alpha }^{1/4}.$$

$$\Sigma (r)={\Sigma }_{\mathrm{gap}}(\tilde{q}),$$

$$\tilde{q}(r)\equiv \frac{q}{(1+D^3((r/a{)}^{1/6}-1{)}^6{)}^{1/3}}.$$

$$\Sigma (r)=\frac{{\Sigma }_0}{1+\frac{0.45}{3\pi }\frac{\tilde{q}(r{)}^2{\mathcal{M}}^5}{\alpha }\delta (\tilde{q}(r))}.$$

Full text

1. Introduction

When a sufficiently massive planet is embedded in a gaseous disk, that planet carves out a low-density annulus in the disk, also known as a “gap”. Gaps are important for a number of reasons; they might change the torques and accretion experienced by the planet (Ward, 1997; Kanagawa et al., 2018), but more importantly, gaps can be directly observed. In particular, the Atacama Large Millimeter Array (ALMA) has generated very detailed, resolved images of protoplanetary disks with distinct gap-like features (Huang et al., 2018a, b, c; Guzmán et al., 2018; Pérez et al., 2018). This provides a potential means of discovering exoplanets during the formation process.

This motivates a complete theory of gap-opening by planets; that is, given a planet mass, local disk temperature, and viscosity, what is the surface density of the gas as a function of radius?

One can address this question numerically, for a single set of disk and planet parameters (e.g. Duffell & MacFadyen, 2013; Fung et al., 2014), but much more useful would be a simple 1D analytical model, so that one can easily find best-fit disk and planet parameters for a given disk observation.

The idea of constructing a 1D gap model is not new; the first such models were constructed assuming torque is deposited in the disk at the same place it is excited (Varnière et al., 2004; Crida et al., 2006). These models did not always agree with numerically calculated gap profiles, especially in the “shallow gap” regime. In particular, they predicted that the planet needed to be above the thermal mass in order for it to open a gap, even in an inviscid disk. Additionally, they predicted an exponential scaling of gap depth with planet mass, which is also not seen in numerical studies.

Duffell & MacFadyen (2012) showed that low-mass planets can open gaps in sufficiently low-viscosity disks, and later found an empirical scaling for the gap depth as a function of planet mass, disk viscosity, and Mach number (Duffell & MacFadyen, 2013). Fung et al. (2014) confirmed this result, and presented a “zero-dimensional” interpretation for the scaling of Duffell & MacFadyen (2013). Kanagawa et al. (2015b) presented a disk model intended to improve upon past methods and agree with the scaling of Duffell & MacFadyen (2013).

Duffell (2015a) built a simple analytical disk model assuming all torque was deposited by shocks, using the semi-analytical scalings found by Goodman & Rafikov (2001) and Rafikov (2002) to determine where the shock forms and how much angular momentum it deposits at different locations in the disk. The gap model of Duffell (2015a) exhibited an unprecedented level of agreement between the model and numerical calculations for shallow gaps, but unfortunately the model was inaccurate for deep gaps.

Subsequently, disk models by Kanagawa et al. (2017) attempted to fix this model in the deep gap regime, also incorporating an empirical result for the scaling of gap widths (Kanagawa et al., 2015a). Ginzburg & Sari (2018) calculated analytical gap profiles for wider gaps, taking into account that torque can be excited at the gap walls.

In this work, we present an approach similarly empirical in nature to Kanagawa et al. (2017), but using a novel method to derive more information from the numerical calculations than any previous disk-planet study. Rather than keeping the planet mass fixed, or growing it self-consistently with some accretion prescription, the planet-to-star mass ratio is very slowly increased by hand from $q=10^{-6}$ to a few times $10^{-3}$. By increasing the mass ratio adiabatically, the planet always experiences an effectively steady-state disk, and when measuring the gap properties as a function of time, one effectively measures the gap properties as a continuous function of mass ratio. In other words, parameter space is densely sampled for all values of $q$ in this range.

The new calculations make it possible to continuously sample gap depth and width as a function of mass ratio, greatly increasing the ability to extract analytical scalings. From this, a new gap model is constructed, which agrees with numerical calculations for all planet masses up through a few times Jupiter’s mass.

The numerical method will be presented in Section 2, including an explanation of the diagnostics measured from the numerical calculations. Analysis of gap depths and widths will be presented in Section 3, and a 1D gap model will be constructed from this analysis. Finally, in Section 4 this analytical gap model will be compared with steady-state numerical calculations from previous studies. In Section 5 the results will be discussed, and future improvements suggested. Additionally, a simple-to-use code will be presented for rapidly generating synthetic disk profiles.

2. Methods

This study uses the moving-mesh code Disco (Duffell, 2016) to numerically integrate the equations of isothermal hydrodynamics in 2D:

[TABLE]

[TABLE]

[TABLE]

where $\Sigma$ is surface density, $P$ is pressure, $\vec{v}$ is velocity, $\nu$ is the kinematic viscosity, $c$ is the local sound speed, and $\phi$ is the gravitational potential from the planet and central star.

The numerical method and treatment of planetary potential is the same as described in previous studies using Disco (e.g. Duffell & MacFadyen, 2012; Duffell, 2015b). An important difference in the present study is that a novel approach is employed for the planetary potential: rather than choosing a single mass ratio and running the calculation until the disk is in steady-state, the mass ratio is very slowly increased from $q=10^{-6}$ to $q=2\times 10^{-3}$, in a sufficiently slow, adiabatic fashion such that at any single time the planet can be considered to be in a steady-state disk:

[TABLE]

with $q_{0}=10^{-6}$ and $\tau=10^{4}/\Omega$.

The benefit of this is that one can measure the gap depth and shape continuously as a function of mass ratio, rather than performing many separate numerical calculations to sample a discrete number of points. Measuring a continuous curve rather than discretely sampled points can be essential for measuring scalings, especially “break points” where the scalings change.

A happy by-product of this method is that it is much more well-behaved in a numerical sense; it has been previously shown that introducing the planetary potential too quickly can artificially excite vortices in the disk (Hammer et al., 2017).

The main diagnostics being measured are the gap depth, and the inner and outer radii at several different depths. The gap depth $\Sigma_{\rm gap}$ is simply the azimuthally-averaged surface density at the bottom of the gap. The contribution to this density due to high-density material very close to the planet is excised by only considering fluid elements in the half of the disk opposite the planet when computing this azimuthal average. The inner and outer radii of the gap $R_{\pm}^{0.9}$, $R_{\pm}^{0.5}$, and $R_{\pm}^{0.1}$ are also measured, defined by

[TABLE]

This will make it possible to measure the complete shape of the gap, using the position of the edges of the gap at different thresholds, for different values of $d$.

3. Analysis of Numerical Results

3.1. Gap Depth

Let us first consider gap depth as a function of mass ratio. Duffell & MacFadyen (2013) provided an empirical formula which was later explained via the analytical gap models of Kanagawa et al. (2015b) and Duffell (2015a):

[TABLE]

This scaling applies to low-mass planets, whose planetary torques can be reasonably approximated in linear theory. Deviations from the scaling of Duffell & MacFadyen (2013) begin at mass ratio of a few times $10^{-4}$ (See Figure 1). This can be most easily demonstrated by defining the function $\delta(q)$ such that

[TABLE]

In other words, $\delta(q)$ is defined as

[TABLE]

The reason for invoking this particular expression is that $\delta(q)$ should be unity as long as the Duffell & MacFadyen (2013) scaling applies, so one can simply measure the break point where it deviates from unity, and find the new scaling after that point. This is shown in Figure 1, and again in Figure 2 for several choices of disk mass ratio and viscosity, to demonstrate the scaling of the break point with $\mathcal{M}$ and $\alpha$.

If one can develop an empirical formula for $\delta(q)$, then one immediately has a complete formula for gap depth as a function of disk and planet parameters. Figure 1 shows the deviation of $\Sigma_{\rm gap}$ from the scaling of Duffell & MacFadyen (2013). This deviation occurs at a particular mass ratio which will be called “$q_{\rm NL}$” as it will be shown to coincide with the threshold for strong nonlinearity in the spiral wave.

The lower panel of Figure 1 plots $\delta(q)$, which shows this deviation more clearly. It is possible to fit $\delta(q)$ with some power-laws; this gives the function

[TABLE]

To determine $\delta(q)$ generally, we must first see how this scaling changes when $\mathcal{M}$ and $\alpha$ are varied.

Note that when $q$ becomes sufficiently large, the $\Sigma_{\rm gap}(q)$ scaling changes again, but the curve becomes much more erratic. At this point, instabilities have been excited in the disk by the planet. This will not be taken into account in our formula for $\delta(q)$ but a correction will be suggested later.

Figure 2 shows how $\delta(q)$ depends on $\mathcal{M}$ and $\alpha$. The break point $q_{\rm NL}$ very simply scales as

[TABLE]

with no dependence on $\alpha$.

This is consistent with the picture that the break point occurs when the spiral wave becomes strongly nonlinear; $q\sim\mathcal{M}^{-3}$ is also the correct scaling for the thermal mass, or the mass ratio at which the Hill radius becomes larger than the disk scale height.

Apparently this non-linearity causes the torque deposited by shocks to be somewhat weaker than would be predicted by analytical theory. This was also noted by Kanagawa et al. (2017), who suggested a nonlinear adjustment to the Duffell & MacFadyen (2013) scaling by a factor of $0.4$ (referred to as $f_{\rm NL}$ in the Kanagawa paper). However, it was not clear at what point their adjustment should be applied.

Additionally, even with this adjustment it is clear that the scaling will be incorrect as the gap deepens; the $q^{3}$ term in $\delta(q)$ might be due to the fact that the wave is no longer being predominantly excited at the bottom of the gap, but at the higher-density gap walls (Ginzburg & Sari, 2018). Regardless of the reason, this effect becomes important around $q\sim q_{w}$, whose scaling we can also determine empirically.

Ginzburg & Sari (2018) predicted that the gap depth scaling should change due to excitation of torque at the gap walls. This transition was predicted to occur when

[TABLE]

or, in terms of $q_{\rm NL}$,

[TABLE]

In practice, instead of a $1/5$ power, a steeper scaling is found:

[TABLE]

Interestingly, this value of $q_{w}$ is the mass ratio at which the gap attains the depth $\Sigma_{\rm gap}\approx 0.02\Sigma_{0}$ according to the Duffell & MacFadyen (2013) formula. So it makes some sense that excitation at the gap walls would dominate when the gap becomes sufficiently deep.

Additionally, Ginzburg & Sari (2018) predicted the scaling $\delta(q)\sim(q/q_{w})^{2}$ for this regime, which is shallower than the scaling observed here, hence the power of three in Equation (9).

The present study only investigated disks for which $q_{\rm NL}<q_{w}$. That is, for $\alpha>8.7\times 10^{-4}/\mathcal{M}$. For sufficiently low viscosity, excitation of torque at the gap walls becomes important before nonlinearity in the wave affects the scaling. It is possible that the Ginzburg & Sari (2018) scalings become valid in this low-viscosity regime, as that is the regime in which they were calculated.

Finally, with this formula for $\delta(q)$, it is now possible to plot $\Sigma_{\rm gap}(q,\alpha,\mathcal{M})$, comparing numerical results to the new analytical formula, for all mass ratios and disk models considered. This is shown in Figure 3. This figure also shows how this scaling breaks down when the gap becomes too deep and unstable. This regime will not be characterized in this study, except to note that it appears possible to describe the unstable regime as hitting a floor in the density of

[TABLE]

indicated by the red curve in the figure. This “floor” did not work for all models (curiously, the lowest viscosity model in the lower panel achieved a much deeper gap) but it appears to work for many disks.

3.2. Gap Width

Gap width in this study is measured at three different surface densities: $\Sigma/\Sigma_{0}=0.1$, $0.5$, and $0.9$. The inner and outer radii with this surface density are measured and plotted in the top panel of Figure 4 as a function of mass ratio for a particular disk model $\mathcal{M}=20$, $\alpha=0.01$.

The gap is not symmetric; $R_{+}$ and $R_{-}$ are not the same distance from the planet. However, through experimentation it was found that $(R_{+})^{1/6}$ and $(R_{-})^{1/6}$ are roughly equidistant from $a^{1/6}$. The middle panel of Figure 4 shows this, by plotting $|(R_{\pm}/a)^{1/6}-1|$ for a range of planet masses and for the three different gap depths measured.

After viewing the central panel of Figure 4, it becomes clear that the three curves corresponding to gap widths at different depths appear to obey a very similar scaling to one another, and that they should approximately lie on top of one another if one re-scales the mass ratio by $q_{\Sigma}$, the mass ratio at which the gap attains the depth $\Sigma$.

This final curve is fit by a power-law. For $\mathcal{M}=20$ and $\alpha=0.01$, this curve is found to be well-fit by the function

[TABLE]

where $D=2000$ for this disk.

This provides an (apparently) universal scaling for the gap width with mass ratio. However, we have yet to establish how the quantity $D$ varies with $\mathcal{M}$ and $\alpha$.

In practice, we find this by fitting the same curves for different disk models (see Figure 5). The scaling found was

[TABLE]

One can compare this to the measured scaling found by Kanagawa et al. (2015a). Note, in the regime where the width $\Delta\ll a$ but $q\gg q_{\Sigma}$, equation (15) becomes:

[TABLE]

Now, apply the Duffell & MacFadyen (2012) scaling for $q_{\Sigma}$, i.e. the scaling of $q$ with $\alpha$ and $\mathcal{M}$ to fix a specific gap depth $\Sigma$:

[TABLE]

and therefore,

[TABLE]

Using the measured scaling for $D$, one arrives at:

[TABLE]

Comparing this to the Kanagawa et al. (2015a) scaling:

[TABLE]

The scalings are reasonably close, off by a factor of $(\alpha/\mathcal{M}^{2})^{1/8}$. Any discrepancies could potentially be explained by the fact that the gap width does not really follow a power-law.

Figure 5 shows that our scaling is close to the numerical curves $R_{\pm}(q,\alpha,\mathcal{M})$ for all disk models considered. It is worth noting that this scaling fails once the disk becomes unstable and turbulent; the gap widens significantly in that case. It is possible that this is related to the different scalings found by Kanagawa et al. (2015a). A future study can address this discrepancy.

4. Comparing the Model to Steady-State Disks

Equipped with the scaling (15), along with the formula for gap depth (7), it becomes possible to build a complete formula for $\Sigma(r)$.

First, note that “$q_{\Sigma}$” is just the inverse function of $\Sigma_{\rm gap}(q)$. Therefore, one can simply write

[TABLE]

where $\tilde{q}$ is given by solving equation (15) for $q_{\Sigma}$:

[TABLE]

Then,

[TABLE]

Equations (23) and (24), along with equation (9) for $\delta(q)$ and equation (16) for $D$, provide a complete formula for $\Sigma(r)$, given any $q$, $\alpha$, and $\mathcal{M}$.

We can put this model to the test by comparing with steady-state numerical calculations of gap-opening planets. For this, we use steady-state models generated from the disk-planet study of Duffell (2015b) (Figure 6). Excellent agreement is found between the numerical and analytical gap profiles for a wide range of gap-opening planets. In particular, the model is able to capture both shallow and deep gaps, whereas the analytical model of Duffell (2015a) was only able to capture shallow gaps.

Figure 7 shows how the new formula for gap depths compares with other formulas in the literature. By design, the model agrees with Duffell & MacFadyen (2013) in the low-mass regime, but deviates in the high-mass regime (as also shown in Figure 1). The scaling measured empirically by Fung et al. (2014) appears to precisely recover the tangent line near Jupiter’s mass. The scaling of Kanagawa et al. (2017) with the nonlinear correction $f_{\rm NL}=0.4$ reasonably captures gap depths up to Jupiter’s mass, but fails above this mass. The scaling here is the only one to recover the gap depth in all regimes, up to a few times Jupiter’s mass (after which point instabilities in the disk may take over, depending on the disk properties).

5. Discussion

The calculations presented here can be used to model gaps in resolved images of protoplanetary disks. The formula for $\Sigma(r)$ presented here, while being largely empirical in nature, may help guide analytical studies to come up with an explanation for these formulas. Understanding the various break points in $\delta(q)$, for example, will require additional analytical work, which can be tested against this empirical formula rather than running full numerical hydrodynamics calculations.

In order to make these scalings as easy-to-use as possible, this 1D gap model has been built into a short, simple code to rapidly generate synthetic disk profiles. This code is available for download at https://github.com/duffell/Duffell_Gap. The code is written in both C and Python and is freely available under the MIT license.

The scalings found here break down when the density gradients become steep enough to trigger instabilities in the disk, leading to the production of vortices. It was found that the onset of instability was suppressed by extending the boundary further from the planet and increasing the planet mass slowly. More care must be taken to determine whether the numerical calculation is faithfully capturing the instability. Future work will require a meticulous survey to characterize the unstable regime. The analytical model is accurate enough in the stable regime that it may be possible to use it to analytically calculate the point where the disk should become unstable to different instabilities such as the Rayleigh instability or Rossby wave instability. This idea is left for a future study.

Another obvious deviation between the empirical model and the numerical results is the “tails” in the gap located many scale heights away from the planet (see for example the top panel of Figure 6, at $r=0.5a$ and $r=1.5a$). These deviations may be reasonable to model empirically as well (in fact, a very similar feature was predicted by Ginzburg & Sari, 2018), but modeling this structure was not attempted in this study because of a lack of confidence in the numerical results for these “tails”; they may be affected by the inner and outer boundary. This detail may be resolvable in a future study.

Finally, a very clear feature that occurs for many of the shallow gaps in Figure 6 is an asymmetry in the gap; the outer side of the gap is consistently more depleted than the inner side. The likely cause of this is the same torque imbalance that leads to Type I migration; the outer Lindblad resonances consistently excite stronger torques than the inner Lindblad resonances. This could potentially used as an indicator to observers as to whether the gap is caused by a planet; if a gap were observed whose asymmetry went in the opposite direction, it might be difficult to model the gap using a planet. However, modeling this asymmetry was not attempted in this study; it appears to occur mainly for shallow gaps where $\Sigma_{\rm gap}\sim 0.5\Sigma_{0}$.

I am grateful to Karin Öberg, Jane Huang, Sivan Ginzburg, Andrew MacFadyen, Daniel D’Orazio, Zoltan Haiman and Andrea Derdzinski for helpful comments and discussions.