Multiple Disk Gaps and Rings Generated by a Single Super-Earth. II. Spacings, Depths, and Number of Gaps, with Application to Real Systems

The Atacama Large Millimeter/submillimeter Array (ALMA) has found multiple dust gaps and rings in a number of protoplanetary disks in continuum emission at millimeter wavelengths. The origin of such structures is under debate. Recently, we documented how one super-Earth planet can open multiple (up to five) dust gaps in a disk with low viscosity (α ≲ 10−4). In this paper, we examine how the positions, depths, and total number of gaps opened by one planet depend on input parameters, and apply our results to real systems. Gap locations (equivalently, spacings) are the easiest metric to use when making comparisons between theory and observations, as positions can be robustly measured. We fit the locations of gaps empirically as functions of planet mass and disk aspect ratio. We find that the locations of the double gaps in HL Tau and TW Hya, and of all three gaps in HD 163296, are consistent with being opened by a sub-Saturn mass planet. This scenario predicts the locations of other gaps in HL Tau and TW Hya, some of which appear consistent with current observations. We also show how the Rossby wave instability may develop at the edges of several gaps and result in multiple dusty vortices, all caused by one planet. A planet as low in mass as Mars may produce multiple dust gaps in the terrestrial planet-forming region.

At first glance, multiple gaps seem to imply multiple planets, with one planet embedded within each gap. Recently, however, it has been appreciated that a single planet can generate multiple gaps. Using gas+dust simulations,  showed that one super-Earth planet (with a mass between Earth and Neptune) can produce up to five dust gaps in a nearly inviscid disk, and that these gaps would be detectable by ALMA (see also Muto et al. 2010;Duffell & MacFadyen 2012;Zhu et al. 2013;Bae et al. 2017;Chen & Lin 2018;Ricci et al. 2018). In particular, a pair of closely spaced narrow gaps (a "double gap") sandwiches the planet's orbit. This structure is created by the launching, shocking, and dissipation of two primary density waves, as shown in Goodman & Rafikov (2001, see also Rafikov 2002a, 2002b. Additional gaps interior to the planet's orbit may also be present, possibly opened by the dissipation of secondary and tertiary density waves excited by the planet (Bae & Zhu 2018a, 2018b. Motivated by the boom of multi-gap structures discovered by ALMA, and the mounting evidence of extremely low levels of turbulence at the disk midplane (from direct gas line observations, e.g., Flaherty et al. 2015Flaherty et al. , 2017Teague et al. 2018, and from evidence for dust settling, e.g., Pinte et al. 2016;Stephens et al. 2017), we investigate here what properties of the planet and disk can be inferred from reallife gap observations. We study systematically how gap properties depend on planet and disk parameters, and propose generic guidelines connecting ALMA observations with diskplanet interaction models. We introduce our models in Section 2, and present our main parameter survey of gap properties in Section 3. These results are applied to three actual disks in Section 4, discussed in Section 5, and summarized in Section 6.

Simulations
We use two-fluid (gas + dust) 2D hydrodynamics simulations in cylindrical coordinates (radius r and azimuth f) to follow the evolution of gas and dust in a disk containing one planet. The code used is LA-COMPASS (Li et al. 2005(Li et al. , 2009). The numerical setup, largely adopted from , is summarized here (more details can be found in Fu et al. 2014 andMiranda et al. 2017).
We adopt a locally isothermal equation of state: pressure = S P c gas s 2 , where S gas and c s are the surface density and sound speed of the gas. Dust is modeled as a pressureless second fluid. The following Navier-Stokes equations are solved: t gas gas gas gas gas where D dust is the dust diffusivity defined as = D dust n t + ( ) 1 stop 2 (Takeuchi & Lin 2002), F G is the gravitational potential (including both the gravity of the planet and the self-gravity of the disk; a smoothing length of 0.7h is applied to the gravitational potential of the planet), is the drag force on dust by gas, t stop is the dimensionless stopping time (see below), and n f is the Shakura-Sunyaev viscous force with n a = W h 2 K , h the disk scale height, W K the Keplerian angular frequency, and a < 1 a dimensionless constant. Other symbols have their usual meanings. The last term in Equation (3) accounts for the dust back-reaction on gas.
The initial surface densities of gas and dust are (we use the subscript 0 to indicate initial conditions): is the radius of the planet's circular orbit, which is fixed in our simulations. We note that if a planet migrates in the radial direction fast enough, gap properties can be affected (Dong et al. 2017, Figure 14). However, low-mass planets in near-invicsid environments as studied here are not expected to migrate significantly (Li et al. 2009;Yu et al. 2010;Fung & Chiang 2017). We ramp up the planet mass with time, t, as µ ( ) ( ) M t t sin p over 10 orbits. In   Figure 6) we experimented with different planet growth timescales up to a few thousand orbits, and found no significant effects on the rings and gaps ultimately produced by the planet.
The disk aspect ratio h/r is assumed radially invariant. We adopt a constant instead of a flared h/r profile to suppress the dependences of gap properties on the h/r gradient, and to increase the timestep at the inner boundary. We note that fractional gap spacings, which we will use as our primary diagnostic of models and observations (see our Equation (8)), depend only weakly on radial gradients in h/r, Σ gas,0 , and the slope of the background S ( ) r gas profile. Dong et al. (2017, Figure 8) showed that as the background S ( ) r gas profile varies from S µ r gas 0 to S µr gas 1.5 , all gaps and rings shift inward, while their relative separations are unchanged at the level of a few percent.
The dimensionless stopping time of dust-the particle momentum stopping time normalized to the dynamical time-is, where r = 1.2 particle g cm −3 is the internal density of dust particles, and s=0.2 mm is the particle size. Such particles are probed by continuum emission at mm wavelengths (e.g., Kataoka et al. 2017). At r=0.2-r 1 p , t stop is on the order of 0.01. We note planet. The green plus and dot mark the locations of the star and the planet, respectively. The inner region is manually masked out. Bottom: azimuthally averaged radial profiles of surface density perturbations in dust (solid line) and gas (dashed line) with gaps and rings labeled: IRn (the nth inner ring from the planet), IGn (the nth inner gap), MR (middle ring), HR (the horseshoe ring), ORn (the nth outer ring), and OGn (the nth outer gap). Density fluctuations in gas are much smaller than in dust. See Section 2 for details. that for mm-to-cm sized particles often traced by VLA cm continuum observations (possibly having t~-0.1 1 stop ), gaps and rings similar to those in 0.2 mm-sized dust are generated. Qualitatively, as particle size increases, gaps widen and deepen with their locations largely unchanged; the amplitude of the middle ring (MR) co-orbital with the planet drops; particle concentrations at the triangular Lagrange points L 4 and L 5 become stronger; and an azimuthal gap centered on the planet's location emerges and widens (Lyra et al. 2009;Zhu et al. 2014;Rosotti et al. 2016;Ricci et al. 2018).
The simulation domain runs from 0 to 2π in f, and 0.1 to r 2.1 p in r. The initial gas disk mass is  M 0.015 . We set up a wave damping zone at = r r 0.1 0.2 p as described in de Val-Borro et al. (2006), and an inflow boundary condition at the outer boundary. The resolution is discussed in Section 2.1.
We summarize our models in p t h , we use even lower viscosities α=10 −6 -10 −5 to speed up gap opening (cf., Fung et al. 2014). Figure 1 shows the 2D gas and dust distributions and their azimuthally averaged 1D radial profiles for Model H003MP02 at 1600 orbits (0.26 Myr at 30 au). A = Å M M 0.2 1.8 th planet produces multiple gaps and rings in both gas and dust, while the surface density perturbations in dust are an order of magnitude larger than they are in gas. In each of the dust rings, a few Earth masses of dust are trapped. We refer to the gaps and rings inside r p as "inner" gaps and rings (IG and IR), and the ones outside r p as "outer" gaps and rings (OG and OR), and number them according to their distances from the planet.

Numerical Resolution in Hydro Simulations
Sufficient spatial resolution is needed to properly capture dust gaps, as resolution affects numerical diffusivity. In inviscid or nearly inviscid disks, the angular momentum carried by a wave launched by a planet should be conserved until the wave shocks and deposits its angular momentum locally, and opens a gap (see the weakly nonlinear theory of density waves; Goodman & Rafikov 2001). Insufficient resolution leads to excessive numerical diffusion; gaps either appear too close to the planet and are too shallow (Yu et al. 2010;Dong et al. 2011bDong et al. , 2011a, or fail to be seen at all (Bae et al. 2017). Figure 2 compares the radial profiles of S S dust dust,0 from five runs of Model H004MP02, labeled A, B, C, D, and E in order of increasing spatial resolution. Runs A, C, and E successively double the number of radial cells (n r ) used, while runs B, C, and D successively double the number of azimuthal cells (n f ). The three runs with  n 4608 r and  f n 3072 (C, D, and E) achieve visual convergence at  r r 0.4 p . Note that insufficient resolutions may also result in a reduced number of detectable density waves ( Figure 3, see also Bae et al. 2017). A higher resolution is needed in the radial direction because density waves have higher spatial frequencies in the radial direction, and gaps are radial structures. At n r ×n f =4608×3072 for Model H004MP02, the Hill sphere of the planet is resolved by 75 and 16 cells in the radial and azimuthal directions, while h is resolved by 96 and 20 cells, respectively.
For all models shown in Table 1, we adopt resolutions necessary to ensure visual convergence of results for all disk radii beyond at least r 0.4 p (2× the outer radius of the wave damping zone). For runs with h/r>0.05, the resolution requirements are less severe, and we achieve convergence for radii outside r 0.3 p .

Results
In this section, we study how the properties of gaps and rings vary with input parameters. Most simulations are run for ∼10 3 -10 4 orbits (∼0.1-1 Myr at 30 au), after which time dust gaps are depleted by a factor of ∼2-10 in surface density, The relative separation (spacing) between two radii (as drawn from r gap , r ring , or r p ) is defined as In real observations, r p is usually unknown, and the spacings of gaps and rings are the direct observables.

General Comments on Model-observation Comparisons
There are at least three observable gap properties: gap depth (contrast with adjacent rings), gap width, and gap location. In model-observation comparisons, gap depth and width are not easy metrics to work with. First, narrow gaps are often  Figure 2. Each run corresponds to one horizontal stripe, and adjacent runs are separated by a dotted line. Note that each map only covers 1/4 of the simulation domain in the azimuthal direction (1/12 for each run). The planet is located at (r,f)=(r , 0 p ) outside the box to the right. In the dust maps (bottom), annular gaps and rings are vertical stripes. If the adjacent two runs achieve convergence, their wave and gap patterns should smoothly join each other at the boundary. The green arrows in the upper panel mark a density wave present in runs C and E but not in A. Simulations with insufficient spatial resolutions produce a smaller number of density waves in gas, and significantly different dust distributions. See Section 2.1 for details. in Model H003MP02 at four epochs. Local maxima and minima on each profile are marked by dots. S dust at OG1 is depleted by 30%, 50%, 70%, and 90% at 700, 1100, 1600, and 2600 orbits, respectively. The locations of all gaps stay roughly the same, while IR1 and OR1 shift away from the planet, as indicated by arrows. See Section 3.2 for details. unresolved in observations, sabotaging measurements of true depth and width-in the ALMA Partnership et al. (2015) observations of HL Tau, the observed gap widths are comparable to the beam size (Akiyama et al. 2016 , Table 2). Second, spectral index analyses have shown that bright rings between gaps may be optically thick (e.g., ALMA Partnership et al. 2015;Liu et al. 2017;Huang et al. 2018), which implies that dust surface density does not translate into a unique surface brightness. Third, gap spacings are less sensitive to the background S gas profile, often hard to determine, than gap depths , Figure 3). Finally, gap depth and width evolve with time (Section 3.2), and how long a (hypothetical) planetary system has perturbed a real disk is usually unknown.
In contrast, it is simpler and more robust to use gap locations (spacings) when comparing models with observations. Gap locations evolve less with time (Section 3.2), and can be determined relatively accurately in observations even when gaps are not well resolved and adjacent rings are optically thick. The total number of gaps may also provide a diagnostic. Bae & Zhu (2018b, Figure 3) showed that more density waves can be discerned in colder disks, or in disks perturbed by lower mass planets. Figure 4 shows the time evolution of a representative Model H003MP02. The radial profiles of S S dust dust,0 at 700, 1100, 1600, 2600 orbits are plotted. At these times, S dust at OG2 (outer gap 2) is depleted by 30%, 50%, 70%, and 90%, respectively. Gap depths are still evolving at 2600 orbits (0.43 Myr at 30 au). Rings become more optically thick and gaps widen with time, making gap contrast and width problematic metrics in model-observation comparisons.

The Time Evolution of r gap and r ring
While gaps deepen noticeably with time, their locations do not change as much. From 700 to 2600 orbits, IG1 (inner gap 1) moves away from the planet by 15%; all other gap locations drift by less than 2%. In comparison, the two most prominent dust rings, OR1 (outer ring 1) and IR1 (inner ring 1), move away from the planet by 17% and 25%, respectively. The fractional separation (a direct observable) D OG1,IG1 increases by 6%, while D OR1,IR1 increases by 21% over the course of the simulation.
We conclude that gap locations (spacings) serve as the most robust metrics for model-observation comparisons. Tables 2 and 3 in the Appendix list the locations and values of S S dust dust,0 , respectively, of rings and gaps in selected models at a time when OG1 is ∼50% depleted (for the exact times, see Table 1, second to last column). Figure 5 plots the radial profiles of S S dust dust,0 in models with varying h/r (top) and M p (bottom). The corresponding 2D dust maps can be found in the Appendix Figure 13. All models are shown when surface densities inside gaps are ∼50% of their original values.

The Dependence of Gap Spacing on h/r and M p
As h/r increases, or M p decreases, gaps move away from the planet and become more widely separated. Here we explain these trends in the framework of the weakly nonlinear density wave theory (Goodman & Rafikov 2001;Rafikov 2002a), and provide empirical fitting functions for r OG1 , r IG1 , and r IG2 .
Under the local approximation (constant S gas and h/r in the planet's vicinity), Goodman & Rafikov (2001) found that the primary density waves launched to either side of the planet's orbit shock after traveling a radial "shocking length" equal to where γ is the adiabatic index (γ=1 for isothermal gas). The shocks mark the onset of wave dissipation and gap opening in gas. For planets withM 0.1 p -M 0.5 th ,~-l h 1 2 sh . We expect the shock locations  r l p s h to mark the locations of dust gaps IG1 and OG1, though the identification is not perfect-the shock locations correspond to the edges of gas gaps, while dust gaps have finite radial widths, and have different profiles from gas gaps. Another issue is that Equation (9)  2 5 -fit the behaviors of r OG1 −r p and r p −r IG1 well, as shown in Figure 6. Consequently, r OG1 −r IG1 , a.k.a.the "double gap" separation, obeys the same scalings. Empirically, we find that Compared with those of the double gap, the dependences of r p −r IG2 on both h/r and M p are weaker. Empirically fitting the locations of gaps further inward of IG2 is possible, but only a small number of models robustly reveal such gaps and so we have not performed this exercise.

The Number and Depths of Gaps
In contrast to gap spacings, the total number of gaps opened by a planet evolves more strongly with time. As they deepen, more gaps become visible (e.g., Figure 4). How fast a gap deepens (d S dust /dt) depends on a variety of factors. The following are a few general observations. 1. Different gaps deplete at different rates. For example, while in Model H003MP02 OG1 and IG1 deepen at roughly the same rates, gaps farther from the planet's orbit deplete more slowly (Figure 4). 2. At fixed M M p th , gaps open faster in disks with higher h/r. Figure 5(a) shows that in such disks it takes less time to deplete OG1, and IG1 and IG2 deepen even faster. Figure 7 (a) shows S S dust dust,0 at OG1, IG1, and IG2 for these models. 3. More massive planets open gaps faster. Figure 5(b) shows that a planet 10× more massive depletes OG1 ∼100× faster. This effect is weaker for gaps farther from the planet. Figure 7(b) shows S S dust dust,0 at OG1, IG1, and IG2 for the same models in Figure 5(b).
Gaps at small radii are difficult to count in both models and observations because of limited resolution. We list the number of dust gaps between r 0.4 p and OG1 (inclusive; OG2 is never prominent) in Table 1 (last column). A gap is counted as long as its S d dt dust is consistently negative (i.e., no minimum threshold in S S dust dust,0 is imposed). Not all gaps are deep enough to be detectable in realistic observations.
In general, a planet opens more gaps in disks with lower h/r-gaps in such a disk are narrower, allowing more gaps to fit within a given radius range (e.g.,   r r r 0.4 p OG1 ). The trend is less clear with planet mass. There are two competing factors. Bae & Zhu (2018b) discerned more density waves with decreasing M p , potentially resulting in more gaps. On the other hand, gaps opened by a less massive planet are more widely separated, decreasing the maximum number of gaps within a given radius interval. It is unclear which effect dominates. Tentatively, models H003MP01-H003MP08 each have five gaps at  r r 0.4 p , suggesting that the number of gaps may be insensitive to M p (note that Model H003MP005 has only three gaps).

Comparison with Real Disks
In this section, we make a quick comparison between our models and observed disks, focusing on those having three or more narrow gaps in mm continuum emission (HL Tau, TW Hya, and HD 163296). These comparisons do not involve detailed fits. Rather, the models are taken directly from our parameter survey grid, and are not specifically tailored to the observations. We focus on gap locations and not on gap depths for reasons discussed in Section 3.1. We begin by picking h/r consistent with disk midplane temperatures, T midplane , as derived from the literature. We then identify a pair of gaps as the "double gap." For HL Tau and TW Hya, this is the most compact pair, and for HD 163296, it is the pair farthest from the star. The midline of the double gap locates the planet position r p . We then estimate M p using Equation (10), selecting the one model from dust dust,0 in 5 models with h/r=0.03 and M p successively doubled. All models are shown at the time when OG1 reaches 50% depletion. OG1, IG1, and IG2 are marked by red diamonds, green squares, and blue circles, respectively. As h/r increases, or M p decreases, gaps move away from the planet and become more widely separated. In addition, gaps further from the planet become deeper relative to OG1 (dashed arrows). See Section 3.3 for details.
our grid that appears to best reproduce the double gap locations. A non-trivial test of the model is to see whether it can reproduce the locations of observed gaps interior to the double gap.
The observed surface brightness radial profiles ( ) I r mm are shown in Figure 8, with modeled planet and gap locations indicated. We mark all modeled gaps at  r r 0.4 p regardless of their depths, but do not show the model surface brightness radial profile, again because our intent is not to fit gap depths. We attempt to remove the dependence of I mm on dust temperature by scaling I mm by r . Figure 9 shows side-by-side comparisons of observations with synthetic model ALMA continuum images, produced using radiative transfer simulations using the Whitney et al. (2013) code. The procedures are described in Dong et al. (2015, see also Dong et al. 2017. Briefly, 2D hydro disk structures are inflated into 3D, assuming hydrostatic equilibrium with each model's h/r. Starlight is reprocessed to longer wavelengths at the disk surface by a population of "small" (sub-μm-sized) interstellar medium dust (Kim et al. 1994). These small grains are assumed to be well-coupled to gas with a volumetric mass density equal Figure 6. Left: gap spacings as a function of h/r for the 7 models in Figure 5(a). The locations of OG1, IG1, and IG2 are plotted (|r gap − | r ; p symbols are the same as in Figure 5), in addition to the distance between OG1 and IG1, r r OG1 IG1 (black stars). The dotted and dashed lines are the fitting functions (10) and (11), respectively. Right: similar to (a), but for gap spacings as a function of M p for the 5 models in Figure 5  ). Right: the dependence of gap depth on M p , showing the 5 models in Figures 5(b) and 6(b). As h/r increases, or M p decreases, the gaps inside r p deepen relative to OG1 (more so for the inner ones). See Section 3.4 for details.
to 10 −2 that of the gas. A population of "big" dust particles (representing the second fluid in hydro simulations) is assumed to have settled to the midplane with a Gaussian vertical density profile and a scale height equal to h dust =0.1h, in line with h dust estimated in real disks (e.g., HL Tau, ∼1% at 100 au; Pinte et al. 2016). The optical properties of the small dust can be found in Dong et al. (2012, Figure 2). The big dust is assumed to have an opacity κ=3.5 cm 2 g −1 at 870 μm. For the central source we assume a  M 1 pre-main-sequence star with R å =2.4 R e and T å =4400 K. Radiative transfer simulations produce synthetic images of continuum emission at 870 μm, which are processed using the simobserve and simanalyze tools in Common Astronomy Software Applications (CASA) to simulate ALMA images with 4 hr integration times and default settings. . The corresponding Model H005MP08 has four gaps at  r r 0.4 p : r IG3 = 32 au, r IG2 =48 au, r IG1 =64 au, and r OG1 =78 au for = r 71 au p . Figure 9 (top) shows an image comparison. These four gaps roughly coincide with D2, D4, D5, and D6, respectively. D1, D3, and D7 in HL Tau find no ready correspondence with our model. Bae et al. (2017)  according to Bae et al. 2017), mainly due to our lower h/r. TW Hya. Continuum observations at 0.87 mm reveal three gaps at 25, 41, and 49 au (G1, G2, and G3, Andrews et al. 2016; these distances have been rescaled using the Gaia Collaboration et al. 2018 distance to TW Hya of 60 pc). Gap G3 is better revealed in the spectral index map than in the  (Loomis et al. 2017). These profiles are azimuthally averaged (except in HL Tau, where the profile is measured along the disk major axis), and scaled by r to suppress (partially) the dependence of I mm on dust temperature. The vertical solid line in each panel marks the location of the planet r p in our models, and the vertical dashed lines mark the locations of all modeled gaps at  r r 0.4 p (regardless of depth). Models are selected from our parameter grid such that the locations of gaps IG1 and OG1 (comprising the "double gap") are approximately reproduced in each case. In HD 163296, we are able to match a third gap IG2. In TW Hya, none of the predicted gaps (apart from IG1 and OG1, which fit by construction) appear to match. In HL Tau, IG3 and IG2 may correspond to observed features, while the observed gaps at 14 and 42 au (D1 and D3) are not reproduced (D1 lies too close to the inner boundary of our simulation to be accurately modeled; see Figure 2). In Elias 24, AS 209, and AA Tau, we introduce simple models to explain the two gaps revealed by existing observations, assuming they compose the double gap and need to be updated if future observations reveal additional substructures. The two gaps in the V1094 Sco disk (not shown) are too far apart to be explained by models with  h r 0.08. In closing this section, we note that gas gaps in our models are shallow (dS S~-1 10% gas gas,0 ), and so are not expected to be easily detectable in near-infrared (NIR) imaging, which traces (sub-)μm-sized dust typically well-coupled to gas (Dong & Fung 2017). This expectation generally accords with observations. In HD 163296, a broad gap is present at ∼0  - . 2 0. 5 in scattered light (Garufi et al. 2014;Monnier et al. 2017), but no counterparts to the three narrow ALMA dust gaps at   r 0. 5 are evident in the NIR (however note the detection of line intensity depressions in CO isotopologues at locations similar to the dust gaps; Isella et al. 2016). In AA Tau, no scattered light gaps are seen at the locations of mm-dust gaps (Cox et al. 2013). In TW Hya, G2 and G3 are not detected in scattered light; however, G1 is revealed at a slightly smaller radius (Akiyama et al. 2015;Rapson et al. 2015;van Boekel et al. 2017;Debes et al. 2017), indicating a large depletion in gas. High fidelity NIR disk imaging is not available for HL Tau (but note the tentative detection of gas gaps in HCO + , Yen et al. 2016), AS 209, and Elias 24.

The Observed Double Gap
To explain the double gap in the HL Tau and TW Hya disks, two mechanisms have been proposed: dust sintering and evolution at snowlines (Zhang et al. 2015;Okuzumi et al. 2016), and planet-disk interactions along the lines of our models. The former needs the snowlines of two relevant volatiles to be rather close (separation/radius ratio ∼17%), while the latter posits a sub-Saturn planet at dozens of au. While accurate assessments of snowline locations are needed to test the former hypothesis, the slow variation of disk temperature with radius may be difficult to make work.
Recently, an eccentric double gap withẽ 0.1 was tentatively discovered in the MWC 758 disk (Dong et al. 2018). A non-zero eccentricity may be compatible with a planet on an eccentric orbit, but is not natural to the snowline mechanism, which has no explicit azimuthal dependence. Evidence for little-to-no turbulence at the disk midplane at tens of aus is accumulating both theoretically (e.g., Perez-Becker & Chiang 2011;Bai & Stone 2013;Bai 2015;Schlaufman 2018) and observationally (e.g., Pinte et al. 2016;Flaherty et al. 2017, see also the discussion in Fung & Chiang 2017, final section). If α10 −4 is common, narrow dust gaps should mostly appear in pairs or multiples.

Low Viscosity and the Rossby Wave Instability
In low viscosity environments, gas gaps sufficiently deep are prone to the Rossby wave instability (RWI; e.g., Li et al. 2001Li et al. , 2005 at their edges. Gaps in this paper are generally too shallow to develop the RWI. Experiments (not shown) suggest that S gas needs to be depleted by 50% to trigger the RWI when α=5×10 −5 . The RWI may develop at the edges of more than one gap, forming multiple dust-trapping vortices. Figure 11 shows an example Model RWI, in which the RWI is triggered at both the inner edge of IG1 and the outer edge of OG1, forming two vortices. Dust are also collected at the triangular Lagrange points L 4 and L 5 . In total, four dust clumps are produced.

Implications for Planet Formation
An emerging trend from recent ALMA high resolution disk surveys is that multi-gap structures at 10-100 au are common in disks around a variety of host stars (S. Andrews & F. Long 2018, private communication). This demands a robust gap opening mechanism insensitive to disk and host star properties. We argue that such structures can be produced by one or several planets ranging in mass from 0.1 to 10s of Earth masses, located at 1-100 au. Figure 12 shows a "Model Mars" where a Å M 0.1 (´- M 3 10 7 ) planet is seen to generate multiple dust gaps over 10 5 orbits. The associated animation shows the gas and dust surface density evolution of this model over 20,000 simulated orbits. The low h/r=0.02 of this model may characterize disk regions within a few au from the star, where terrestrial planets in the Solar System reside. Mars-mass planets are capable of being formed by the streaming instability + pebble accretion (Johansen & Youdin 2007;Johansen & Lambrechts 2017;Ormel 2017;Lin et al. 2018).
If most multi-gap structures at 10-100 au are produced by planet-disk interactions, the ubiquity of disk structures implies a ubiquity of planets analogous to the ice giants (Uranus and Neptune) in our Solar System, and a short time for the Figure 11. Dust surface density map for Model RWI at 1200 orbits in log scale with an aggressive color stretch. The green dot and the green plus symbol mark the locations of the planet and the star, respectively. Two vortices triggered by the RWI form at the outer edge of OG1 and the inner edge of IG1. Two dust clumps at the triangular Lagrange points L 4 and L 5 are also present. In total, there are four dust clumps in this disk with one planet. See Section 5.2 for details. formation of their cores (e.g., HL Tau and GY 91 are believed to be 1 Myr old; ALMA Partnership et al. 2015;Sheehan & Eisner 2018). Microlensing surveys might be able to uncover such a population of "cold Neptunes" (e.g., Poleski et al. 2014). Analyzing microlensing survey data available at the time, Suzuki et al. (2016, Figure 15; see also Zhu et al. 2017) concluded that cold Neptunes are common, with an order unity occurrence rate around main-sequence stars.

Summary
We carry out 2D, two-fluid hydrodynamical simulations of protoplanetary disks with low viscosities (α  10 −4 ) to study the spacings, depths, and total number of dust gaps opened by one sub-thermal mass planet ( 3 ). Our results provide basic guidelines for interpreting multi-gap structures seen in mm continuum emission. Our main findings are: 1. Among the observable properties of gaps and rings, gap location is the least affected by observing conditions, finite optical depth, and time evolution (Sections 3.1 and 3.2; Figure 4). 2. Gap spacings increase with increasing h/r and decreasing M p (Section 3.3 and Figure 5). We provide empirical fitting functions for the double gap OG1+IG1 (defined in Figure 1 and Equation (10)) and IG2 (Equation (11) and Figure 6). 3. Gaps are opened faster by planets with higher masses, and in disks with higher h/r (Section 3.4). 4. More gaps are opened within a given radius range in disks with lower h/r. The number of gaps is less sensitive to M p . 5. The spacings of the double gap in HL Tau (ALMA Partnership et al. 2015) and TW Hya (Andrews et al. 2016), and all three gaps in HD 163296 ) match those of modeled gaps produced by a sub-Saturn mass planet (Section 4; Figures 8 and 9). 6. A low midplane viscosity (α  10 −4 ) is needed for a subthermal mass planet to open multiple gaps ( Figure 10). The Rossby wave instability may develop at the edges of gaps, producing multiple dusty vortices (Figure 11). A planet as low mass as Mars may significantly perturb the dust disk in terrestrial planet-forming regions (Figure 12).   Figure 13. Two-dimensional dust maps for two series of models, showing the evolution of dust gap morphologies with varying h/r (top) and M p (bottom). See Section 3.3 for details.