On the formation of multiple concentric rings and gaps in protoplanetary disks

TL;DR

This paper demonstrates through simulations that multiple gaps in protoplanetary disks can be formed by a single planet's spiral arms, and these gaps' positions can help estimate planetary mass.

Contribution

It reveals that secondary and tertiary spiral arms can also open observable gaps, extending understanding beyond the primary arm in disk-planet interactions.

Findings

Low-mass planets can open multiple gaps over time.

Simulations match observed gaps in HL Tau with a single planet.

Gap positions depend on planetary mass and disk temperature.

Abstract

As spiral waves driven by a planet in a gaseous disk steepen into a shock, they deposit angular momentum, opening a gap in the disk. This has been well studied using both linear theory and numerical simulations, but so far, only for the primary spiral arm -- the one directly attached to the planet. Using two-dimensional hydrodynamic simulations, we show that the secondary and tertiary arms driven by a planet can also open gaps as they steepen into shocks. The depths of the secondary/tertiary gaps in surface density grow with time in a low viscosity disk ($\alpha = 5\times10^{-5}$), so even low-mass planets (e.g., super-Earth or mini-Neptune-mass) embedded in the disk can open multiple observable gaps, provided that sufficient time has passed. Applying our results to the HL Tau disk, we show that a single 30 Earth-mass planet embedded in the ring at 68.8~au (B5) can reasonably well reproduce the positions of the two major gaps at 13.2 and 32.3 au (D1 and D2), and roughly reproduce two other major gaps at 64.2 and 74.7 au (D5 and D6) seen in the mm continuum. The positions of secondary/tertiary gaps are found to be sensitive to the planetary mass and the disk temperature profile, so with accurate observational measurements of the temperature structure the positions of multiple gaps can be used to constrain the mass of the planet. We also comment on the gaps seen in the TW Hya and HD~163296 disk.

Equations

$$\frac{\partial \Sigma }{\partial t}+\nabla \cdot (\Sigma v)=0,$$

$$\Sigma \left(\frac{\partial v}{\partial t}+v\cdot \nabla v\right)=-\nabla P-\Sigma \nabla ({\Phi }_{\ast }+{\Phi }_p).$$

$${\Phi }_p(r,\varphi )=-\frac{G{M}_p}{(|\mathbf{r}-{\mathbf{r}}_{\mathbf{p}}{|}^2+b^2{)}^{1/2}},$$

$$\Sigma (r)={\Sigma }_p{\left(\frac{r}{r_p}\right)}^{-p}$$

$$T(r)=T_p{\left(\frac{r}{r_p}\right)}^{-q},$$

$$\frac{h}{r}={\left(\frac{h}{r}\right)}_p{\left(\frac{r}{r_p}\right)}^{(-q+1)/2}.$$

$$\zeta \equiv \frac{\nabla \times \mathbf{v}}{\Sigma }.$$

$$|l_{\mathrm{sh}}|\approx 0.93{\left(\frac{\gamma +1}{12/5}\frac{M_p}{M_{\mathrm{th}}}\right)}^{-2/5}h,$$

$$\Sigma (r)={\Sigma }_c{\left(\frac{r}{r_c}\right)}^{-\gamma }\exp \left[-{\left(\frac{r}{r_c}\right)}^{2-\gamma }\right]$$

Full text

hydrodynamics, planet-disk interaction, stars: individual (HL Tau, TW Hya, HD 163296)

1 INTRODUCTION

Recent high-resolution observations at (sub)-mm and $\mu$m wavelengths have revealed interesting structures in protoplanetary disks, including rings and gaps (e.g., HL Tau, ALMA Partnership et al. 2015; TW Hya, Andrews et al. 2016, Tsukagoshi et al. 2016, van Boekel et al. 2017; HD 163296, Isella et al. 2016, Zhang et al. 2016; HD 97048, Ginski et al. 2016, Walsh et al. 2016), spiral waves (e.g., MWC 758, Grady et al. 2013, Benisty et al. 2015; SAO 206462, Muto et al. 2012, Garufi et al. 2013, Stolker et al. 2016, AB Aug, Tang et al. 2017), and vortex-like asymmetries (e.g., Oph IRS 48, van der Marel et al. 2013). Among these structures, particularly surprising is rings and gaps, not simply because of their existence but because of their concentric and axisymmetric shape, seeming prevalence among the young protoplanetary disk population (although based on a small number of objects observed with sufficient spatial resolution), and multiplicity in individual disks.

