Multiple Disk Gaps and Rings Generated by a Single Super-Earth

TL;DR

This study demonstrates that a single super-Earth can create multiple disk gaps and rings in protoplanetary disks, which are observable with ALMA, providing insights into planet formation processes.

Contribution

The paper introduces a novel simulation approach showing how a super-Earth induces multiple gaps and rings in gas and dust distributions, matching observed features in disks like HL Tau.

Findings

Super-Earths open double gaps in low-viscosity disks.

Dust accumulates into rings at gap edges and between gaps.

Double gaps are detectable with ALMA at high resolution.

Abstract

We investigate the observational signatures of super-Earths (i.e., Earth-to-Neptune mass planets) in their natal disks of gas and dust. Combining two-fluid global hydrodynamics simulations with a radiative transfer code, we calculate the distributions of gas and of sub-mm-sized dust in a disk perturbed by a super-Earth, synthesizing images in near-infrared scattered light and the mm-wave thermal continuum for direct comparison with observations. In low viscosity gas ($\alpha\lesssim10^{-4}$), a super-Earth opens two annular gaps to either side of its orbit by the action of Lindblad torques. This double gap and its associated gas pressure gradients cause dust particles to be dragged by gas into three rings: one ring sandwiched between the two gaps, and two rings located at the gap edges farthest from the planet. Depending on system parameters, additional rings may manifest for a single planet. A double gap located at tens of AUs from a host star in Taurus can be detected in the dust continuum by the Atacama Large Millimeter Array (ALMA) at an angular resolution of ~0".03 after two hours of integration. Ring and gap features persist in a variety of background disk profiles, last for thousands of orbits, and change their relative positions and dimensions depending on the speed and direction of planet migration. Candidate double gaps have been observed by ALMA in systems like HL Tau (D5 and D6) and TW Hya (at 37 and 43 AU); we submit that each double gap is carved by one super-Earth in nearly inviscid gas.

Figures

Tables

Table 1. Table 1 Models. In all non-migration simulations, the planet’s orbit is fixed at r p = 30 subscript 𝑟 p 30 r_{\rm p}=30 AU, where h / r = 0.05 ℎ 𝑟 0.05 h/r=0.05 ( M thermal = 40 ​ M ⊕ subscript 𝑀 thermal 40 subscript 𝑀 direct-sum M_{\rm thermal}=40M_{\oplus} ). In the last four models with a migrating planet, the migration rate is fixed at | r ˙ | = 2 × 10 − 3 ˙ 𝑟 2 superscript 10 3 |\dot{r}|=2\times 10^{-3} AU orbit -1 (or about 10 − 5 superscript 10 5 10^{-5} AU year -1 ); the planet reaches r p = subscript 𝑟 p absent r_{\rm p}= 30 AU after 0.25 Myr.

Model Name	$M_p$	$\alpha$	Initial ${\Sigma }_0-r$	Migration
2ME	$2M_{\oplus }$	$5\times {10}^{-5}$	$\propto 1/r$	No
5ME	$5M_{\oplus }$	$5\times {10}^{-5}$	$\propto 1/r$	No
10ME	$10M_{\oplus }$	$5\times {10}^{-5}$	$\propto 1/r$	No
17ME	$17M_{\oplus }$ (Neptune)	$5\times {10}^{-5}$	$\propto 1/r$	No
60ME	$60M_{\oplus }$	$5\times {10}^{-5}$	$\propto 1/r$	No
10ME-Beta0	$10M_{\oplus }$	$5\times {10}^{-5}$	Constant	No
10ME-Beta15	$10M_{\oplus }$	$5\times {10}^{-5}$	$\propto 1/r^{1.5}$	No
17ME-LowVis	$17M_{\oplus }$	$5\times {10}^{-6}$	$\propto 1/r$	No
17ME-HighVis	$17M_{\oplus }$	$5\times {10}^{-4}$	$\propto 1/r$	No
10ME-Beta0-MigrateIn	$10M_{\oplus }$	$5\times {10}^{-5}$	Constant	$\dot{r}=-2\times {10}^{-3}$AU orbit-1
10ME-Beta0-MigrateOut	$10M_{\oplus }$	$5\times {10}^{-5}$	Constant	$\dot{r}=2\times {10}^{-3}$AU orbit-1
17ME-MigrateIn	$17M_{\oplus }$	$5\times {10}^{-5}$	$\propto 1/r$	$\dot{r}=-2\times {10}^{-3}$AU orbit-1
17ME-MigrateOut	$17M_{\oplus }$	$5\times {10}^{-5}$	$\propto 1/r$	$\dot{r}=2\times {10}^{-3}$AU orbit-1

Equations

$$M_{\mathrm{thermal}}\simeq M_{\star }{\left(\frac{h}{r}\right)}^3\simeq 40M_{\oplus }\left(\frac{M_{\star }}{M_{\odot }}\right){\left(\frac{h/r}{0.05}\right)}^3$$

$$l_{\mathrm{sh}}\approx 0.8{\left(\frac{\gamma +1}{12/5}\frac{M_{\mathrm{p}}}{M_{\mathrm{thermal}}}\right)}^{-2/5}h$$

$$\frac{{\Sigma }_{\mathrm{gap}}}{{\Sigma }_0}\simeq 0.4{\left(\frac{3\times 10^{-5}}{M_{\mathrm{p}}/M_{\star }}\right)}^{2.16}{\left(\frac{h/r}{0.05}\right)}^{6.61}{\left(\frac{\alpha }{10^{-4}}\right)}^{1.41}$$

$${\Sigma }_{\mathrm{gas},0}=5.56{\left(\frac{r_{\mathrm{p}}}{r}\right)}^{\beta }\mathrm{g}\mathrm{c}{\mathrm{m}}^{-2},$$

$$\frac{h}{r}=0.05{\left(\frac{r}{r_{\mathrm{p}}}\right)}^{0.25};$$

