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

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 in debate. Recently, we documented how one super-Earth planet can open multiple (up to five) dust gaps in a disk with low viscosity ($\alpha\lesssim10^{-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 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. 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,b).
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 real-life gap observations. We study systematically how gap properties depend on planet and disk parameters, and propose generic guidelines connecting ALMA observations with disk-planet interaction models. We introduce our models in §2, and present our main parameter survey of gap properties in §3. These results are applied to three actual disks in §4, discussed in §5, and summarized in §6.

SIMULATIONS
We use two-fluid (gas + dust) 2D hydrodynamics simulations in cylindrical coordinates (radius r and azimuth φ) 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 P = Σ gas c 2 s , where Σ 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: where D dust is the dust diffusivity defined as D dust = ν/(1 + τ 2 stop ) (Takeuchi & Lin 2002), Φ 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), f drag = ΩK τstop (v gas − v dust ) is the drag force on dust by gas, τ stop is the dimensionless stopping time (see below), and f ν is the Shakura-Sunyaev viscous force with ν = αh 2 Ω K , h the disk scale height, Ω K the Keplerian angular frequency, and α < 1 a dimensionless constant. Other symbols have their usual meanings. The last term in Eqn. (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): where r p = 30 AU 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, Fig. 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 p (t) ∝ sin (t) over 10 orbits. In Dong et al. (2017, Fig. 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 Σ gas (r) profile. Dong et al. (2017, Fig. 8) showed that as the background Σ gas (r) profile varies from Σ gas ∝ r 0 to Σ gas ∝ r −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, in the Epstein drag regime, where ρ particle = 1.2 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-1 r p , τ stop is on the order of 0.01. We note that for mm-tocm sized particles often traced by VLA cm continuum observations (possibly having τ stop ∼ 0.1 − 1), 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 φ, and 0.1 to 2.1r p in r. The initial gas disk mass is 0.015M . We set up a wave damping zone at r = 0.1-0.2r p as described in de Val-Borro et al. (2006), and an inflow boundary condition at the outer boundary. The resolution is discussed in §2.1.
We summarize our models in Table 1. The model name indicates the values of h/r and M p /M th adopted (e.g., in Model H003MP02, h/r = 0.03 and M p = 0.2M th ). In all cases the planet mass M p is lower than the disk thermal mass M th = M (h/r) 3 (M = M is the stellar mass); the size of the planet's Hill sphere is comparable to h when M p = M th . An extremely low viscosity is necessary to enable low-mass planets to open gaps. We adopt α = 5×10 −5 for models with M p > 0.1M th , where α is the Shakura & Sunyaev parameter characterizing the disk viscosity; in those models with M p ≤ 0.1M th , 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 0.2M th = 1.8M ⊕ 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),

Models
Purpose h/r  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. 2011b,a), or fail to be seen at all (Bae et al. 2017). Figure 2 compares the radial profiles of Σ dust /Σ dust,0 from five runs of Model H004MP02, labeled A, B, C, D, & 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 φ ). The three runs with n r ≥ 4608 and n φ ≥ 3072 (C, D, and E) achieve visual convergence at r 0.4r 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 φ = 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 0.4r 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 0.3r 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, reminiscent of observed gap depths (e.g., ALMA Partnership et al. 2015;Andrews et al. 2016, of course, observed gap depths are measured in units of surface brightness). Dust gaps in our simulations do not have detectable eccentricities. The locations of gaps and rings, r gap and r ring , are defined as the radii at which Σ dust /Σ dust,0 attains local minima and maxima, respectively; only robust (non-transient) features are counted. 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 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 Σ gas profile, often hard to determine, than gap depths (Dong et al. 2017, Fig. 3). Finally, gap depth and width evolve with time ( §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 ( §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 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 n-th inner ring from the planet), IGn (the n-th inner gap), MR (middle ring), HR (the horseshoe ring), ORn (the n-th outer ring), and OGn (the n-th outer gap). Density fluctuations in gas are much smaller than in dust. See §2 for details. provide a diagnostic. Bae & Zhu (2018b, Fig . 3) showed that more density waves can be discerned in colder disks, or in disks perturbed by lower mass planets.
3.2. The Time Evolution of r gap and r ring Figure 4 shows the time evolution of a representative Model H003MP02. The radial profiles of Σ dust /Σ dust,0 at 700, 1100, 1600, 2600 orbits are plotted. At these times, Σ 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 modelobservation comparisons.
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) ∆ OG1,IG1 increases by 6%, while ∆ 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 A.1 and A.2 in the Appendix list the locations and values of Σ 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).
3.3. The Dependence of Gap Spacing on h/r and M p Figure 5 plots the radial profiles of Σ 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 A.1. All models are shown when surface densities inside gaps are ∼50% of their original values.
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 Σ 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 with M p ∼ 0.1-0.5M th , l sh ∼ 1-2h. We expect the shock locations r p ± l sh 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 Eqn. (9) is only exact in the limit M p M th . Nevertheless, we find the scaling relations implied by Eqn. (9) -l sh ∝ h and l sh ∝ M −2/5 p -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 We expect the agreement to deteriorate as M p approaches M th . The gaps inward of IG1 appear opened by the dissipation of additional (secondary, tertiary, etc.) density waves excited by the planet (Bae et al. 2017). No theoretical calculations quantifying how these waves dissipate have been carried out. We find that r p − r IG2 does not scale with M p and h/r as in Eqns. (9)-(10); instead, the location of IG2 can be fitted by: 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).  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, φ)=(rp, 0) 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 §2.1 for details.
How fast a gap deepens (dΣ dust /dt) depends on a variety of factors. 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).