The origin of the observed concentric rings and gaps is yet to be confirmed. Possible mechanisms include gravitational interaction between planets and disks (e.g., Krist et al., 2000; Calvet et al., 2002; van der Marel et al., 2015; Dong et al., 2017), rapid changes in dust properties at the condensation front of volatile species (Zhang et al., 2015; Debes et al., 2016), variation of magnetic activity within disks (Johansen et al., 2009; Ruge et al., 2013; Flock et al., 2015; Pinilla et al., 2016; Ruge et al., 2016) or in MHD disk wind (Suriano et al., 2017), dust evolution (Birnstiel et al., 2015; Okuzumi et al., 2016), photoevaporation (Gorti et al., 2011; Owen et al., 2012), secular gravitational instability (Takahashi & Inutsuka, 2014, 2016), and self-induced dust traps (Gonzalez et al., 2015, 2017).

Assuming planets are responsible for the disk gaps, numerical simulations generally use as many planets as the number of gaps to explain multi-ring/gap structures (e.g., Dipierro et al. 2015, Isella et al. 2016, Jin et al. 2016). In some simulations, however, it appears that multiple gaps and rings can form from one planet (see for example Figure 1 and 4 of Zhu et al. 2014, and Figure 5 of Bae et al. 2016c). Interestingly, the additional gaps open far from the planet, $\sim 0.6~{}r_{p}$ in Zhu et al. (2014) and $\sim 0.45$ and $0.6~{}r_{p}$ in Bae et al. (2016c), where $r_{p}$ denotes the planet’s semi-major axis. Dong et al. (2017) recently showed that an observable “double gap” can be opened at the vicinity of a super-Earth mass planet, about a disk scale height from the planet, in a low-viscosity disk. In addition to the “double gap” adjacent to the planet, in Dong et al. (2017) an additional gap opens at $\sim 0.6~{}r_{p}$ which seems to share the same origin with the ones seen in Zhu et al. (2014) and Bae et al. (2016c). The detailed cause of the additional sets of rings and gaps that form far from a planet, however, has not yet been investigated in the literature.

In this paper, we study the origin of multiple rings and gaps created by a single planet. When a planet forms, the gravitational interaction between the planet and its host disk excites spiral arms. Even for a sufficiently low-mass planet, for which the excitation and initial propagation of spiral density waves launched by the planet can be explained well by linear theory, it is known that spiral waves non-linearly steepen to spiral shocks as the waves propagate (Goodman & Rafikov, 2001). When a spiral wave shocks, angular momentum is transferred to the disk, which in turn leads the redistribution of disk material, opening a gap (Rafikov, 2002). This phenomenon has been well studied using both linear theory (Goodman & Rafikov, 2001; Rafikov, 2002) and numerical simulations (Muto et al., 2010; Dong et al., 2011; Duffell & MacFayden, 2012; Zhu et al., 2013; Dong et al., 2017), but so far, only for the primary arm – the one directly attached to the planet.

Previous protoplanetary disk simulations have shown that a secondary spiral arm can be excited by a planet (e.g., Kley, 1999; Fung & Dong, 2015; Juhász et al., 2015; Zhu et al., 2015; Bae et al., 2016a, c; Lee, 2016). In addition, some simulations have revealed that a tertiary arm can also be excited in the inner disk (Zhu et al., 2015; Bae et al., 2016c). We will present the mechanism for multiple spiral arms in a forthcoming paper (Bae & Zhu 2017, in preparation). In the present work, we show that secondary and tertiary spiral arms excited by a planet can create secondary and tertiary gaps in the inner disk when they shock, in the same way that the primary arm opens a gap. We propose that spiral shocks associated with secondary/tertiary arms can be a possible origin for the observed multiple sets of rings and gaps in some protoplanetary disks, as we will show for the HL Tau disk as an example later in the paper.

This paper is organized as follows. In Section 2 we introduce our numerical setup. In Section 3 we present our results, describing the formation of multiple rings and gaps from a low-mass planet in Section 3.1 and a high-mass planet in Section 3.2. We then study a long-term evolution of the rings and gaps in Section 4. In Section 5, applying our results we explain the multiple rings and gaps seen in the HL Tau disk using a single planet. Finally, in Section 6 we discuss implications of the present work and comment on some individual disks with multiple rings and gaps.

2 NUMERICAL METHODS

We solve the hydrodynamic equations for mass and momentum conservation in two-dimensional polar coordinates $(r,\phi)$ using FARGO 3D (Benítez-Llambay & Masset, 2016):

[TABLE]

[TABLE]

In the above equations, $\Sigma$ is the gas surface density, $v$ is the velocity, $P=\Sigma c_{s}^{2}$ is the pressure where $c_{s}$ is the isothermal sound speed, and $\Phi_{*}=GM_{*}/r$ is the gravitational potential from the central star. The planetary potential $\Phi_{p}$ is

[TABLE]