$$\mathrm{St}=\frac{\pi s{\rho }_{\mathrm{dust}}}{2{\Sigma }_{\mathrm{gas}}},$$

$${\rho }_{\mathrm{small}-\mathrm{dust}}(z)$$

$${\rho }_{\mathrm{big}-\mathrm{dust}}(z)$$

$$\frac{\partial {\Sigma }_{\mathrm{dust}}}{\partial t}+\frac{1}{r}\frac{\partial (r{\Sigma }_{\mathrm{dust}}{v}_{r,\mathrm{dust}})}{\partial r}\approx 0,$$

$$v_{r,\mathrm{dust}}=v_{r,\mathrm{gas}}+\frac{t_{\mathrm{s}}}{{\Sigma }_{\mathrm{gas}}}\frac{\partial p}{\partial r}$$

Full text

1. Introduction

Practically every star in our Galaxy harbors planets (e.g., Fressin et al., 2013; Burke et al., 2015; Dressing & Charbonneau, 2015). The most common planets discovered so far are “super-Earths,” having radii $R_{\rm p}$ intermediate between those of Earth and Neptune (for reference, Neptune’s $R_{\rm p}=4R_{\oplus}$ and its mass $M_{\rm p}=17M_{\oplus}$).111Planets at the upper end of this size range, $2R_{\oplus}\leq R_{\rm p}\leq 4R_{\oplus}$, are sometimes called “sub-Neptunes” or “mini-Neptunes”; here we use the term “super-Earth.” The Kepler mission showed that more than half of FGKM dwarfs have super-Earth companions out to orbital periods of $\sim$100 days (e.g., Fressin et al., 2013; Dressing & Charbonneau, 2015). These planets generally weigh less than $20M_{\oplus}$, as determined from transit-timing analyses (e.g., Jontof-Hutter et al., 2014; Wu & Lithwick, 2013) and radial velocity measurements (e.g., Weiss & Marcy, 2014; Marcy et al., 2014; Dumusque et al., 2014; Dressing et al., 2015). Their compositions, as inferred from their bulk densities, span a wide range; some are consistent with pure rock while others have 1–10% of their mass in voluminous H/He envelopes (e.g., Lopez & Fortney, 2014; Rogers, 2015; Lee & Chiang, 2016).

Nearly all such planets are detected orbiting Gyr-old host stars. To better understand how super-Earths form, we need to find them in younger systems, while they are still emerging from their parent circumstellar disks. Unfortunately, the two most prolific planet detection methods, the radial velocity and transit techniques, are not suited to detecting super-Earths in disks. Pre-main-sequence stars tend to be too spectroscopically and photometrically noisy because of magnetic activity and ongoing, time-variable accretion. An exception is provided by K2-33b, an $R_{\rm p}=6R_{\oplus}$ planet discovered transiting a 5–10 Myr old star (David et al., 2016). This object may be a super-Earth that will eventually cool and contract to a radius $R_{\rm p}<4R_{\oplus}$. Radial velocity measurements of its mass are not, however, strongly constraining ($M_{\rm p}<3.6M_{\rm J}$, where $M_{\rm J}=318M_{\oplus}$ is the mass of Jupiter). Another example of a young planet, discovered using radial velocities, is V830 Tau b, but it is a comparatively rare gas giant, not a super-Earth ($M_{\rm p}=0.77\pm 0.15M_{\rm J}$; Donati et al. 2016). As far as direct imaging goes, gas giants a few dozen Myr old have been detected with high-contrast near-infrared (NIR) adaptive optics systems (Gemini/GPI and VLT/SPHERE): see 51 Eri b (Macintosh et al., 2015) and HD 131399Ab (Wagner et al., 2016). The current detection limit of such instruments is $\sim$1$M_{\rm J}$ for planets a few Myr old (e.g., Bowler, 2016). Directly imaging super-Earths is still out of reach.

On the other hand, these same NIR adaptive optics technologies, in combination with the Atacama Large Millimeter Array (ALMA), have enabled us to detect structures in disks produced by gravitational interactions with embedded planets, and thus to observe, indirectly, planet formation in action. Scattered light images in the NIR and thermal emission maps at millimeter (mm) wavelengths have achieved angular resolutions less than $0\arcsec.05$, corresponding to length scales of a few AUs at the locations of nearby star-forming regions such as Taurus. Disk-planet interactions can produce a variety of structures, depending on the planet mass $M_{\rm p}$, the disk scale height $h$, and the disk viscosity as parameterized by $\alpha$, the ratio of the viscous stress to the thermal pressure (Shakura & Sunyaev, 1973). So far, three types of possible planet-induced structures have been possibly detected in disks:

Gaps (in, e.g., HD 169142; Quanz et al. 2013; Momose et al. 2015; Fedele et al. 2017). 2. 2.

Density waves (in, e.g., AB Aur; Hashimoto et al. 2011; Tang et al. 2017). 3. 3.

Vortices created by the Rossby wave instability (RWI; Lovelace et al. 1999; Li et al. 2001, 2005) at the edges of planet-opened gaps (in, e.g., IRS 48 and HD 142527; van der Marel et al. 2013; Casassus et al. 2013).

In this paper, we expand upon the first of these disk features, showing how a single planet can actually carve out two gaps to either side of itself, separated by a lane of co-orbital material. This is the signature of a low mass planet in a low viscosity or inviscid disk. By “low mass,” we mean a planet whose mass is less than the thermal mass

[TABLE]