At fixed M p /M 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 Σ 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 Σ 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 0.4r p and OG1 (inclusive; OG2 is never prominent) in Table 1 (last column). A gap is counted as long as its dΣ dust /dt is consistently negative (i.e., no minimum threshold in Σ 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., 0.4r p ≤ r ≤ r 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 Fig. 4.-The radial profile of Σ dust /Σ dust,0 in Model H003MP02 at four epochs. Local maxima and minima on each profile are marked by dots. Σ 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 §3.2 for details.
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 ≥ 0.4r p , suggesting that the number of gaps may be insensitive to M p (note that Model H003MP005 has only three gaps).

COMPARISON TO REAL DISKS
In this section, we make a quick comparison between our models and observed disks, focussing 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 §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 Eqn. (10), selecting the one model from 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 mm (r) are shown in Figure 8, with modelled planet and gap locations indicated. We mark all modelled gaps at r ≥ 0.4r 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 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, Fig . 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 1M pre-mainsequence star with R = 2.4R and T = 4400K. 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-hour integration times and default settings.  Kwon et al. (2011) andJin et al. (2016) found T midplane ∼ 13 K and 25 K, respectively, by performing radiative transfer modeling of continuum emission. We adopt h/r = 0.05, corresponding to T midplane = 17 K (M HL Tau = 1.7M , Pinte et al. 2016). Eqn. (10) then yields M p 0.8M th (57M ⊕ ). The corresponding Model H005MP08 has four gaps at r ≥ 0.4r p : r IG3 =32 AU, r IG2 =48 AU, r IG1 =64 AU, and r OG1 =78 AU for r p = 71 AU. 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) fit D1, D2, D5, and D6 using a 30M ⊕ planet at 69 AU with h/r = 0.074. Our model better reproduces the observed double gap (∆ D5,D6 0.28 according to Bae et al. 2017), mainly due to our lower h/r. We identify G2 and G3 as a double gap, with ∆ G2,G3 = 0.17. Andrews et al. (2016) found T midplane ∼ 15 K at 45 AU using radiative transfer modeling, corresponding to h/r 0.05 (M TW Hya = 0.88M , Huang et al. 2018). These parameters are essentially identical to those in HL Tau; accordingly, Model H005MP08 is appropriate for TW Hya as well. A planet with M p = 0.8M th (29M ⊕ ), The radial profiles of Σ dust /Σ dust,0 in 7 models with Mp = 0.2M th and successively larger h/r. Bottom: The radial profiles of Σ dust /Σ dust,0 in 5 models with h/r = 0.03 and Mp 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 Mp 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 §3.3 for details. located at r p = 45 AU, opens four gaps at r ≥ 0.4r p : r IG3 = 20 AU, r IG2 = 30 AU, r IG1 = 41 AU, and r OG1 = 50 AU. Figure 9 (middle) shows an image comparison. The model roughly reproduces the spacing of the observed double gap, while G1 finds no ready correspondence with our model.  Fairlamb et al. 2015) in Model H008MP02: r IG2 =54 AU, r IG1 =82 AU, and r OG1 =133 AU for r p = 108 AU. Figure 9 (bottom) shows an image comparison. The adopted h/r = 0.08 is roughly consistent with T midplane ∼ 20 K at r p as inferred by Rosenfeld et al. (2013) and Liu et al.  Figure 5(a). The locations of OG1, IG1, and IG2 are plotted (|rgap−rp|; symbols are the same as in Figure 5), in addition to the distance between OG1 and IG1, r OG1 − r 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 Mp for the 5 models in Figure 5 Figures 5(a) and 6(a) (normalized relative to OG1 with Σ dust /Σ dust,0 = 0.5). Right: The dependence of gap depth on Mp, showing the 5 models in Figures 5(b) and 6(b). As h/r increases, or Mp decreases, the gaps inside rp deepen relative to OG1 (more so for the inner ones). See §3.4 for details.  (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 Imm on dust temperature. The vertical solid line in each panel marks the location of the planet rp in our models, and the vertical dashed lines mark the locations of all modelled gaps at r ≥ 0.4rp (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 Mp = 0.2M th . The location of a third gap IG2 in the models (the innermost vertical dashed line) is labeled too. See §4 for details.  Partnership et al. 2015) and Model H005MP08. The model has been scaled such that rp = 0. 51 (71 AU); Mp = 57M⊕; the elapsed time equals 110 orbits (0.05 Myr); its total flux density is 1 Jy; and the viewing geometry matches the HL Tau disk. The green plus and dot symbol mark the star and the planet, respectively. The beam is indicated in the bottom left corner. The innermost model disk regions are manually masked out, and the color stretch is linear. Middle row: Similar to the top row, but for TW Hya (Andrews et al. 2016). The model (also scaled from H005MP08) has a 29M⊕ planet at 0. 75 (45 AU). Bottom row: Similar to the top row, but for HD 163296 ). The underlying model, H008MP02, has a 65M⊕ planet at 1. 07 (108 AU), and has been run for 700 orbits (0.6 Myr). The color stretch for this row only is logarithmic. The locations of some but not all modelled gaps coincide with those of observed gaps. See §4 for details.
(2018). Note that our identification of the double gap, G2+G3, is ambiguous -it is also possible to identify G1+G2 as the double gap and find a corresponding model, though in that case G3 cannot be accounted for.
A few systems have been found to host two dust gaps in existing observations (e.g., AS 209, Fedele et al. 2018;AA Tau, Loomis et al. 2017;V1094 Sco, van Terwisga et al. 2018and Elias 24, Dipierro et al. 2018). Interpreting each pair as a double gap, we can infer the locations of possible additional, as yet unseen inner gaps. Taking M p = 0.2M th , the two gaps in the Elias 24, AS 209, and AA Tau disks coincide with the double gap opened by a 23M ⊕ planet at 80 AU with h/r = 0.07, a 31M ⊕ planet at 80 AU with h/r = 0.08, and a 19M ⊕ planet at 100 AU with h/r = 0.07, respectively (all distances rescaled using data from the Gaia Collaboration et al. 2018; M Elias 24 = 1M , Wilking et al. 2005; M AS 209 = 0.9M , Tazzari et al. 2016; M AA Tau = 0.85M , Loomis et al. 2017). The location of IG2 in each model is marked in Figure 8. In every case, IG2 falls where I mm has a steep radial gradient, hindering the detection of possible gaps to date. We emphasize that these models are not unique, and they may need to be updated if future observations reveal additional sub-structures. 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 (δΣ gas /Σ gas,0 ∼ 1 − 10%), 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. 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 e ∼ 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.

Low Viscosity and the Rossby Wave Instability
A low viscosity (α 10 −4 ) at the disk midplane (where ∼mm-sized dust particles are expected to settle) is required for a planet with M p M th to open multiple gaps. Figure 10 compares two runs of Model H003MP02 with α = 5×10 −5 (left; standard α) and 3×10 −3 (right). No dust gaps are opened in the latter, while five gaps are opened in the former.
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 Σ 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, private comm.). 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 0.1M ⊕ (3 × 10 −7 M ) planet is seen to generate multiple dust gaps over 10 5 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;Ormel 2017;Johansen & Lambrechts 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 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  Fig. 15; see also Zhu et al. 2017) concluded that cold Neptunes are common, with an order unity occurrence rate around main-sequence stars.
6. SUMMARY We carry out two-dimensional, 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 (M p < M th = M (h/r) 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 ( §3.1; §3.2; Figure 4).
3. Gaps are opened faster by planets with higher masses, and in disks with higher h/r ( §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 modelled gaps produced by a sub-Saturn mass planet ( §4; Figures 8 and 9). 6. A low midplane viscosity (α 10 −4 ) is needed for a sub-thermal 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). Fig. 11.-The 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 4 dust clumps in this disk with one planet. See §5.2 for details.