where $M_{p}$ is the planetary mass, ${\bf r}$ and ${\bf r_{p}}$ are the radius vectors of the center of grid cells in question and of the planet, and $b$ is the smoothing length which is assumed to be $10\%$ of the disk scale height at $r=r_{p}$. To provide a clearer picture, we ignore the indirect term that arises from the offset between the central object and the origin of the coordinate system; however, we confirmed that including the indirect term does not affect our conclusion.

We increase the planetary mass over a period time $t_{\rm growth}$ at the beginning of each simulation in order to avoid sudden perturbation to the disk, following $M_{p}(t)\propto\sin[(\pi/2){\rm min}(t/t_{\rm growth},1)]$. The growth timescale $t_{\rm growth}$ is set to $t_{\rm growth}=20~{}t_{p}\times(M_{p}/0.1M_{\rm th})$ for various planetary masses, with a maximum of $100~{}t_{p}$, where $t_{p}=2\pi/\Omega_{p}$ is the planet’s orbital time. In the above growth timescale equation, $M_{\rm th}\equiv c_{s}^{3}/\Omega G=M_{*}(h/r)_{p}^{3}$ is the so-called thermal mass (Lin & Papaloizou, 1993; Goodman & Rafikov, 2001), at which planetary mass the Hill radius of the planet is comparable to the disk scale height. Below the thermal mass, the excitation and the initial propagation of the spiral density wave from a planet can be well approximated in the linear regime (e.g., Goldreich & Tremaine, 1980). When planetary mass becomes greater than the thermal mass, however, the spiral density wave becomes non-linear and shock background disk gas as soon as they launch (Goodman & Rafikov, 2001).

Our initial disk has power-law surface density and temperature distributions:

[TABLE]

and

[TABLE]

where $\Sigma_{p}$ and $T_{p}$ are the surface density and temperature at the location of the planet $r=r_{p}$. The above temperature profile corresponds to the disk aspect ratio profile of

[TABLE]

We adopt $p=1$, $q=0$, and $(h/r)_{p}=0.07$ in our fiducial model.

While the assumed constant temperature ($q=0$) may not be a very realistic temperature structure for global protoplanetary disks, this enables us to directly compare our numerical results with predictions for spiral wave excitation and propagation from linear theory that is developed assuming a constant temperature (e.g., Goodman & Rafikov, 2001). However, we note that the physics behind the gap opening through spiral shocks should remain the same even for different power-law temperature slopes. We add a small viscosity $\alpha=5\times 10^{-5}$ to the gas, where $\alpha$ denotes a Shakura-Sunyaev viscosity parameter (Shakura & Sunyaev, 1973): $\nu=\alpha c_{s}^{2}/\Omega$. As viscosity tends to smooth out structure, we initially choose this low value of $\alpha$ to show the mechanism of gap opening more clearly and to compare the numerical results more directly with linear theory. The effects of larger viscosity are discussed in Section 6.

Our simulation domain extends from $r_{\rm in}=0.2~{}r_{p}$ to $r_{\rm out}=2.0~{}r_{p}$ in radius and from 0 to $2\pi$ in azimuth. Because we use an isothermal equation of state and do not include disk self-gravity, no physical units are necessary and the quantities presented in this work can be scaled. For our fiducial calculations presented in Section 3, we adopt 4096 logarithmically-spaced grid cells in the radial direction and 11160 uniformly-spaced grid cells in the azimuthal direction. For the long-term calculations in Section 4 we use $2048\times 5580$ grid cells. With both choices all grid cells are nearly square-shaped ($\Delta r:r\Delta\phi\simeq 1:1$). We have tested with various numerical resolution and found that this very high spatial resolution is required to properly capture the opening of secondary and tertiary gaps. For the model disk assumed in this work, using $2048\times 5580$ or more grid cells can resolve the spiral shock at all radii but using $1024\times 2290$ grid cells fails to resolve the shock in the inner disk $r\lesssim 0.4~{}r_{p}$. When a spiral shock is under-resolved we find that no gap associated with the spiral arm forms.

At the radial boundaries, we adopt a wave-damping zone (de Val-Borro et al., 2006) to reduce wave reflection. We have tested with other boundary conditions, including outflow boundary condition and reflecting boundary condition, and confirmed that the choice of boundary condition does not affect our conclusion that secondary/tertiary gaps open when secondary/tertiary arms steepen into shocks.

3 Formation of Rings and Gaps

3.1 Results with a $0.1~{}M_{\rm th}$ planet