where $M_{\star}$ is the host stellar mass and $r$ is the planet-star separation. In flared disks, $h/r$ increases with $r$, and by extension so does $M_{\rm thermal}$. Thermal mass planets excite density waves at their Lindblad resonances that are already non-linear (characterized by fractional surface density perturbations $\Delta\Sigma/\Sigma\sim 1$) at launch. Sub-thermal masses—a category to which super-Earths belong—generate weaker waves. By “low viscosity,” we mean that the waves damp less by intrinsic disk viscosity and more by steepening and eventually breaking. Goodman & Rafikov (2001) showed that in an inviscid disk having a two-dimensional (2D) adiabatic index $\gamma$, waves launched to either side of the planet’s orbit travel a radial distance

[TABLE]

before shocking and dissipating (see also Rafikov, 2002a). For planets with $M_{\rm p}/M_{\rm thermal}\sim 0.1$–$0.5$ (i.e., $M_{\rm p}=4$–$20M_{\oplus}$ with $h/r=0.05$), we have $l_{\rm sh}\sim 1$–$2h$. Where the waves break, they deposit the angular momentum they carry to disk gas. As a result, material is repelled away from the planet to form two gaps, at radial separations $\pm l_{\rm sh}$ from the planet’s orbit. Numerical simulations confirm this “double gap” feature in inviscid disks (see, e.g., Fig. 2 in Dong et al. 2011a or Fig. 8 in Duffell & MacFadyen 2012; see also Li et al. 2009; Muto et al. 2010; Dong et al. 2011b; Zhu et al. 2013).

The double gaps are expected to be shallow. From the numerical fitting formula of Fung et al. (2014; see also Duffell & MacFadyen 2013; Kanagawa et al. 2015; Duffell 2015), a super-Earth of mass $10M_{\oplus}$ orbiting a $1M_{\odot}$ star opens a gap having a surface density contrast of

[TABLE]

where $\Sigma_{\rm gap}$ is the surface density inside the gap and $\Sigma_{0}$ is the unperturbed value. The above formula is known to break down as $\alpha\rightarrow 0$ because hydrodynamic instabilities substitute for viscous diffusion in filling up disk gaps; simulations by Fung & Chiang (2017) show that super-Earths open gaps with $\Sigma_{\rm gap}/\Sigma_{0}\sim 0.1$ when $\alpha=0$ and $h/r=0.03$.

Although gas surface density contrasts generated by super-Earths may be modest, aerodynamic interactions between gas and dust can produce higher contrast features in dust (e.g., Rice et al., 2006; Pinilla et al., 2012b). In particular, the edges of planet-opened gaps are regions of high gas pressure, and dust particles drift toward such regions and collect into overdense rings (e.g., Paardekooper & Mellema, 2004, 2006; Zhu et al., 2012; Pinilla et al., 2012a). Dust rings have been detected (e.g., Zhang et al., 2014; Hashimoto et al., 2015; van der Marel et al., 2015; Dong et al., 2017).

In this paper, by combining global 2D two-fluid (gas+dust) hydrodynamical simulations with a radiative transfer code, we study how super-Earths produce characteristic double-gap patterns in mm-sized dust particles in nearly inviscid ($\alpha\lesssim 10^{-4}$) disks. Such structures should be readily detectable in ALMA dust continuum observations of disks 140 pc away with $\sim$0$\arcsec.03$ angular resolution. We will introduce our simulations in Section 2, present results in Section 3, and summarize and discuss in Section 4.

2. Simulations

We calculate the dust and gas distributions in a disk perturbed by a planet using global 2D (radial-azimuthal) hydrodynamics simulations with the LA-COMPASS code (Li et al. 2005, 2009; Fu et al. 2014; Section 2.1). Hydro models are post-processed using the Whitney et al. (2013) Monte Carlo radiative transfer (MCRT) code to produce synthetic images at various wavelengths (Section 2.2).

2.1. Hydrodynamics Simulations

The simulation setup is largely adopted from Fu et al. (2014) and Miranda et al. (2016), and is briefly summarized here. Dust is treated as a pressureless fluid (e.g., Takeuchi & Lin, 2002). An $\alpha$-viscosity is added to the gas dynamics, which introduces an additional diffusion term in the dust continuity equation (see Eqn. 28 in Takeuchi & Lin, 2002). We do not model turbulent mixing. An axisymmetric 2D global disk is initialized in cylindrical coordinates ($r\times\phi=1024\times 1024$) with a gas surface density

[TABLE]

where $\beta$ is the global radial power law index ranging from 0 to 1.5, and $r_{\rm p}=30$ AU is the orbital radius of the planet. The planet’s orbit is coplanar with the disk, circular, and fixed in time for our default simulations; in Section 3.5, we experiment with prescriptions for radial migration. The planet mass ramps up with time $t$ as $M_{\rm p}(t)\propto\sin{(t)}$ over $0\leq t\leq\tau_{\rm growth}$. Our default setting is $\tau_{\rm growth}=10$ orbits; in Section 3.1.4, we experiment with $\tau_{\rm growth}=1000$ and 3000 orbits. The size of the simulation domain is 2$\pi$ in $\phi$ and 0.1–2.1$r_{\rm p}$ in $r$ (3–64 AU). At the end of the hydro simulations, before MCRT post-processing, the inner disk within twice the inner boundary is manually removed to avoid possible boundary condition artifacts. Most models take $\beta=1$, for which the initial gas mass is about $0.01M_{\odot}$ within 60 AU. The gas surface density and velocity are fixed at their initial values at both the inner and outer boundaries.

The equation of state of the gas is locally isothermal. The scale height of the gas disk $h$ is taken as

[TABLE]

the disk is flared. The central star is $1M_{\odot}$ in mass, and the thermal mass at $r_{\rm p}$ is $0.05^{3}\times M_{\odot}=40M_{\oplus}$ (Eqn. 1). The gas has a kinematic viscosity $\nu=\alpha h^{2}\Omega_{\rm k}$ (Shakura & Sunyaev, 1973) with constant $\alpha$ and $\Omega_{\rm k}$ equal to the Keplerian angular velocity. We experiment with $\alpha$ between $5\times 10^{-4}$ and $5\times 10^{-6}$.

