On the Formation of Multiple Concentric Rings and Gaps in Protoplanetary Disks

Recent high-resolution observations at (sub)-mm and μm wavelengths 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, 2016; Flock et al. 2015; Pinilla 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, Figures 1 and 4 of Zhu et al. 2014, and Figure 5 of Bae et al. 2016a). Interestingly, the additional gaps open far from the planet, $\sim 0.6\,{r}_{p}$ in Zhu et al. (2014) and ∼0.45 and $0.6\,{r}_{p}$ in Bae et al. (2016a), where r_p denotes the planet's semimajor 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 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. (2016a). 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 nonlinearly 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, 2017; Duffell & MacFayden 2012; Zhu et al. 2013), but so far, only for the primary arm, which is 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, 2016c; 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. 2016a). We will present the mechanism for multiple spiral arms in a forthcoming paper (J. Bae & Z. 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.

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

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\rm{\Sigma }}}{\partial t}+{\rm{\nabla }}\cdot ({\rm{\Sigma }}v)=0,\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&{\rm{\Sigma }}\left(\displaystyle \frac{\partial v}{\partial t}+v\cdot {\rm{\nabla }}v\right)=-{\rm{\nabla }}P-{\rm{\Sigma }}{\rm{\nabla }}({{\rm{\Phi }}}_{* }+{{\rm{\Phi }}}_{p}).\end{eqnarray} \tag{ 2 }$$

In the above equations, Σ is the gas surface density, v is the velocity, $P={\rm{\Sigma }}{c}_{s}^{2}$ is the pressure where c_s is the isothermal sound speed, and ${{\rm{\Phi }}}_{* }={{GM}}_{* }/r$ is the gravitational potential from the central star. The planetary potential ${{\rm{\Phi }}}_{p}$ is

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{p}({\boldsymbol{r}},\phi )=-\displaystyle \frac{{{GM}}_{p}}{{(| {\boldsymbol{r}}-{{\boldsymbol{r}}}_{p}{| }^{2}+{b}^{2})}^{1/2}},\end{eqnarray} \tag{ 3 }$$