We start with results for a $0.1~{}M_{\rm th}$ planet (corresponding to $11.4~{}M_{\earth}$ with $M_{*}=1M_{\odot}$ and $(h/r)_{p}=0.07$), for which the excitation and the initial propagation of the spiral density wave can be explained with linear theory. In Figure 1 (a), we show the perturbed density distribution $\delta\Sigma/\Sigma_{0}$ at $t=50~{}t_{p}$, where $\Sigma_{0}$ is the initial density at each position of the disk and $\delta\Sigma=\Sigma-\Sigma_{0}$. The planet launches a trailing spiral arm that is directly attached to the planet (i.e., primary arm). The planet starts to launch a secondary arm around $r=0.6~{}r_{p}$. The secondary arm is closely located to the primary in azimuth, only $\sim 50^{\circ}$ apart. There is a hint of another arm at $r\sim 0.3-0.4~{}r_{p}$. However, this tertiary arm-like structure does not produce a shock within our simulation domain. We have tested this with a factor of two higher resolution ($8192\times 22320$ grid cells) but no shock was found to be produced from this arm-like structure. We find that the perturbation from this tertiary arm-like structure gradually increases but also excites closer to the planet when we increase the planetary mass, as we discuss in the following section.

In order to diagnose the radial location at which a planet-driven spiral wave becomes a shock, we compute the potential vorticity (PV) $\zeta$, defined as

[TABLE]

For a barotropic fluid as we have in our models, it is well known that the PV is conserved in the absence of shocks. When a shock is present, on the other hand, the PV experiences a jump at the shock front (Li et al., 2005; Dong et al., 2011). In Figure 1 (b), we show the perturbed PV distribution $\delta\zeta/\zeta_{0}$, where $\delta\zeta=\zeta-\zeta_{0}$ and $\zeta_{0}=(\nabla\times{\bf v}_{0})/\Sigma_{0}$ is the unperturbed PV at the beginning of the simulation. Looking at the perturbed PV distribution along the primary arm first, one can see that it does not produce a shock at the immediate vicinity of the planet but a few scale heights away from the planet as expected from linear theory. Linear waves from a low-mass planet ($M_{p}\ll M_{\rm th}$) steepen into shocks as they propagate away from the planet (Goodman & Rafikov, 2001) and the distance at which linear waves become a shock can be approximated from linear theory (Goodman & Rafikov, 2001):

[TABLE]

where $\gamma$ is the adiabatic index and $h$ is the disk scale height. With $M_{p}=0.1~{}M_{\rm th}$ and $(h/r)_{p}=0.07$, $l_{\rm sh}\approx\pm 2.5h\approx\pm 0.176$, so our numerical results show a good agreement with the prediction from linear theory. The PV perturbation accumulates over time, creating two rings – one interior and the other exterior to the planetary orbit – of high PV region around the radii a spiral becomes a shock ($r\sim 0.84$ and $1.17~{}r_{p}$).

We find that the same happens for the secondary arm. The secondary arm starts to shock background medium around $r=0.5~{}r_{p}$, producing a PV jump at its shock front and creating a ring of high PV.

In Figure 1 (c) and (d), we plot the azimuthally averaged profiles of the perturbed density and perturbed PV. As can be seen from the figures, the locations where spiral arms produce shocks and the locations of density deficits match very well with each other. Again, no apparent shock is seen at the tertiary arm-like structure so no gap forms from it.

In the outer disk, additional spiral arms are not excited within our simulation domain ($r_{\rm out}=2~{}r_{p}$), so only one gap associated with the primary arm opens at $r=1.17~{}r_{p}$. As an aside, for the planetary mass and disk temperature profile considered here ($0.1~{}M_{\rm th}$ and $(h/r)_{p}=0.07$), it is unlikely that additional spiral arms form even beyond $r=2~{}r_{p}$ (Bae & Zhu 2017, in preparation).

3.2 Results with a $3~{}M_{\rm th}$ planet

We now describe results for a more massive planet, $3~{}M_{\rm th}$ (corresponding to $1.03~{}M_{\rm Jup}$ with $M_{*}=1M_{\odot}$ and $(h/r)_{p}=0.07$). This mass is chosen such that three spiral arms and the associated gaps develop well within our simulation domain. As we gradually increase the planetary mass from $0.1~{}M_{\rm th}$ to $3~{}M_{\rm th}$, we find that the secondary and tertiary arms are excited in the inner disk closer to the planet – this is much more prominent for the tertiary as we show below. Increasing planetary mass also increases the surface density perturbations at the shock front and the separation between the arms.