Dust particles with size $s$ are treated as a second fluid in the hydro simulation. The Stokes number St (momentum stopping time normalized to the dynamical time) of dust particles in the Epstein regime is

[TABLE]

where $\rho_{\rm dust}$ is the internal density of a dust particle. At $t=0$ and $r=r_{\rm p}$, ${\rm St}=0.007$ for $\rho_{\rm dust}=1.2$ g cm*-3* and $s=0.2$ mm (“big dust”; see Section 2.2). The two-fluid hydrodynamic equations describing the coupled evolution of gas and dust are solved in LA-COMPASS including both aerodynamic drag on dust, and dust back-reaction on gas. The big dust particles are initialized with a surface density linearly proportional to the gas, $\Sigma_{\rm big-dust,0}=0.5\%\times\Sigma_{\rm gas,0}$. An outflow inner boundary condition and an inflow outer boundary condition are imposed on the dust.

We summarize the models and their parameters in Table 1. Our default is to run for 1500 orbits (0.25 Myr at 30 AU); in Section 3.1.3, we experiment with longer run times.

2.2. Radiative Transfer Simulations and Synthetic Observations

We follow the procedures described in Dong et al. (2015) to translate the hydro disk models to synthetic observations. The main steps are summarized here. Light from the central star is scattered and absorbed by “small dust” particles at the disk surface. The small dust is assumed to resemble sub-$\mu$m-sized, interstellar medium (ISM) dust in its optical properties (Kim et al. 1994; see Fig. 2 in Dong et al. 2012). Because of their short stopping times, small dust grains are assumed to be always well mixed with gas at a mass fraction of 0.5%. Big dust particles, having sizes $s=0.2$ mm as in the hydro simulations, are assumed to have settled to the disk midplane (e.g., Dullemond & Dominik, 2004; Birnstiel et al., 2010). Initially in the hydro simulations, the big-to-small-dust mass ratio is 1:1, so that the initial total-dust-to-gas mass ratio is 1:100; at the end of the hydro simulation, the big-dust-to-gas ratio in a given gas parcel has evolved because of drag forces. The optical properties of the big dust are calculated using the Bohren & Huffman (1983) Mie routine,222The B. T. Draine version, https://www.astro.princeton.edu/$\sim$draine/scattering.html, wrapped into a Python package by T. Robitaille, https://github.com/hyperion-rt/bhmie. assuming an ISM composition. The Whitney MCRT code calculates the equilibrium temperature of all dust particles using iterative methods (Lucy, 1999).

We assume the central source is a luminous blackbody sphere of radius $R_{\star}=2.3R_{\odot}$ and $T=4350$K, appropriate for a 1 $M_{\odot}$ star at 1 Myr (Baraffe et al., 1998). We manually fill the inner disk by using the initial surface density profile (Eqn. 4) to extrapolate surface densities inward from $r=0.2r_{\rm p}$. The hydro simulations are 2D and yield only the surface densities $\Sigma_{\rm small-dust}$ and $\Sigma_{\rm big-dust}$; these are made into 3D density fields for use in the MCRT calculations by puffing them up assuming Gaussian profiles in the vertical direction $z$:

[TABLE]

where $h_{\rm small-dust}$ is set equal to the gas scale height $h$ (Eqn. 5), and $h_{\rm big-dust}=0.1h_{\rm gas}$ to mimic dust settling, as in Miranda et al. (2016).

We produce both $H$-band polarized intensity scattered light images and ALMA band 7 (870 $\micron$; 345 GHz; 7.5 GHz bandwidth) continuum observations, assuming the disk is 140 pc away at Dec=+25∘ (e.g, Taurus). Unless otherwise noted, we assume 4 hours of integration time for simulated ALMA observations. For $H$-band images, we take an inner working angle of $0\arcsec.07$ (10 AU) and convolve the raw MCRT products with a Gaussian point spread function (PSF) having a full-width-half-maximum of $0\arcsec.04$ (6 AU) to simulate the performances of VLT/SPHERE and Gemini/GPI. For mm dust continuum observations, we transform the raw MCRT products to simulated ALMA observations using the simobserve and simanalyze tools in Common Astronomy Software Applications (CASA). A full array of 50 12-m antennas is used. We mostly use the array configuration alma.out24.cfg,333https://casaguides.nrao.edu/index.php?title=Antenna_Configurations_Models_in_CASA. achieving a beam size of $0\arcsec.029\times 0\arcsec.020$ with the major axis of the beam in the north-south direction (by comparison, the smallest beam achievable for a Taurus target at 345 GHz with the full ALMA array is $0\arcsec.017\times 0\arcsec.012$ with configuration alma.out28). Simulated ALMA observations are corrupted by thermal noise assuming the default atmospheric model with 0.9 mm precipitable water vapor ($3^{\rm rd}$ octile). Four hours of integration reaches a root-mean-squared (RMS) sensitivity of 13 $\mu$Jy beam*-1* using the ALMA sensitivity calculator.444https://almascience.eso.org/proposing/sensitivity-calculator Throughout the paper, we use a gray color scheme for surface density, a blue-hot color scheme in the $H$-band, and a red-hot color scheme for ALMA images. All scattered light images are for face-on disks, while for ALMA images we experiment with non-zero disk inclinations.

3. Results

After presenting a fiducial model that illustrates the basic observable disk features produced by a super-Earth (Section 3.1), we then explore how these features vary with planet mass (Section 3.2), disk viscosity (Section 3.3), background surface density profile (Section 3.4), and orbital migration of the planet (Section 3.5).

3.1. A 10-$M_{\oplus}$ Super-Earth at 30 AU

We use model 10ME (see Table 1) to showcase the basic disk features produced by a super-Earth in a low viscosity environment.

3.1.1 Surface Density Distributions