where M_p is the planetary mass, ${\boldsymbol{r}}$ and ${{\boldsymbol{r}}}_{{\boldsymbol{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 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}_{\mathrm{growth}}$ at the beginning of each simulation in order to avoid sudden perturbation to the disk, following ${M}_{p}(t)\propto \sin [(\pi /2)\min (t/{t}_{\mathrm{growth}},1)]$. The growth timescale ${t}_{\mathrm{growth}}$ is set to ${t}_{\mathrm{growth}}=20\,{t}_{p}\times ({M}_{p}/0.1{M}_{\mathrm{th}})$ for various planetary masses, with a maximum of $100\,{t}_{p}$, where ${t}_{p}=2\pi /{{\rm{\Omega }}}_{p}$ is the planet's orbital time. In the above growth timescale equation, ${M}_{\mathrm{th}}\equiv {c}_{s}^{3}/{\rm{\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 with the 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 nonlinear and shock background disk gas as soon as they launch (Goodman & Rafikov 2001).

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

$$\begin{eqnarray}&&{\rm{\Sigma }}(r)={{\rm{\Sigma }}}_{p}{\left(\displaystyle \frac{r}{{r}_{p}}\right)}^{-p}\end{eqnarray} \tag{ 4 }$$

and

$$\begin{eqnarray}&&T(r)={T}_{p}{\left(\displaystyle \frac{r}{{r}_{p}}\right)}^{-q},\end{eqnarray} \tag{ 5 }$$

where ${{\rm{\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

$$\begin{eqnarray}&&\displaystyle \frac{h}{r}={\left(\displaystyle \frac{h}{r}\right)}_{p}{\left(\displaystyle \frac{r}{{r}_{p}}\right)}^{(-q+1)/2}.\end{eqnarray} \tag{ 6 }$$

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 α denotes a Shakura–Sunyaev viscosity parameter (Shakura & Sunyaev 1973): $\nu =\alpha {c}_{s}^{2}/{\rm{\Omega }}$. As viscosity tends to smooth out structure, we initially choose this low value of α 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}_{\mathrm{in}}=0.2\,{r}_{p}$ to ${r}_{\mathrm{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 11,160 uniformly spaced grid cells in the azimuthal direction. For the long-term calculations in Section 4, we use 2048 × 5580 grid cells. With both choices, all grid cells are nearly square-shaped (${\rm{\Delta }}r\,:r{\rm{\Delta }}\phi \simeq {\rm{1\; :\; 1}}$). We 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 × 5580 or more grid cells can resolve the spiral shock at all radii, but using 1024 × 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 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.1. Results with a 0.1 ${M}_{{th}}$ Planet

We start with results for a $0.1\,{M}_{\mathrm{th}}$ planet (corresponding to $11.4\,{M}_{\oplus }$ with ${M}_{* }=1\,{M}_{\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 {\rm{\Sigma }}/{{\rm{\Sigma }}}_{0}$ at $t=50\,{t}_{p}$, where ${{\rm{\Sigma }}}_{0}$ is the initial density at each position of the disk and $\delta {\rm{\Sigma }}={\rm{\Sigma }}-{{\rm{\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\mbox{--}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 × 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 it also excites closer to the planet when we increase the planetary mass, as we discuss in the following section.

Zoom In Zoom Out Reset image size

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

$$\begin{eqnarray}&&\zeta \equiv \displaystyle \frac{{\rm{\nabla }}\times {\boldsymbol{v}}}{{\rm{\Sigma }}}.\end{eqnarray} \tag{ 7 }$$

For a barotropic fluid such 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}=({\rm{\nabla }}\times {{\boldsymbol{v}}}_{0})/{{\rm{\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 it as expected from linear theory. Linear waves from a low-mass planet (${M}_{p}\ll {M}_{\mathrm{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):

$$\begin{eqnarray}&&| {l}_{\mathrm{sh}}| \approx 0.93{\left(\displaystyle \frac{\gamma +1}{12/5}\displaystyle \frac{{M}_{p}}{{M}_{\mathrm{th}}}\right)}^{-2/5}h,\end{eqnarray} \tag{ 8 }$$

where γ is the adiabatic index and h is the scale height. With ${M}_{p}=0.1\,{M}_{\mathrm{th}}$ and ${(h/r)}_{p}=0.07$, ${l}_{\mathrm{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 Figures 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}_{\mathrm{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}_{\mathrm{th}}$ and ${(h/r)}_{p}=0.07$), it is unlikely that additional spiral arms form even beyond $r=2\,{r}_{p}$ (J. Bae & Z. Zhu 2017, in preparation).

3.2. Results with a 3 ${M}_{\mathrm{th}}$ Planet

We now describe results for a more massive planet, $3\,{M}_{\mathrm{th}}$ (corresponding to $1.03\,{M}_{\mathrm{Jup}}$ with ${M}_{* }=1\,{M}_{\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}_{\mathrm{th}}$ to $3\,{M}_{\mathrm{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}_{\mathrm{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 nonlinear effects shift the launching radius of the tertiary arm more than that of the secondary arm (J. Bae & Z. Zhu 2017, in preparation). Also, the arm-to-arm separation increases as the planetary mass increases, consistent with previous studies (e.g., Fung & Dong 2015; Zhu et al. 2015; Lee 2016). At $r=0.6\,{r}_{p}$, the secondary arm is offset by 112° from the primary arm as opposed to 50° 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.

Zoom In Zoom Out Reset image size

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 Figures 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.

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}_{\mathrm{th}}$, corresponding to $11.4\,{M}_{\oplus },0.10\,{M}_{\mathrm{Jup}},0.17\,{M}_{\mathrm{Jup}},0.34\,{M}_{\mathrm{Jup}},0.69\,{M}_{\mathrm{Jup}}$, and $1.03\,{M}_{\mathrm{Jup}}$ when ${M}_{* }=1\,{M}_{\odot }$ and ${(h/r)}_{p}=0.07$ are assumed.

We computed the gap depth ${\rm{\Delta }}{{\rm{\Sigma }}}_{\mathrm{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.

Zoom In Zoom Out Reset image size

The gap depth calculated following this method is plotted in Figures 3(b) and (c). Looking first at planets with ${M}_{p}\leqslant 1{M}_{\mathrm{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}_{\mathrm{th}}$) to a factor of ∼3 ($1\,{M}_{\mathrm{th}}$) in terms of ${\rm{\Delta }}{{\rm{\Sigma }}}_{\mathrm{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}_{\mathrm{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 Figures 3(a) and (c), the RWI is not able to get rid of the rings/gaps, but it only impedes further growth of them. 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 (J. Bae & Z. 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}_{\mathrm{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 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.

Zoom In Zoom Out Reset image size

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):

$$\begin{eqnarray}&&{\rm{\Sigma }}(r)={{\rm{\Sigma }}}_{c}{\left(\displaystyle \frac{r}{{r}_{c}}\right)}^{-\gamma }\exp \left[-{\left(\displaystyle \frac{r}{{r}_{c}}\right)}^{2-\gamma }\right],\end{eqnarray} \tag{ 9 }$$

where $\gamma =-0.23$ and ${r}_{c}=79\,\mathrm{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.⁵ 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}_{\mathrm{in}}=0.05\,{r}_{p}$ and ${r}_{\mathrm{out}}=2\,{r}_{p}$ and use 4096 × 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. Because 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, as the gas flow near the planetary orbit can be extremely complicated (e.g., Fung et al. 2015; Fung & Chiang 2016), thus the dust distribution in the corotation region that 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.

Zoom In Zoom Out Reset image size

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 because 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 ±20% of change in temperature results in a noticeable change in the gap locations: the secondary gap location changes by about ±4 au, which is about the linear resolution for the Band 6 ALMA continuum data, and the tertiary gap location changes by about ±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 (∼47 au) and B7 (∼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).

Observability of the secondary/tertiary gaps. In general, $\sim 10 \%$ of local density enhancement—corresponding approximately to ${\rm{\Delta }}{{\rm{\Sigma }}}_{\mathrm{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}_{\mathrm{th}}$ case have ${\rm{\Delta }}{{\rm{\Sigma }}}_{\mathrm{gap}}\geqslant 1.2$ after 1000 planetary orbits. This suggests that even for the planets with ${M}_{p}\lt 1{M}_{\mathrm{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 μ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 ∼2.4 and ∼5.0. As our $3\,{M}_{\mathrm{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 3D, it is possible that planet-driven spiral waves can become unstable to the spiral wave instability (SWI; Bae et al. 2016a, 2016b). On the other hand, the 3D simulations in Bae et al. (2016a) 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 (J. Bae & Z. 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 × 10⁻³ for the fiducial model presented in Section 3.2 ($3\,{M}_{\mathrm{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 measurements 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.

Zoom In Zoom Out Reset image size

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 ($\lt 1{M}_{\mathrm{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}_{\mathrm{th}}$ model. For the HD 163296 disk, presumably a more massive planet ($\gtrsim 1{M}_{\mathrm{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}_{\mathrm{th}}$ model.

Future work. In a forthcoming paper (J. Bae & Z. 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.

J.B. 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.