In Figure 2 (a), we plot the perturbed density distribution. Three arms are launched by a $3~{}M_{\rm th}$ planet: a primary arm that is attached to the planet, a secondary arm launched around $r=0.6~{}r_{p}$, and a tertiary arm launched around $r=0.5~{}r_{p}$. While the secondary arm launches at a similar radius to the low-mass planet case, the tertiary arm starts to appear closer to the planet, presumably because non-linear effects shift the launching radius of the tertiary arm more than that of the secondary arm (Bae & Zhu 2017, in preparation). Also, the arm-to-arm separation increases as the planetary mass increases, consistent with previous studies (e.g., Zhu et al. 2015; Fung & Dong 2015; Lee 2016). At $r=0.6~{}r_{p}$, the secondary arm is offset by $112^{\circ}$ from the primary arm, as opposed to $50^{\circ}$ in the low-mass planet case. Another noticeable difference is that the planet opens a deep gap around itself and therefore it is no more embedded in the disk.

In Figure 2 (b) we plot the perturbed PV distribution to find the shock locations. The secondary and tertiary arms start to shock at $r\sim 0.56$ and $0.37~{}r_{p}$. The arms accumulate PV as they shock, creating rings of high PV. The azimuthally averaged distributions of the perturbed density and PV are presented in Figure 2 (c) and (d), respectively. As in the case of the low-mass planet, the locations of density gaps and the local PV maxima match well with each other, supporting that the idea of spiral shocks opening gaps is also valid for more massive gap-opening planets.

4 Long-term Evolution of Rings and Gaps

One of the interesting questions regarding the secondary/tertiary gaps is how deep could they be. To address the question, we ran simulations for 1000 planetary orbits and measured the depths of gaps for various planetary masses: 0.1, 0.3, 0.5, 1, 2, and $3~{}M_{\rm th}$, corresponding to $11.4~{}M_{\earth},0.10~{}M_{\rm Jup},0.17~{}M_{\rm Jup},0.34~{}M_{\rm Jup},0.69~{}M_{\rm Jup}$, and $1.03~{}M_{\rm Jup}$ when $M_{*}=1M_{\odot}$ and $(h/r)_{p}=0.07$ are assumed.

We computed the gap depth $\Delta\Sigma_{\rm gap}$ as follows. We first generated a power-law fit for the adjacent density rings to a gap, illustrated with dotted lines in Figure 3 (a). For the inner and outer primary gaps (indicated with red and blue arrows in the figure), we generated the fit using the density rings at the inner and outer gap edges, instead of the local density maximum at the corotation region ($r\sim 1~{}r_{p}$). We then calculated the expected surface density at the gap center (i.e., local density minimum) using the power-law fit. The gap depth is defined as the ratio between the power-law surface density at the gap center from the fit and the actual surface density there. This ratio should be unity in a power-law disk when there is no gap, and gradually increase as a gap opens.

The gap depth calculated following this method is plotted in Figure 3 (b) and (c). Looking first at planets with $M_{p}\leq 1M_{\rm th}$, neither the primary nor the secondary gaps reach a steady state until $1000~{}t_{p}$. We find that the two primary gaps have a comparable depth, although the gap exterior to the planetary orbit is slightly deeper than the gap interior to the planetary orbit in all cases. The secondary gaps are shallower than the primary gaps by a few percent ($0.1~{}M_{\rm th}$) to a factor of $\sim 3$ ($1~{}M_{\rm th}$) in terms of $\Delta\Sigma_{\rm gap}$ value. Secondary gaps are shallower than primary gaps presumably because secondary arms are weaker than primary arms when compared at their shock locations.

For the $3~{}M_{\rm th}$ case, we find that the gap depths reach to a quasi-steady state as the Rossby wave instability (RWI) is triggered at the rings. This happens at $\sim 200~{}t_{p}$ for the primary rings, $\sim 500~{}t_{p}$ for the secondary ring, and $\sim 600~{}t_{p}$ for the tertiary ring. On the other hand, as inferred from Figure 3 (a) and (c), the RWI is not able to get rid of the rings/gaps but it only impedes further growth of the rings/gaps. It is also interesting to note that the tertiary gap can be as deep as or even deeper than the primary/secondary gaps. We conjecture that this is because (1) the tertiary arm can be as strong as the secondary arm at least locally in radius (Bae & Zhu 2017, in preparation) and/or (2) the gaps are (partially) filled back via the RWI in such a way that the resulting tertiary gap depth is comparable to or deeper than the primary/secondary gaps.

In Figure 4 (a) we show the surface density distribution at $t=1000~{}t_{p}$ for the $3~{}M_{\rm th}$ case. Multiple sets of rings and gaps are immediately noticeable: one gap exterior to the planetary orbit opened by the primary arm, three gaps interior to the planetary orbit each of which is opened by the primary, secondary, and tertiary arms (from large to small radii). As depicted in the figure all the gaps and rings are concentric, similar to the ones seen in recent high resolution observations. In Figure 4 (b) we display the density perturbation above/below the azimuthal average at each radius in order to highlight any non-axisymmetric structure present in the disk. One can see that the rings are nearly axisymmetric, which is again similar to the observed multiple rings. Note also that the radial locations where the PV experiences a jump (indicated by the arrows in the figure) match well with the gap locations.