Figure 1 shows the surface density maps of the gas and of the $0.2$-mm-sized big dust, as well as their azimuthally averaged radial profiles. Henceforth, for convenience, we sometimes omit “big” when referring to “big dust.” In the gas, the 10-$M_{\oplus}$ super-Earth gradually opens two narrow gaps about $2h$ away from its orbit, as predicted by Rafikov (2002b) in this nearly inviscid environment. The gas in the coorbital region is less depleted than in the two gaps, and forms an annular ring. We identify six features in the gas radial profile: apart from the IG (inner gap), OG (outer gap), and MR (middle ring), we observe an IR (inner ring), an OR (outer ring), an IR2 (inner ring 2), and an IG2 (inner gap 2). The gas originally in the gap regions is pushed away and piled up at the gap edges to form the IR and OR. Because of the low viscosity, viscous diffusion is unable to smooth away the two rings at the time shown ($t=1500$ orbits). These structures in the gas surface density have been seen in previous hydro simulations of low mass planets in low viscosity disks (see, e.g., Figure 2 of Li et al. 2009, Figure 14 of Muto et al. 2010, and Figure 8 of Duffell & MacFadyen 2012). The IR2 feature (and an accompanying annular gas depletion IG2) has been seen in simulations with inner boundaries smaller than $\sim$0.4$r_{\rm p}$ (see, e.g., Figure 2 of Li et al. 2009, Figure 9 of Yu et al. 2010, and Figure 1 of Zhu et al. 2013).

Dust drifts relative to gas. In a 2D disk, ignoring azimuthal gradients and a small diffusion term, the continuity equation for the dust reads

[TABLE]

where

[TABLE]

is the radial drift velocity of dust, $p$ is gas pressure, and $t_{\rm s}={\rm St}/\Omega_{\rm k}$ is the dust stopping time. Dust will pile up at locations where the second term in (9) is most negative. This term is plotted at the bottom of Figure 1. In total, 5 local minima can be identified, each corresponding to a dust concentration. The big-dust-to-gas mass ratio, initially equal to 0.5% everywhere, increases at the locations of dust rings by factors of up to 4. The horseshoe dust ring HR (seen before in, e.g., Figure 4 of Zhu et al. 2014) has no obvious counterpart in $\Sigma_{\rm gas}$ and is caused by a different mechanism — 1:1 gravitational resonant interaction with the planet.

The dust rings are not necessarily located at the positions of low-order mean motion resonances relative to each other, although they may appear to be close to resonances by chance. The locations of the rings can change with planet mass, background disk profile, or planet migration, as will be seen later.

The behavior of dust in a disk with a sub-Jupiter mass planet has been investigated by many others (e.g., Paardekooper & Mellema, 2004, 2006; Lyra et al., 2009; Zhu et al., 2014; Rosotti et al., 2016) using two-fluid simulations. Our results generally agree with these works in the low $M_{\rm p}$ and low $\alpha$ limit (see, e.g., the multiple gaps and rings in the 0.05$M_{\rm J}$ model in Figure 8 of Paardekooper & Mellema 2006, and the 8$M_{\oplus}$ model in Figure 4 of Zhu et al. 2014). Note that our 10$M_{\oplus}$ planet is not massive enough to trap dust effectively at its 2:1 resonance at $r=1.6r_{\rm p}$, by contrast to the higher mass planet considered by Paardekooper & Mellema (2004, 2006).

Our planet excites density waves both inside and outside its orbit. The pressure structures in the density waves are not able to concentrate dust particles via dust-gas coupling because the waves corotate with the planet and not with the local disk. Dust particles only spend a small fraction of time inside the pressure overdensities associated with waves. By contrast, density waves excited by gravitational instability do rotate at approximately the local Keplerian velocity and are able to concentrate particles of certain sizes (Rice et al., 2004).

3.1.2 Synthetic Observations and Detectability

Figure 2 shows simulated $H$-band polarized intensity and ALMA Band 7 images. At $H$-band, the depletion of the gas (and thus of the entrained small dust) in the two gaps results in two shallow gaps in the full resolution image. After PSF convolution, however, the two gaps merge to form an apparent single gap. Depending on the noise level and the absolute brightness of the disk, such a shallow gap may or may not be detectable in a real system. In the full resolution mm dust continuum image, four dust rings can be distinguished: IR2, IR, OR, and a merged MR+HR. The simulated ALMA observation recovers nearly 100% of the total flux, and generally preserves the main features of the $\Sigma_{\rm dust}$ radial profile. In particular, in both the image and the radial profile, the signature of a super-Earth in a low viscosity disk — a pair of gaps (a “double gap”) separating three rings — is distinctive.

Figure 3 illustrates the detectability of gaps and rings under various observational conditions using ALMA. Assuming the disk is at Taurus (dec. $\sim 25^{\circ}$) and using array configuration alma.out24.cfg (angular resolution $0\arcsec.03\times 0\arcsec.02$), the distinctive double-gap feature can be clearly resolved in 2 hours of integration time (panel (b); RMS noise level 18 $\mu$Jy beam*-1*; S/N$\sim$22 at the OR), and is marginally detectable with 1 hour of integration time (panel (a)). With 4 hours of integration (panel (c)) all 4 rings can be detected with S/N$>$20 (although the separation of IR2 and IR is marginal). If the object is in the southern sky (dec. $-23^{\circ}$, panels (d)–(e)), the double-gap feature is more easily resolved because of better uv coverage. For array configuration alma.out21.cfg which yields a coarser angular resolution of $0\arcsec.05\times 0\arcsec.04$, the double-gap structure becomes difficult to discern (panel (f)).

Figure 4 shows the morphology of the rings and gaps in ALMA observations at various disk inclinations. In all cases shown, the disk has a position angle of $45^{\circ}$, in between the directions of the major and minor axes of the beam. At $45^{\circ}$ disk inclination, the structures in the disk can still be detected, particularly the double gap; at $60^{\circ}$ inclination, structures along the minor axis of the disk threaten to be lost, while rings and gaps can still be seen along the major axis.

Our hydro simulations are 2D and are thus not suitable for simulating the appearance of density waves in scattered light, as these waves have substantive vertical motions (Zhu et al., 2015). However, as judged from the 3D hydro and radiative transfer simulations of Dong & Fung (2017a), the density waves excited by low mass planets seem generally too weak to be detected by current NIR imaging facilities.

3.1.3 Time Dependence

Figure 5 shows the time evolution of the disk (azimuthally averaged radial profiles are shown in the last column of Figure 8). The observational signatures noted above are transient; in particular, the central ring (MR+HR) seen in the mm continuum gradually diffuses, merging the two gaps into one and eventually killing the signature of a double gap. Nevertheless, the double-gap feature can still last long enough to be detectable, about 2000 orbits (0.3 Myr for $r_{\rm p}=30$ AU) for the model shown.

3.1.4 Lengthening the Planet Growth Timescale

In this section, we explore the effect of lengthening the planet growth timescale. Figure 6 compares our fiducial model 10ME ($\tau_{\rm growth}$=$10$ orbits) at 1500 orbits, with two other models having $\tau_{\rm growth}=1000$ and 3000 orbits, sampled at 2000 and 3500 orbits, respectively. The three models show similar surface density profiles in both gas and dust. Growing the planet in a few thousand orbits does not materially affect the multiple rings and gaps ultimately produced by the planet.555Hammer et al. (2017) found that more slowly growing planets trigger weaker vortices at gap edges.

3.2. Dependence on Planet Mass

In this section, we explore how the morphologies of the rings and gaps vary with planet mass. Surface density maps and simulated images are shown in Figure 7, radial profiles are shown in Figure 8 (first column), and azimuthal profiles of selected models are shown in Figure 9. We find the following behavior when increasing $M_{\rm p}$ from $2M_{\oplus}$ to $60M_{\oplus}$:

Gas gaps deepen with increasing $M_{\rm p}$. As more gas is pushed to gap edges, $\Sigma_{\rm gas}$ in the IR and in the OR increases, leading to an increase in the gap-ring contrast in $H$-band images. In simulated ALMA images, the double-gap feature first appears at $M_{\rm p}=5M_{\oplus}$ and grows more prominent with increasing $M_{\rm p}$. 2. 2.

Each dust ring varies in its own way with $M_{\rm p}$. For the MR, $\Sigma_{\rm gas}$ and by extension $\Sigma_{\rm dust}$ decrease with increasing $M_{\rm p}$. By comparison, $\Sigma_{\rm dust}$ in the HR rises. As a result, as $M_{\rm p}$ increases in the simulated ALMA images, the location of the (practically merged) MR+HR shifts from slightly inward of the planet’s orbit to being more nearly coincident with the planet’s orbit. The IR, IR2, and OR all become more prominent in $\Sigma_{\rm gas}$ and $\Sigma_{\rm dust}$ as $M_{\rm p}$ increases. In simulated ALMA images, IR and IR2 are distinguishable only for $M_{\rm p}\gtrsim$ Neptune mass. 3. 3.

The OR is pushed outward. 4. 4.

As $M_{\rm p}$ increases, the HR becomes more azimuthally asymmetric. Figure 9 shows azimuthal profiles of $\Sigma_{\rm gas}$, $\Sigma_{\rm dust}$, and ALMA flux density in the coorbital region $r=r_{\rm p}-2h$ to $r=r_{\rm p}+2h$. In $\Sigma_{\rm gas}$, a 60 $M_{\oplus}$ planet starts to evacuate material in its vicinity to produce an azimuthal dip (compare with Figure 9 of Yu et al. 2010). A similar structure can also be found in the dust distribution and starts at a smaller $M_{\rm p}\sim$ Neptune mass, widening in angle with increasing $M_{\rm p}$. In addition, dust starts to accumulate around the triangular Lagrange points L4 and L5 $\pm 60^{\circ}$ away from the planet, as seen in previous works (see, e.g., the 40$M_{\oplus}$ planet case in Zhu et al., 2014, Figure 5). The simulated ALMA image with $M_{\rm p}=60M_{\oplus}$ shows two Trojan arcs (bottom right panel in Figure 7). Eventually, as $M_{\rm p}$ approaches a Jupiter mass (Lyra et al., 2009), the dust HR retreats into two azimuthally confined spots around L4 and L5, and the double-gap feature vanishes.

In general, dust gaps are more prominent than gas gaps, consistent with Paardekooper & Mellema (2006). We note that the minimum planet mass for the double gap to appear in observations ($5M_{\oplus}$ in our models) depends on disk parameters; for example, the minimum $M_{\rm p}$ decreases with decreasing scale height (data not shown).

3.3. Dependence on Viscosity

Viscous diffusion tends to smooth out density perturbations in the gas induced by planets. In Figures 8 (second column) and 10, we compare three simulations with a Neptune-mass planet and different values for $\alpha$: 5$\times$10*-4*, 5$\times$10*-5* (our fiducial model of Section 3.1), and 5$\times$10*-6*.

In the $\Sigma_{\rm gas}$ snapshots and the NIR images, as viscosity decreases the gaps deepen, as expected (e.g., Fung et al., 2014; Kanagawa et al., 2015; Duffell, 2015). A second concomitant effect is that the surface densities in the IR and OR, located at the gap edges, increase with decreasing viscosity: when $\alpha$ drops, more gas is pushed by Lindblad torques into the gap edges, and low viscosity makes it harder to diffuse away this excess. A third effect is that at the highest viscosity ($\alpha=5\times 10^{-4}$), the HR+MR in the gas disappears and the two gaps merge. High viscosity enables gas to more easily diffuse away from the co-orbital region, which is a potential maximum. In addition, for planets with $M_{\rm p}\lesssim M_{\rm thermal}$, viscous diffusion can prevent waves from shocking and damping and can thereby suppress gap formation (Goodman & Rafikov, 2001; Dong et al., 2011a). All three effects have been seen in previous hydro simulations (e.g., Li et al., 2009; Yu et al., 2010).