5 Parameter Study and Possible Application to HL Tau

Now that we have shown that multiple gaps can be opened by multiple spiral waves from a single planet, we explore the dependence of additional gap positions on the background disk temperature profile and the planetary mass, and consider a possible application of our results to the case of HL Tau. Our goal here is not to produce a detailed model of the disk but to explore whether some of the dust gaps might be explained with a single planet.

The dust emission of the HL Tau disk shows remarkable multiple, concentric rings and gaps (ALMA Partnership et al., 2015). The ALMA observation revealed four major gaps at 13.2 (D1), 32.3 (D2), 64.2 (D5), and 73.7 au (D6), and rings in between (ALMA Partnership et al., 2015). Also, there are additional narrower and shallower rings and gaps (B3/D3, B4/D4, and B7/D7; ALMA Partnership et al. 2015).

To model the HL Tau disk, we adopt the gas surface density profile constrained by the mm continuum observations of Kwon et al. (2011):

[TABLE]

where $\gamma=-0.23$ and $r_{c}=79$ au. For the temperature profile, we make use of the brightness temperature obtained in ALMA Partnership et al. (2015): using 59 K at 20 au and 24 K at 81 au we obtain the temperature power-law index $q=0.65$111Since the HL Tau disk is likely to be (marginally) optically thick at 1.3 mm (ALMA Partnership et al., 2015; Carrasco-González et al., 2016; Pinte et al., 2016), this temperature profile may not be a completely accurate estimate of the midplane temperature.. Assuming a central stellar mass of $1.3~{}M_{\odot}$ (ALMA Partnership et al., 2015) and a mean molecular weight of 2.4, the disk aspect ratio profile follows $h/r=0.074(r/r_{p})^{0.175}$, consistent with irradiation heating by the central object in a moderately-flared disk. For the simulations presented in this section, we adopt $\alpha=3\times 10^{-4}$ as estimated for the HL Tau disk by Pinte et al. (2016). We adopt $r_{\rm in}=0.05~{}r_{p}$ and $r_{\rm out}=2~{}r_{p}$ and use $4096\times 6974$ grid cells.

Assuming that the planet has a circular orbit and does not migrate, we place a planet at B5 (68.8 au; ALMA Partnership et al. 2015), the ring in between D5 and D6, which is also suggested by the “double gap” feature in Dong et al. (2017). We then vary the planetary mass and search for the one that opens secondary and tertiary gaps at the locations of the observed D1 and D2. Since we have only the gas component in our simulations while the ALMA continuum emission traces solid particles, we do not attempt to reproduce the depth of each gap in the observed intensity profile. As we show below, we note that the planetary mass that best reproduces the two innermost gaps (D1 and D2) does not result in a perfect match for the gaps near the planet (D5 and D6). However, we do not aim to reproduce the locations of D5 and D6, since the gas flow near the planetary orbit can be extremely complicated (e.g., Fung et al., 2015; Fung & Chiang, 2016) and thus the dust distribution in the corotation region, which produces the observed emission, cannot be accurately predicted with gas-only calculations. More thorough comparisons between numerical models and the observed data will require multi-fluid calculations that can simultaneously compute gas and dust dynamics; we defer this to a future study. On the other hand, we note that the locations of the secondary/tertiary gaps do not evolve over time when the temperature profile and the location of the planet is fixed as we assumed, and the gap depths are likely to grow over time as Figure 3 suggests.

Based on a suite of calculations varying the planetary mass we find that the two innermost gaps can be reasonably well reproduced with a 30 Earth-mass planet at B5. In Figure 5 (a) we present the azimuthally averaged perturbed density distribution, together with the normalized intensity profile of the ALMA continuum observation at Band 6 ALMA Partnership et al. (2015). The planet opens three gaps interior to its orbit, each of which is associated with the primary, secondary, and tertiary spiral arm (from large to small radii), and one gap exterior to its orbit. As shown, with the 30 Earth-mass planet at 68.8 au the locations of the secondary and tertiary gaps from our simulation match very well with the observed gaps at 13.2 and 32.3 au.

In Figure 5 (b), we present the perturbed density distributions with 20, 30, and 40 Earth-mass planets to show how a change in planetary mass affects the locations of secondary/tertiary gaps. As seen in the figure, the secondary/tertiary gaps open closer to the planet with a larger mass planet, and open farther from the planet with a lower mass planet. This is because the stronger the spiral arms are, they start to shock disk gas closer to the planet.