In the $\Sigma_{\rm dust}$ snapshots and simulated ALMA observations, the two low viscosity runs produce qualitatively similar ring and gap structures. By contrast, in the $\alpha=5\times 10^{-4}$ run, the two gaps merge and only two broad overdensities remain at the gap edges. A wide single gap may not be distinguishable from the gap opened by a more massive planet in a disk with higher viscosity. Thus, low viscosity ($\alpha<10^{-3}$) is a necessary condition for super-Earths to induce their signature double gaps. The two low viscosity models differ, however, in that for $\alpha=5\times 10^{-6}$, a vortex is created, as indicated in Figure 11 using a more dramatic color stretch. The development of vortices at edges of planet-opened gaps requires low viscosity (e.g., Zhu & Stone, 2014), and is probably driven by the RWI (e.g., Li et al., 2000, 2001). As with dust trapping by radial pressure structures, pressure peaks in azimuth associated with vortices also trap dust particles of certain sizes (e.g., Lyra & Lin, 2013). The vortex in 17ME-LowVis is quite visible in the simulated ALMA image (Figure 11, right panel).

3.4. Dependence on the Global Background Disk Profile

The dust distribution, a direct observable from ALMA, is determined by gas drag forces, which depend on local gas pressure gradients. These local gradients in turn depend on the global background gradient. Figures 8 (third column) and 12 compare models 10ME-Beta15 ($\Sigma_{0}\propto 1/r^{1.5}$), 10ME ($\Sigma_{0}\propto 1/r$), and 10ME-Beta0 ($\Sigma_{0}\propto$ constant). The three models share similar radial profiles for $\Sigma_{\rm gas}/\Sigma_{\rm gas,0}$, but the locations and contrasts of the dust rings depend on $\beta$, the power-law slope of the background surface density. As $\beta$ decreases, the IR, HR, MR, and OR in $\Sigma_{\rm dust}/\Sigma_{\rm dust,0}$ expand outward. Overall, in the simulated ALMA observation, the 4 rings in model 10ME-Beta0 appear to be more distinctive and well separated than in the other two models.

3.5. Effects of Planet Migration

Finally, we experiment with the planet’s radial migration. We manually move the planet in models 10ME-Beta0 and 17ME across the disk at a constant user-controlled migration rate to explore how planet-induced structures depend on migration.

Figures 13 and 14 show the results. In models 10ME-Beta0-MigrateIn and 17ME-MigrateIn, the planet migrates inward at a rate $\dot{r}=-2\times 10^{-3}$AU orbit*-1* (or about $-10^{-5}$AU year*-1*). In models 10ME-Beta0-MigrateOut and 17ME-MigrateOut, the planet migrates at the same rate but outward. In all four models, the planet crosses $r_{\rm p}=30$ AU at the end of the simulation after 1500 orbits. We note that it takes 1500 orbits for the planet to migrate 3 AU (= 2 scale heights), about the distance between the nearest gap and the planet’s orbit. We find that when the planet migrates (either inward or outward), those features radially behind the planet lag further behind, while features radially ahead of the planet are squeezed closer. The one exception is the dust HR, which simply follows the planet. The radial shifts of the IR2, IR, MR, and OR could, in principle, be used to infer if a planet is migrating. For example, if the HR+MR ring in ALMA images is dominated by the HR, the gap ahead of a migrating planet narrows, while the gap behind the migrating planet widens, as shown in the left panel of Figure 14. However, reading the signature of migration would require knowledge of the background disk profile to determine the “intrinsic” locations of these features in the absence of migration.

3.6. Future Improvements

Our investigations of observational signatures of super-Earths can be extended for greater realism.

Simulations with a more realistic dust size distribution are a natural extension (e.g. Ruge et al., 2016). In this work we chose 200 $\micron$ as a representative size for the “big” dust particles. But dust of other sizes can have a non-trivial opacity at $\sim$mm wavelengths. 2. 2.

Relatedly, grain growth and evolution can be included. Once dust is concentrated into rings, coagulation can be facilitated (e.g., Youdin & Shu, 2002; Johansen et al., 2007; Pinilla et al., 2016), and instabilities such as the streaming instability (Youdin & Goodman, 2005; Chiang & Youdin, 2010) can be triggered. Our simulations focus on the first stage, the drift of $\sim$mm-sized dust into gas pressure perturbations. A natural extension would be to follow the subsequent evolution of dust in the resulting concentrations. 3. 3.

Our simulations assume a non-zero $\alpha$ viscosity. Truly inviscid disks can behave differently (Fung & Chiang, 2017). 4. 4.

Our simulations are 2D. For planets with $M_{\rm p}>M_{\rm thermal}$, Fung & Chiang (2016) have shown that gap opening proceeds in 3D much the same way that it does in 2D. It will be interesting to explore 3D effects of sub-thermal-mass planets in low viscosity disks.

4. Summary and Discussion

In this work, we combined two-fluid hydrodynamics and radiative transfer simulations to calculate the distributions of gas and of $\sim$mm-sized dust in a low viscosity disk perturbed by a super-Earth, and created synthetic observations of such systems in both near-infrared (NIR) scattered light and ALMA mm dust continuum emission. Our main findings are:

A $\sim$10 $M_{\oplus}$ planet in a low viscosity disk initially opens two narrow annular gaps, each displaced radially by about a vertical pressure scale height to either side of the planet’s orbit. Perturbations in the gas pressure cause dust particles to drift differentially. Millimeter-sized dust particles concentrate into three distinctive rings: two rings at the far edges of the two gaps, and one ring sandwiched between the two gaps (Figure 1). The 1:1 resonance may produce a fourth dust ring exactly coincident with the planet’s orbit. A fifth dust ring interior to the double-gap region is also possible, depending on parameters. These dust rings are not necessarily located at the positions of low-order mean motion resonances relative to each other, although they may appear to be so by chance. For our assumed disk parameters, a minimum planet mass of $\sim$5$M_{\oplus}$ is needed to produce detectable rings on $\sim$1000 orbits timescale. Planets more massive than Neptune ($\gtrsim 17M_{\oplus}$) generate azimuthal dust concentrations at the triangular Lagrange points L4 and L5. 2. 2.

The dust rings and gaps produced by a super-Earth at 30 AU in a typical disk with \sim$$10 $M_{\oplus}$ of mm-sized dust at 140 pc distance can be readily detected by ALMA in mm dust continuum observations with an angular resolution of $0\arcsec.03$ and 2 hours of integration. In NIR scattered light, the gaps and rings are less prominent and may not be detectable with current NIR capabilities. 3. 3.

To produce the signature double-(or-multiple)-gap feature in mm continuum observations, the disk’s $\alpha$ viscosity parameter has to be lower than $\sim$10*-4*. At higher viscosity, only one gap is opened at the planet’s orbit. If $\alpha$ is lower than $\sim$10*-5*, the Rossby wave instability at gap edges can be triggered, forming dust-trapping vortices. 4. 4.

The double-gap signature of a single super-Earth manifests in disks with a variety of background profiles and tolerates modest planet migration (10 AU / Myr at 30 AU). Migration changes the locations of the rings in potentially observable ways. Although the double gap can be discerned for thousands of orbits, it weakens as material diffuses away from the co-orbital region.

At tens of AUs, super-Earths can produce a characteristic double-gap feature detectable in high resolution mm continuum observations by ALMA. A low disk viscosity of $\alpha\lesssim 10^{-4}$ is critical for observability. Super-Earths, unlike their Jovian siblings, induce only weak perturbations in the gas disk, which must be nearly inviscid for such perturbations to manifest observationally (in particular, to facilitate the development of the Goodman & Rafikov 2001 wave nonlinearity), and for aerodynamic drag on dust to amplify these features. The low viscosity at the disk midplane that is advocated here finds support from both observations and theory. The leading explanation for lopsided disks observed by ALMA (e.g., van der Marel et al., 2013; Casassus et al., 2013; Pérez et al., 2014)—vortices formed at planet-opened gap edges—requires $\alpha\lesssim 10^{-4}$ (e.g., Zhu & Stone, 2014). Similarly, if most cavities in transitional disks are opened by multiple planets, $\alpha$ needs to be lower than $10^{-3}$ in order to mobilize enough low-mass planets to account for the prevalence of transitional disks (Dong & Dawson, 2016). See also Fung & Chiang (2017) for how multiple super-Earths in inviscid disks can both survive orbital migration and drive disk accretion. Magnetohydrodynamics (MHD) simulations have shown than non-ideal MHD effects, in particular ambipolar diffusion, can significantly suppress the magnetorotational instability in disks at tens of AUs, resulting in practically zero turbulence at disk midplanes (Bai & Stone, 2011; Perez-Becker & Chiang, 2011b, a; Bai, 2013; Turner et al., 2014; Bai, 2015). Not even in disk surface layers has turbulence been detected (Flaherty et al. 2015; Flaherty et al. 2017, in prep.).

A double gap observed in nature could, naively, be explained with two planets in a standard viscous disk, each opening their own gap. But this interpretation is problematic. The two gaps in the double-gap feature that we found are typically $\sim$2–4 scale heights apart, or 0.1–0.2 $r_{\rm p}$ for $h/r=0.05$. In a disk with $\alpha>10^{-4}$, the width of a single planet-opened gap is greater than $0.2r_{\rm p}$ (e.g., Kanagawa et al. 2016, Figure 3; Dong & Fung 2017b, Figure 7). Two closely spaced gaps opened by two planets would merge into one in a standard viscous disk. The two-planet configuration may also be prone to orbital instability, particularly if the planets are not in mean-motion resonance (Tamayo et al., 2015).

We close with a few remarks on real disks. To date, two systems, HL Tau (ALMA Partnership et al., 2015) and TW Hya (Andrews et al., 2016; Tsukagoshi et al., 2016), have been observed by ALMA with sub-$0\arcsec.05$ resolution. Intriguingly, in each case a candidate double gap can be identified: D5 and D6 at 64 and 74 AU in HL Tau, and the two gaps at 37 and 43 AU in TW Hya. We plan to discuss these two systems and explore the possibility that their ring/gap features are produced by super-Earths in a forthcoming paper. In addition, a Large Program (PI: S. Andrews) was approved in ALMA Cycle 4 to resolve 20 nearby disks at high angular resolution. We look forward to seeing whether this sample also features multiple rings, and ultimately to inferring whether super-Earths commonly form at large distances, as they appear to have done at small ones.

Acknowledgments

We are grateful to an anonymous referee for constructive suggestions that improved our paper. R.D. thanks Bekki Dawson, Jeffrey Fung, Eve Lee, Paola Pinilla, Roman Rafikov, and Zhaohuan Zhu for insightful discussions; Fumi Egusa at the ALMA help desk on using CASA; and Sean Andrews and David Wilner for motivating this work. Part of the numerical calculations were performed on the SAVIO cluster provided by the Berkeley Research Computing program, supported by the UC Berkeley Vice Chancellor for Research and the Berkeley Center for Integrative Planetary Science. H. Li and S. Li gratefully acknowledge the support by LANL’s LDRD program and a CSES project. E.C. acknowledges support from the NSF and NASA.