Another important factor that determines secondary/tertiary gap depths is the background disk temperature profile since it governs the propagation of waves as well as wavelength of the spiral arms. In Figure 5 (c), we present the perturbed density distributions with $20\%$ higher ($(h/r)_{p}=0.081$) and $20\%$ lower temperatures ($(h/r)_{p}=0.066$) than our fiducial choice ($(h/r)_{p}=0.074$). Note that even with a $\pm 20\%$ of change in temperature results in a noticeable change in the gap locations – the secondary gap location changes by about $\pm 4$ au, which is about the linear resolution for the Band 6 ALMA continuum data, and the tertiary gap location changes by about $\pm 2$ au. We also test with different temperature power-law slopes, $q=0.8$ and $q=0.5$, fixing the temperature at the planet’s location. As shown in Figure 5 (d), we find that the temperature slope can also affect the gap locations although the effect may not be as significant as changing the background temperature.

Because of this degeneracy in gap locations, there are other possible sets of planetary masses and disk temperature profiles that can produce similar ring/gap positions as seen in the HL Tau disk. For example, when we increase the disk temperature by $20~{}\%$ we find that the two gap locations can be well reproduced by a 60 Earth-mass planet (see the black dashed curve in Figure 5 (a)). Having more accurate observational measurements of the disk temperature structure and other disk properties, such as the disk surface density and viscosity which can affect spiral wave excitation and propagation, will therefore help constrain the mass of the currently unseen planet. Also, numerical simulations that can simultaneously compute gas and dust dynamics are warranted as the present calculation simulates only the gas component.

With only one planet we were not able to reproduce finer structures (e.g., B3/D3, B4/D4, B7/D7 in ALMA Partnership et al. 2015). It is interesting to speculate that two additional newly formed planets in the rings B3 ($\sim 47$ au) and B7 ($\sim 97$ au) have just started to open the shallow gaps, possibly triggered by the enhancement of solid population at the local pressure maxima created by a hypothesized planet at B5 (and the resulting secondary gap).

6 DISCUSSION

Observability of the secondary/tertiary gaps. In general, $\sim 10~{}\%$ of local density enhancement – corresponding approximately to $\Delta\Sigma_{\rm gap}\sim 1.2$ – is known to be able to induce trapping of solid particles in the density rings (e.g., Pinilla et al., 2012; Zhu et al., 2014; Dong et al., 2017), although the exact number should be dependent upon the unperturbed density and temperature profiles. As shown in Figure 3, we find that all secondary/tertiary gaps but the $0.1~{}M_{\rm th}$ case have $\Delta\Sigma_{\rm gap}\geq 1.2$ after 1000 planetary orbits. This suggests that even for the planets with $M_{p}<1M_{\rm th}$, sub-Jovian mass in typical protoplanetary disks, the secondary and tertiary rings may collect enough solid particles to produce observable signatures in (sub-)mm continuum, when sufficient time has passed. It is also possible that the secondary and tertiary gaps are detectable in gas, or scattered light which traces small $\mu$m-sized particles that are well mixed with gas, for high-mass planets. For instance, Yen et al. (2016) have detected two gaps in HCO+ emission from the HL Tau disk with gap depths, measured in a similar way to ours, of $\sim 2.4$ and $\sim 5.0$. As our $3~{}M_{\rm th}$ case in Figure 3 (c) suggests, the secondary and tertiary gap depths for massive planets can be well above the detected gap depths in HL Tau.

Planet formation. Solid particles (with appropriate sizes) can be efficiently collected in density rings via aerodynamic drag (e.g., Pinilla et al., 2012). The concentration of particles in the multiple rings formed by planets will provide a favorable condition to expedite the growth of solid bodies within the bumps through coagulation between solid particles, gravitational instability (Youdin & Shu, 2002), the streaming instability (Youdin & Goodman, 2005), and more efficient pebble accretion onto planetesimals (Johansen & Lacerda, 2010; Ormel & Klahr, 2010; Lambrechts & Johansen, 2012; Morbidelli & Nesvorny, 2012). It can also prevent solid particles from rapidly drifting inward, which would happen in a disk with a smooth radial density distribution, thus potentially overcoming the radial drift barrier (Whipple, 1972; Weidenschilling, 1977). We also note that when the first planet forms in a disk, it can have an influence on the position of subsequent planet formation in the disk.

Three dimensional effects. In three dimensions, it is possible that planet-driven spiral waves can become unstable to the spiral wave instability (SWI; Bae et al., 2016b, c). On the other hand, the three-dimensional simulations in Bae et al. (2016c) indicate that the formation of secondary/tertiary rings and gaps are not completely suppressed even when the SWI is in action – in their simulation the three gaps open in the inner disk at $R\sim 0.85,0.6$, and 0.45 (see their Figure 5) and these correspond to the radii at which each of the three arms start to become apparent in the density plot (see their Figure 6). It is likely though that the rings and gaps become less prominent than when the SWI is absent.

Outer disk. With a planetary mass range and the disk temperature profiles considered in this work, planets excite only one arm in the outer disk. On the other hand, more than one spiral arm can excite in the outer disk depending on the planetary mass and the disk temperature profile (Bae & Zhu 2017, in preparation), in which case more than one gap can open in the outer disk. Also, as noted in Zhu et al. (2014), it is possible that the vortex at the primary ring in the outer disk can excite density waves. The waves then can steepen to shocks, leading to the formation of additional sets of rings/gaps at larger radii.

Effect of viscosity. The fiducial models in Section 3 used $\alpha=5\times 10^{-5}$. To test the effect of larger viscosity, we adopt $\alpha=5\times 10^{-4}$ and $5\times 10^{-3}$ for the fiducial model presented in Section 3.2 ($3~{}M_{\rm th}$ planet case). As shown in Figure 6, we find that the secondary and tertiary gaps still open with $\alpha=5\times 10^{-4}$ but their formation is suppressed with $\alpha=5\times 10^{-3}$. If the multiple concentric gaps observed in young protoplanetary disks are due to spiral shocks driven by planets, the existence of those gaps might thus imply low viscosities in the disks ($\alpha\lesssim 10^{-3}$), in potential support of recent theoretical (e.g., Bai & Stone, 2013; Gressel et al., 2015; Hartmann & Bae, 2017) and observational studies (e.g., Flaherty et al., 2015, 2017; Pinte et al., 2016; Mulders et al., 2017). We also note that the tertiary gap forms at a slightly smaller radius with $\alpha=5\times 10^{-4}$ than with $\alpha=5\times 10^{-5}$, suggesting that having accurate observational measurments on the disk viscosity is also desired to determine the mass of unseen planets. A more complete exploration of the parameter space is out of the scope of the paper and will be presented in a future paper; quantitatively, however, additional gaps can open when the “effective” viscosity (e.g., Goodman & Rafikov, 2001) arising from spiral shocks is larger than the background disk viscosity.

Orbital eccentricity and migration. When a planet has a non-zero orbital eccentricity, it will give a rise to new family of spiral waves having differing pattern speeds around the mean motion of the planet (Goldreich & Tremaine, 1980). As these spiral waves will become shocks at different radial locations, additional gaps will become broader and shallower. Orbital migration can also affect the properties of the gaps (e.g., location, depth) when migration is fast (e.g., Malik et al., 2015). In particular, the gap depths at any instance can be shallower as the gap-opening occurs over an extended period of time (see Figure 3).

Comments on individual disks with multiple gaps. We have shown HL Tau as a possible example of a single planet opening multiple sets of rings/gaps. To help facilitate future numerical modeling of the observed disks with multiple gaps, here we comment on some individual disks. For the TW Hya disk, it may well be that a low-mass planet ($<1M_{\rm th}$) embedded in the ring at 40 au opens the double gap at 37 and 43 au via the primary spiral shock, as suggested by Dong et al. (2017), and opens a secondary gap at 22 au through a shock from the secondary arm. This resembles our $0.1~{}M_{\rm th}$ model. For the HD 163296 disk, presumably a more massive planet ($\gtrsim 1M_{\rm th}$) opens a broad gap around 160 au and the secondary and tertiary spiral arms open two gaps in the inner disk at 60 and 100 au, similar to what is seen in our $3~{}M_{\rm th}$ model.

Future work. In a forthcoming paper (Bae & Zhu 2017, in preparation), we discuss the origin of multiple spiral arms driven by a single planet. Better understanding of the formation of secondary/tertiary arms will enable comparison between linear theory and numerical results. It will also help make predictions for secondary/tertiary gap locations that can be used as a prescribed formula, for example, in the pebble accretion calculations where the evolution of the gas surface density often relies on simplified viscous evolution (plus a prescribed primary gap formula in some models). Applications to individual disks with multiple rings/gaps can shed light on the existence of planets in the disks, but also provide critical insights into studies of planet formation and planet-disk interactions.

JB is grateful to Ke Zhang for providing the radial intensity profile of the ALMA HL Tau data and for helpful discussion. This work was supported in part by NASA grant NNX17AE31G. We acknowledge the following: computational resources and services provided by Advanced Research Computing at the University of Michigan, Ann Arbor; the XStream computational resource, supported by the National Science Foundation Major Research Instrumentation program (ACI-1429830); the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation grant number ACI-1548562; and the NASA High-End Computing (HEC) Program through the NASA Advanced Supercomputing (NAS) Division at Ames Research Center. We made use of the following ALMA data: ADS/JAO.ALMA $\#$2011.0.00015.SV. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), NSC and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc.