Applying a Particle-only Model to the HL Tau Disk

Recent high-resolution observations of a proplanetary disk around the young (∼1 Myr) T Tauri star HL Tauri by the Atacama Large Millimeter/submillimeter Array (ALMA) have revealed unprecedented detailed structures, which are considered likely to be the signatures of planets in the making. This image was taken as part of ALMA's science verification phase in 2014 October and was released a month later (see NRAO 2014). The disk was observed in dust continuum emission at 233 GHz (1.28 mm) using 25–30 antennas and a maximum baseline of 15.24 km as part of ALMA's Long Baseline Campaign, and achieved an angular resolution of 35 mas, equivalent to 5 au at HL Tau's distance of 130 pc. It reveals a series of concentric gaps that have become the subject of many studies, shedding light on the properties of the planets that are believed to be carving out these gaps, and providing a better understanding of the processes involved in the formation and evolution of planets and planetary systems. Table 1 lists orbital radii and the masses of the potentially hidden planets in five of the seven gaps that can be identified in the ALMA image, derived from past studies that used hydrodynamic and numerical simulations as well as analytic estimates which we briefly review here. We will compare the planetary masses derived in the literature to the results of our "particle-only" model in Section 3.

Table 1.  Estimated Orbital Radii and Masses of Possible Planets Forming HL Tau's Five Major Gaps from the Literature

${r}_{p}$ (au)					M $({M}_{{\rm{J}}})$					M⋆ (M⊙)	Method	References
${11.2}_{-0.1}^{+0.2}$	${29.7}_{-2.9}^{+2.9}$	${64.2}_{-0.3}^{+0.3}$			${0.40}_{-0.00}^{+0.02}$	${0.02}_{-0.02}^{+0.03}$	${0.21}_{-0.01}^{+0.02}$			0.55	Particle-only numerical sim. (Wisdom–Holman)	This work
13.1	33.0	68.6			0.35	0.17	0.26			0.55	HD & radiative transfer	Jin et al. (2016)
11.8	32.3	82			0.77	0.11	0.28			0.55	2D HD sim., Gap width measurements (Equation (2))	Kanagawa et al. (2016)
${13.5}_{-0.4}^{+0.4}$	${32.4}_{-0.4}^{+0.6}$	${65.2}_{-0.9}^{+1.3}$	${77.2}_{-0.7}^{+0.8}$		>0.85	>0.61	>0.62	>0.51		0.55	Gap depth measurements (Equation (1))	Akiyama et al. (2016)
${13.6}_{-0.2}^{+0.2}$	${33.3}_{-0.2}^{+0.2}$	${71.2}_{-0.5}^{+0.5}$	${93.0}_{-0.9}^{+0.9}$				≲0.30 No MMR ≳0.72 MMR	≲0.30 No MMR ≳0.72 MMR		0.55	N body sim., REBOUND package (Rein & Liu 2012)	Tamayo et al. (2015)
${13.6}_{-0.2}^{+0.2}$	${33.3}_{-0.2}^{+0.2}$	${65.1}_{-0.6}^{+0.6}$	${77.3}_{-0.4}^{+0.4}$	${93.0}_{-0.9}^{+0.9}$			≲0.11 No MMR ≲0.30 MMR	≲0.11 No MMR ≲0.30 MMR	≲0.11 No MMR ≲0.30 MMR
${11.3}_{-0.1}^{+0.2}$	${29.4}_{-3.6}^{+2.6}$	${63.7}_{-0.4}^{+0.5}$			${0.81}_{-0.01}^{+0.02}$	${0.04}_{-0.04}^{+0.04}$	${0.37}_{-0.03}^{+0.01}$			1.3	Particle-only numerical sim. (Wisdom–Holman)	This work
13.2	32.3	68.8			0.2	0.27	0.55			1.3	3D dust+gas SPH	Dipierro et al. (2015)
11.8	32.3	82			1.4	0.2	0.5			1.0	2D HD sim., Gap width measurements (Equation (2))	Kanagawa et al. (2016)
${13.2}_{-0.2}^{+0.2}$	${32.3}_{-0.1}^{+0.1}$	${64.2}_{-0.1}^{+0.1}$	${73.7}_{-0.1}^{+0.1}$	91	$\gt {0.02}_{-0.01}^{+0.01}$	$\gt {0.07}_{-0.01}^{+0.01}$	$\gt {0.03}_{-0.01}^{+0.00}$	$\gt {0.08}_{-0.05}^{+0.03}$	$\gt {0.11}_{-0.06}^{+0.03}$	1.7	Surface density measurements	Pinte et al. (2016)

Note. If the two gaps at D5 and D6 are formed by a single planet at 71.2 au, Tamayo et al. (2015) estimated the masses of the two outermost planets to be ≲0.30 M_J if the two planets at D5+D6 and D7 are not in MMR and ≳0.72 M_J if they are. In the four-planet scenario, Pinte et al. (2016) determined the mass of a single planet at 69.0 au forming the gap at D5+D6 to be at least ${0.44}_{-0.09}^{+0.05}\,{M}_{{\rm{J}}}$, though they used a much larger stellar mass to derive those planetary masses. Kanagawa et al. (2016) used two different stellar masses in their measurements of planet masses based on gap widths, both of which are listed here. Their mass measurement for the planet orbiting around a 1 M_⊙ star at ∼30 au is consistent with their earlier paper (Kanagawa et al. 2015) where they had determined the mass from gap depth to be at least 0.3 M_J using the same stellar mass. The table also lists the results of our work, determined for two stellar masses of 0.55 and 1.3 M_⊙, explained in subsequent sections.

Download table as:  ASCIITypeset image

ALMA Partnership et al. (2015) identified seven pairs of distinct dark and bright rings in the ALMA image of the HL Tau disk which they labeled D1...D7 and B1...B7 (more on this in Section 4). They approximated the radial distance of the center of each ring by making a cross-cut along the disk's major axis and found the dark rings to be at 13.2 ± 0.2, 32.3 ± 0.1, ∼42, ∼50, 64.2 ± 0.1, 73.7 ± 0.1, and ∼91.0 au, placing the first four dark rings in a chain of MMRs, specifically 1:4:6:8. Pinte et al. (2016) measured the missing dust mass in each of the seven gaps by integrating the dust surface density of each gap and comparing it to its surrounding bright rings. They argued that these provide the mass of the rocky core of the possibly embedded planet in each dark ring.

Other authors have also attempted to constrain the masses of the planets that are believed to be shepherding the HL Tau gaps. Based on the depth of the gap seen ∼30 au from a central star, Kanagawa et al. (2015) estimated the mass of its embedded planet. They did so by using the relationship between the depth of a gap formed by a planet in its feeding zone in a protoplanetary disk and the mass of the planet as well as the disk's viscosity and scale height (Duffell & MacFadyen 2013; Fung et al. 2014; Kanagawa et al. 2015, see), given by (Kanagawa et al. 2015):

$$\begin{eqnarray}&&\displaystyle \frac{{M}_{p}}{{M}_{\star }}=5\times {10}^{-4}{\left(\displaystyle \frac{1}{{{\rm{\Sigma }}}_{p}/{{\rm{\Sigma }}}_{0}}-1\right)}^{1/2}{\left(\displaystyle \frac{{h}_{p}}{0.1}\right)}^{5/2}{\left(\displaystyle \frac{{\alpha }_{\mathrm{ss}}}{{10}^{-3}}\right)}^{1/2},\end{eqnarray} \tag{ 1 }$$

where M_p is the planet's mass in stellar mass units M_*, Σ_p/Σ₀ is the gap depth which is the ratio of the surface density of the planet-induced gap to that of the unperturbed disk, h_p is the disk's aspect ratio at the planet's orbital radius (h/r, with h being the scale height), and α_ss is the Shakura–Sunyaev kinematic viscosity parameter (Shakura & Sunyaev 1973).

Adopting a stellar mass of 1.0 M_⊙, a viscosity parameter of 10⁻³, and estimating the gap depth and the disk's aspect ratio to be ∼1/3 and ∼0.07 respectively, Kanagawa et al. (2015) were able to determine that the mass of the planet at 30 au is at least 0.3 M_J.

Also, using the gap depth (Equation (1)) and a method based on angular momentum transfer analysis in gas disks, Akiyama et al. (2016) estimated the masses of the planets in the HL Tau system to be comparable to or less than 1 M_J.

Estimating the dust mass deficits in the gaps as done by Pinte et al. (2016), Kanagawa et al. (2015), and Akiyama et al. (2016) provides a lower limit for the planet masses since an accurate measurement of the gap depth requires high signal-to-noise ratio data, otherwise the gaps cannot be fully resolved and seem to be partially "filled in" (Pinte et al. 2016). For this reason, in a follow-up paper, instead of using the gap depth to measure the masses of the planets, Kanagawa et al. (2016) derived an empirical relationship between the width of a planet-induced gap and planet mass, disk aspect ratio, and viscosity. Using two-dimensional hydrodynamic simulations, and assuming that the dust particles are strongly tied to the gas (i.e., dust filtration is not a major concern), they determined this relationship to be:

$$\begin{eqnarray}&&\displaystyle \frac{{M}_{p}}{{M}_{\star }}=2.1\times {10}^{-3}{\left(\displaystyle \frac{{{\rm{\Delta }}}_{\mathrm{gap}}}{{r}_{p}}\right)}^{2}{\left(\displaystyle \frac{{h}_{p}}{0.05{r}_{p}}\right)}^{3/2}{\left(\displaystyle \frac{{\alpha }_{\mathrm{ss}}}{{10}^{-3}}\right)}^{1/2},\end{eqnarray} \tag{ 2 }$$

with r_p and Δ_gap being the orbital radius of the planet and the width of the gap it creates in the disk, respectively, where the gap edges are defined by regions where the surface density drops to less than half the unperturbed surface density.

Using Equation (2) to determine planetary masses probably results in more accurate estimates than simply using the size of each planet's Hill radius, which tends to predict rather large planetary masses.³ For the planets in the HL Tau disk, Akiyama et al. (2016) measured the gap widths to be 5.0, 4.1, 6.2, and 4.5 au, at r_p = ∼13.5, ∼32.4, ∼65.2, and ∼77.2 au and used the size of each planet's Hill sphere to calculate its mass. This resulted in planetary masses of ${88.8}_{-5}^{+5}$, ${3.6}_{-0.6}^{+0.7}$, ${1.5}_{-0.5}^{+0.5}$, and ${0.3}_{-0.1}^{+0.1}$ M_J, much larger than other mass estimates for the HL Tau planets, especially the innermost planet whose Hill radius suggests a stellar-mass body. Although the possibility of low-mass stellar companions in this system is not ruled out, high-sensitivity direct imaging in the mid-infrared by Testi et al. (2015) did not reveal any point sources in the HL Tau disk. Although their observations were more focused around the gaps at ∼70 au (for which the contrast level reached was ∼7.5 mag.), their search for point sources in the HL Tau disk was not exclusive to the outer disk. Nevertheless, to examine the possibility of stellar/substellar companions in the HL Tau disk, further studies are needed to determine the stability of the system under such conditions. Therefore, we exclude mass measurements from planetary Hill radii when we later compare our results to those of others (see Section 3.1).

It is worth noting here that dust filtration by the planet's pressure bump as well as dust migration under radiation and drag forces can also cause gaps to be filled in temporarily and yield inaccurate measurements of the gap width and depth. Thus the mass of a gap-opening planet derived from the width and depth of the gap must be taken with caution, particularly in cases where the disk is massive and hence there is high rate of collisional fragmentation down to grain sizes that are affected by radiation and drag forces. However, such effects are more important when planetary masses are estimated from the depths of the gaps than their widths (see Dong et al. 2015).

Using hydrodynamic simulations and radiative transfer models, Jin et al. (2016) attempted to match the width and depth of the three prominent gaps in the HL Tau disk, located at at 13.1, 33.0, and 68.6 au and constrained the masses of the planets that are believed to be in those gaps to be 0.35, 0.17, and 0.26 M_J, respectively, assuming no planet migration through the disk. The model assumes a disk mass of ∼7.35 × 10⁻² M_⊙ and the same α_ss parameter as Kanagawa et al. (2015) while the dust-to-gas ratio is taken to be 1%. Furthermore, the authors also tried to match the eccentricities of the gaps where they placed the three planets and found them to be 0.246, 0.274, and 0.277, respectively. On the other hand, smoothed particle hydrodynamic (SPH) models by Dipierro et al. (2015) constrained the masses of the planets embedded in the HL Tau disk to be 0.2, 0.27, and 0.55 M_J with planets at 13.2, 32.3, and 68.8 au.

Gas and dust interact differently with planets. Numerical simulations by Jin et al. (2016) showed that the three gaps formed by tidal interaction with the embedded planets in the HL Tau disk are shallower in gas distribution and deeper in dust, though both have similar morphologies (see their Figure 1). The difference in the gap's gas and dust surface density arises from the fact that submillimeter dust is pushed toward the edges of the gap as it starts to open since gas drag tends to accumulate dust particles in high-pressure regions as suggested by the enhanced dust emission near gap edges (see Haghighipour & Boss 2003; Fouchet et al. 2007; Maddison et al. 2007).

Most authors place three planets in the HL Tau disk; however, the possibility of additional planets in this disk has also been discussed in the literature. For instance Tamayo et al. (2015) considered the possibility of up to five planets in the HL Tau disk at nominal radii of 13.6, 33.3, 65.1, 77.3, and 93.0 au. This places the outer three planets nearly in a chain of 4:3 MMR. The authors determined the masses of the five planets under two different scenarios: if the planets are not in MMR, they found a maximum mass of ∼2 Neptune masses for the outer three bodies; however, if the outer three planets are in resonance, as suggested by the locations of the gaps, they can grow to larger masses via resonant capture as they migrate through the disk during which their masses can reach at least that of Saturn. The masses of the two inner planets were not well constrained in this study since these planets are dynamically decoupled from the other three.

Planets forming in a multi-planet system can grow to where the system becomes unstable simply because of the growth in the sizes of the planets' Hill spheres. Planets whose orbits around the star are separated by less than several of their mutual Hill spheres are unstable: this stability criterion is defined by Gladman (1993) who suggested planets that are separated by less than 3.46 r_H destabilize on a timescale that is roughly their conjunction period. On the other hand, Tamayo et al. (2015) showed through numerical simulations that planets can still survive well beyond the above stability criterion if they capture in MMR at low masses and grow together. This is because resonance mitigates the effect of close encounters. For the system to be Hill stable, the maximum masses for the three outermost planets in the HL Tau disk were found by Tamayo et al. (2015) using Equation (3), taking the stellar mass to be M_* = 0.55 M_⊙:

$$\begin{eqnarray}&&M\lesssim {M}_{\mathrm{crit}}=8\,\left(\displaystyle \frac{{M}_{\star }}{{M}_{\odot }}\right){\left(\displaystyle \frac{{\rm{\Delta }}a/a}{0.1}\right)}^{3}\,{M}_{\oplus },\end{eqnarray} \tag{ 3 }$$

where Δa is the planet separation and M_crit is the maximum mass to ensure stability.

Therefore, according to numerical simulations of the HL Tau disk by Tamayo et al. (2015), if the outer three planets are not in MMR they become unstable at conjunction timescale once they exceed the mass threshold beyond which their separation becomes less than ∼3.5 r_H. However, if they are captured at resonances while they migrate through the disk, they can grow well past the above limit until they become so massive (∼40% beyond mass of Saturn or 0.44 M_J) that their mutual gravitational perturbation at conjunctions brings them out of resonance at which point swift instability ensues (Tamayo et al. 2015).

Moreover, their numerical simulation suggested that the system would be substantially more stable if not all the gaps were made by planets, particularly the more closely spaced gaps at 64.2 and 73.7 au, or D5 and D6 according to ALMA Partnership et al. (2015). (Note that these two gaps are at 65.1 and 77.3 au in Tamayo et al. 2015). They suggested that these two gaps might not be made by two different planets; there might instead be a single planet at 71.2 au that has shaped a horseshoe-like gap in the disk of HL Tau. If four planets are considered instead of five in the HL Tau system, their numerical simulations put a final mass limit of at least 230 M_⊕ for the outer two planets if they are in MMR while they can reach ∼1 Saturn mass if they are not.

Other authors have also suggested that the double gap at D5 and D6 in the HL Tau disk could be made by a single planet exciting Lindblad torques with the bump in the middle being co-orbital horseshoe material where the planet is possibly hiding. Using 2D gas+dust hydrodynamic simulations combined with radiative transfer modeling, Dong et al. (2017) showed how super-Earths placed in a low-viscosity disk can produce characteristic double gaps in mm-dust distribution and argued that the D5+D6 gaps in the HL Tau disk could be carved by a single planet. Besides the double gap on either side of the planet's orbit, their simulations also suggested that additional gaps could arise in the disk for a single planet, depending on disk and planet parameters. In a more recent paper by Bae et al. (2017), using 2D hydrodynamic simulations of the gas component of the HL Tau disk, the authors were able to reproduce not only D5 and D6 with a single planet at 68.8 au, but also noticed that the same planet can reproduce the D1 and D2 gaps at 13.2 and 32.3 au, though their model did not reproduce the finer structures such as D3, D4, and D7. They also argued that the mass of the planet can be constrained from the positions of multiple gaps, provided that the disk temperature profile can be accurately measured.

Table 1 summarizes the masses and locations of the possible planets in the HL Tau system obtained by the studies mentioned above. It is important to note that the masses derived for the HL Tau planets in the literature depend on the mass of the central star, which is not well constrained. Estimates of HL Tau's stellar mass range from 0.55 M_⊙ (e.g., Tamayo et al. 2015) to 1.7 M_⊙ (e.g., Pinte et al. 2016). Therefore, in order to be able to compare our results to those of others, we only focus on two previously used values of 1.3 and 0.55 M_⊙. Nevertheless, even for the same stellar mass, Table 1 shows that, despite many attempts to constrain planetary parameters in the HL Tau disk, much work is still needed to determine the number and parameters of its potentially embedded planets. This is the primary motivation of this work: can we reproduce the key features of the HL Tau disk using the computationally inexpensive model of a particle-only disk to address whether some parameters of its planets, specifically their mass and orbital radii, can be determined without the need for sophisticated models which are, nevertheless, required to fully describe gas-rich disks? In Table 1, we also list our results for comparison with others but will explain how we arrived at these values in the next two sections.

The observed profile of the HL Tau disk used here is extracted by the following method. First, we obtained the FITS image of the HL Tau disk available publicly at the ALMA website and observed in dust continuum emission in band 7 (the highest resolution). We made a cross-cut across the disk's major axis to extract HL Tau's radial brightness profile in Jy/beam per radial distance from the star. The extracted profile is 186 pixels long over a physical distance of 115 au. However, the resolution of the image is only 35 mas or ∼5 au at 130 pc (see NRAO 2014) and so we assess that we really only have 115/5 ≈ 23 bins for the purposes of determining our degrees of freedom (see Section 2.4) and 186/23 ≈ 8 pixels per bin.

Our simulations are performed with a symplectic integrator based on the Wisdom–Holman algorithm (Wisdom & Holman 1991). A fixed timestep of 150 days is used for all simulations. Only point particles are simulated, without any gas drag, radiation pressure, or Poynting–Robertson drag. These effects are likely to be important in sculpting the HL Tau disk but our purpose here is to determine what, if any, of the planetary parameters can be recovered by the simplest possible model.

Simulations are run for 10,000 yr (∼1000 inner orbits) and recorded at 100 yr intervals. Three planets and 1000 particles are placed within the disk on circular orbits around a 1.3 or 0.55 solar-mass central star. Particles are removed if they reach a distance less than ∼500 solar radii or greater than 220 au. The planets are placed nominally at 11.7, 29.1, and 64.5 au based on the locations of the gaps in the HL Tau disk, but the planets' locations will be varied as part of the fitting process, described in Section 2.3.

Simulated disk profiles are created from the last five snapshots of the disk. The use of several snapshots increases our signal-to-noise without the computational expense associated with simulating additional particles, though it assumes that the disk is in a quasi-steady state. Examination of the disk during the final stages confirms that indeed the disk structures are well-established.

For plotting purposes, the simulation data are extracted into a histogram with 186 bins to match the observations. The bins are weighted by the blackbody emission of their particles assuming a dust albedo of 0.5 and emissivity of 1.0 at mm wavelengths to calculate the equilibrium temperature of the disk particles. The stellar luminosity and effective temperature are also uncertain but are taken to be 8.3 L_⊙ and 4000 K, respectively (Ruge et al. 2016). For calculation of the χ² of our fits, the data are box-car smoothed down to the effective resolution of the observations (8.3 bin box-car). On the basis of the χ² value, new parameter values are chosen and a new simulation is initialized. The whole process is iterated until convergence is achieved.

Best-fit parameters are established on the basis of the χ² between the observational profile and a simulated profile normalized to the first bin in the observed profile. This normalization reduces our degrees of freedom by one. Minimization of the χ² parameter is accomplished using Interactive Data Language (IDL) and the Amoeba package, which is a multi-dimensional derivative-free optimization algorithm based on the downhill simplex method of Nelder and Mead (1964). Typical Amoeba runs require 900–1000 simulations and a total of 10 hr to complete on a single CPU. Amoeba requires the tolerance to be at least equal to the machine's double precision, so we set the tolerance to 10⁻¹². This is the decrease in the fractional value of the χ² in the terminating step.

Chi-squared minimization using the Amoeba algorithm does not require calculating derivatives. Furthermore, each iteration only takes one or two function evaluations and therefore Amoeba converges faster than some other minimization routines such as nonlinear least-squares fitting using the Levenberg–Marquardt algorithm (Marquardt 1944; Levenberg 1944) which takes several calculations per iteration. Amoeba is also more robust for problems with stochastic components such as what we are dealing with here (e.g., the particle positions are chosen randomly for each simulation, which introduces some statistical noise to the radial profiles), and we chose the Amoeba algorithm for these reasons.

A downside to using Amoeba is that it can get to a point where the changes in the parameter values become insignificant before a minimum is reached. Thus it is generally recommended to restart Amoeba from the point where it claims to have found a minimum (see Press et al. 1992) and this is what we do 10 times until the routine converges again. Our procedure was to first perform initial minimization runs using parameter values chosen arbitrarily, except for the orbital radii of the three planets (these were estimated from the locations of the major gaps in the HL Tau disk) and each parameter was allowed to vary by ±50% by the minimization routine. From the lowest χ² obtained from these initial runs (our "initial solution"), in order to ensure as much as possible that the minimum χ² achieved is the global minimum, we performed 10 additional minimization runs where we changed the initial conditions such that each parameter fell randomly within 10% of that obtained from the initial solution. At the end, we recorded the parameters that produced the lowest χ² from the 10 + 1 Amoeba runs. Our restarting process helps avoid terminating at a spurious local minimum, but we cannot exclude the possibility of a true global minimum that might exist far away from our final result in parameter space.

For our simulations here, we fit 10 parameters of the planets and disk (in our model with the broken power-law but seven when we use a single power-law for disk density distribution, see Section 3). We assume that there are three planets on circular orbits. In addition to the masses and orbital radii of these three planets, we also fit a power law to the disk surface density. The surface density of circumstellar disks is generally taken to have a profile of the form Σ ∝ R^−α with the power-law index, α, between 0 and 1 depending on the mass of the protoplanetary disk (Andrews & Williams 2007). (Note that the power-law index derived from Minimum Mass Solar Nebula is 1.5; Weidenschilling 1977). However, the use of a single power law does not well reproduce the radial profile of HL Tau's flux density. A much lower χ² value is obtained by selecting a different power-law index beyond the location of the outermost planet (see Section 3). Yen et al. (2016) also used a broken power law in their measurements of gap widths and depths in the HL Tau disk where the slopes of the dust distribution based on the column density of HCO⁺, assuming that gas and dust are well coupled thermally, were found to be 0.5 ± 0.2 at ∼20 au and 0.9 ± 0.3 at ∼60 au, suggesting a steep decline in dust continuum emission beyond where the outer major gap lies. Jin et al. (2016) proposed that the deficit in dust in the outer part of the HL Tau disk is due to the inward drift of dust caused by gas drag and the absence of a source to supply the dust at large radii (also see Birnstiel & Andrews 2014). In fact, disks are found to have exponentially tapered edges and an exponential decrease in dust surface density has also been observed for a number of other circumstellar disks (see for instance, McCaughrean & O'dell 1996) with power-law indices beyond the above-mentioned range, suggesting that the pure power-law relation (i.e., Σ ∝ r^−α) does not accurately represent a disk's intensity profile and must be replaced by an exponentially truncated density distribution with Σ ∝ r^−γ, where γ is the exponent in the viscosity dependence on distance from the star (e.g., Hartmann et al. 1998). For simplicity (that is, to avoid adding additional parameters to our fit), we assume a standard power-law slope without an exponential term. Therefore, we argue that our fit to the radial profile of the outer disk would be improved if we adopted the above surface density profile and incorporated dust re-generation and gas drag in our model, which we leave to future work.

Uncertainties in the fitted parameters are estimated based on the χ² values. The number of degrees of freedom, ν, will be the effective number of bins (23, see Section 2.1) minus one for the normalization discussed in Section 2.3, and minus one for each free parameter. We have 10 free parameters, giving us a total of 12 degrees of freedom.

The uncertainties to be at the locations in phase-space can be approximated as where the χ² value is increased over its minimum value by an amount Δχ² dependent on ν and the stringency of the uncertainty bounds desired. Here we choose a p = 0.95 (nominally 2σ) confidence region, which means that our uncertainties correspond to the locations for which (Press et al. 1992)

$$\begin{eqnarray}&&Q\left(\displaystyle \frac{\nu }{2},\displaystyle \frac{{\rm{\Delta }}{\chi }^{2}}{2}\right)=1-p,\end{eqnarray} \tag{ 4 }$$

where Q is the incomplete gamma function, and Δχ² gives the increase in χ² corresponding to our uncertainty.

Note that we compute our uncertainties from χ² values with all the parameter values except the one in question held constant. This implicitly assumes that the parameters are uncorrelated, which we assume here for reasons of simplicity and practicality. Our χ² is derived by a process with inherent stochasticity (i.e., the initial conditions of particles within the disk have a random component), thus we have too many free parameters and too noisy a system to determine the covariance between them all effectively. This will be more apparent when the uncertainty results are discussed in Section 3.

A comparison between planetary parameters (masses and orbital radii) that we obtained for the M_* = 0.55 M_⊙ case using our simple model of the HL Tau disk shows that, except for the mass of the middle planet which is underestimated by our model, these parameters are comparable to what some authors found using models of higher complexity. In particular, our mass measurements for planets 1 and 3 are similar to those found by Jin et al. (2016) who used hydrodynamic gas+dust simulations coupled with 3D radiative transfer calculations. Planet 3's mass is also comparable to what is suggested by Kanagawa et al. (2016) through their hydrodynamic simulations. However, those authors note that if the gap they measured was narrower by ∼2 au, the mass of the innermost planet would be the same as what Jin et al. (2016) found (which is similar to our result).

Planetary masses derived using Equation (2) depend on r_p since gap widths are scaled by the location of their centers. Kanagawa et al. (2016) identified the three prominent gaps in the HL Tau disk from the radial profile of the optical depth in band 6 which is offset from the radial profile of the dust brightness temperatures in continuum emission at the locations of the first (D1) and third (D5+D6) planets (see their Figures 1(a) and (b)). Compared to their plot of the temperature profile of dust, D1 is shifted inward while D5+D6 is shifted outward in optical depth. This means that if r_p was determined from dust temperature, the mass of the planet in D1 would be less than what they report while that of the planet in D5+D6 would be larger. Also they found r_p by taking (r_in + r_out)/2 where (r_in and r_out) are the inner and outer edges of each gap. This assumes that the gaps are symmetric, but in fact they are slightly asymmetric. This could affect r_p and therefore their calculation of planetary mass from the gap width. Furthermore, Kanagawa et al. (2016) also noted that their mass estimates depend strongly on the disk scale height (and hence temperature) as well as dust opacity spectral index, both of which need to be well constrained for the planet mass to be determined accurately. Determining the viscosity parameter is also important when using the formula given by Kanagawa et al. (2016) to measure planetary masses, although the dependence is not as strong since ${M}_{p}\propto {\alpha }_{\mathrm{ss}}^{1/2}$.

Compared to our results for the planet masses, the masses derived by Akiyama et al. (2016) are overestimated. This could be due to the fact that they used gap depths (i.e., Equation (1)) to find planet masses and, as discussed in Section 1.1, high signal-to-noise ratio data are required to measure the emission at the bottom of the gap and determine Σ_p/Σ₀ (see Kanagawa et al. 2016). Note that both Equations (1) and (2) apply to gap depth and width in gas emission but assume that they are similar for dust gaps, which is true if gas and dust are well mixed and dust filtration is not strong. However, studies show that even if dust filtration is weak (which is the case for relatively massive disks such as the HL Tau: M_disk = 0.07–0.17 M_⊙; see Dong et al. 2015; Jin et al. 2016), gas gaps are shallower than dust gaps though the widths remain comparable in gas versus dust. This is because filtration affects gap depths more than their widths (see Dong et al. 2015; Yen et al. 2016). Thus planetary masses are more accurately measured using gap widths (i.e., Equation (2)) than gap depths (i.e., Equation (1)). According to Equation (1), shallower gaps result in overestimating planetary mass, which is likely why planetary masses found by Akiyama et al. (2016) are larger than those found by ourselves and Jin et al. (2016).

Tamayo et al. (2015) did not constrain the masses of the two planets in D1 and D2 since they are dynamically decoupled from the other planets they placed in their simulations. They did, however, determine the limit for the mass of the planet in D5+D6 by letting it grow together with a fourth planet which they placed in D7 (∼90 au) under two scenarios: M_p ≲ 0.30 M_J if the two planets are not in MMR and M_p ≳ 0.72 M_J if they are. We did not put any planet in D7 and leave the investigation of the possibility of additional planets in the HL Tau disk to a future paper, so we shall not comment on how our results compare to theirs other than naively mentioning that in both their four- and five-planet simulations, the masses derived are sub-Jovian (except perhaps where they suggest a lower limit of 0.72 M_J for the outer two planets in MMR in the four-planet case), which is consistent with the other studies mentioned here, ours included. However, we acknowledge that placing additional planets in the system may in fact affect our results, which we defer to future work.

The mass of the central star in the HL Tau system is not well known. Estimates based on the Keplerian velocity of gas (e.g., Sargent & Beckwith 1991; Pinte et al. 2016) or protostellar evolutionary tracks (e.g., Beckwith et al. 1990; Güdel et al. 2008) suggest a star of mass 0.55–1.7 M_⊙. Therefore, we tried our simulations with a higher stellar mass to see how well our results match those of Dipierro et al. (2015) for a stellar mass of 1.3 M_⊙. Again, the mass of our second planet is significantly lower than theirs while they found the mass of the first planet to be much less than what we did, though the mass of our third planet is comparable to theirs. Since planet masses scale with the stellar mass, which is also true in our simulations when our results for the two different stellar masses are compared against each other, it is not clear why Dipierro et al. (2015) found the mass of the first planet to be half that of, for instance, Jin et al. (2016) even though the stellar mass is more than doubled.

Our best-fit parameters for the power-law indices of the disk's surface density profile are ${\alpha }_{1}={0.26}_{-0.02}^{+0.02}$ and ${\alpha }_{2}={4.96}_{-0.05}^{+0.07}$ where the break occurs at ${r}_{b}={69.70}_{-0.12}^{+0.33}$. We note that these two values are very different from each other and from what Yen et al. (2016) report (see Section 2.3), partly owing to the fact that we introduced a sudden increase in the disk's surface density by almost a factor of 3 (i.e., f = 2.69) where we broke the intensity profile of the disk into the two segments. Furthermore, the Yen et al. (2016) values are derived from the the column density of HCO⁺ gas (though they assumed that dust is tightly coupled to the gas) whereas our model only includes solid dust particles. Nevertheless, our surface density slope in the outer disk is close to the value obtained by Pinte et al. (2016) derived from the missing dust mass in each gap. According to their model, the surface density profile of the HL Tau disk has a slope of −3.5 out to about 75 au but falls off faster in the outer part of the disk. They find the power-law slope in the outer disk to be −4.5, which is similar to what we obtained from our model. They attributed the change in the surface density of dust to two possible reasons: lack of efficient grain growth in the outer disk or the removal of a significant fraction of mm-sized grains from the outer disk via radial migration of dust.

Using a larger stellar mass in our simulations (i.e., 1.3 instead of 0.55 M_⊙), we find a similar value for α₂ ($={4.68}_{-0.05}^{+0.06}$) but α₁ is reduced by a factor of 2 (${\alpha }_{1}={0.18}_{-0.01}^{+0.01}$). In this case, when using a single power law for the disk's surface density, we find α to be ∼0.20. We therefore conclude that we would need a more complicated model for the disk surface density to better estimate the power-law indices in the two segments.

Nevertheless, our results explain the three prominent gaps at D1 and D2, and the double gap at D5+D6 and the planetary masses found are similar to the results of others, especially when the stellar mass is 0.55 M_⊙, while also reproducing some finer gaps, particularly at D3 and D7. In the next section we explain the possible nature of the two narrower gaps seen both in the ALMA image of the HL Tau disk and in our simulations. Therefore, our model is successful in reproducing the observed intensity profile of the HL Tau disk without the need to include certain elements that are necessary to fully study a gas disk.

The orbital radii of the planets found by our fitting procedure represent the locations of the three major gaps in the HL Tau image to within uncertainties (where the two gaps made by the outermost planet are considered to be a double gap separated by particles in 1:1 MMR with a planet at ∼64 au). A closer look at the observed intensity profile of the HL Tau disk reveals a few other narrower gaps which have motivated some authors to include more planets in their modeling of the HL Tau disk. However, our earlier studies, Tabeshian & Wiegert (2016, 2017), have shown that not all disk gaps need to contain planetary bodies and that some gaps can, in fact, be made via MMR with a planet that is located outside the gaps and can be used to learn about the hidden planets. We note that our model of the HL Tau disk, which only includes three planets, is able to reproduce some of those narrower gaps as well. In fact, given the locations of the second and the third planets, we argue that the gaps seen at ∼38 and ∼84 au, roughly corresponding to the locations of the D3 and D7 dark gaps in ALMA Partnership et al. (2015), are made by exterior 3:2 MMR with those two planets, respectively. These are shown by the vertical dotted lines in Figures 2 and 3. Furthermore, the locations of the first and second planets place them in a 4:1 MMR with each other.

HL Tau is considered a relatively massive disk in which rapid pericenter precession rates alter the location of resonances by an amount roughly given by the ratio of the disk mass to star mass (Tamayo et al. 2015). Taking M_disk to be = 0.13 M_⊙ (Kwon et al. 2011), this means that the disk is 25% the mass of the star which, according to Equation (10) of Tamayo et al. (2015), causes the location of the 3:2 resonance to move by <10%. However, they point out that, due to the uncertainty in calculating the precession rate, the exact locations of resonances in massive disks are uncertain to within ∼M_disk/M_⋆ as well.

In order to visually compare the result of our simulation with ALMA's image of the HL Tau disk, we make a simulated image using the Common Astronomy Software Applications (CASA) for simulating ALMA observations (McMullin et al. 2007) based on the disk produced with our best-fit parameters. To make the CASA simulated image, we assume that our disk is placed at the HL Tau distance of 130 pc and therefore has the same radial size on the sky. We also assume that the particles are perfect blackbodies at local thermal equilibrium and take the disk's total flux to be 700 mJy at 1.3 mm (Kwon et al. 2011). Stellar radius and effective temperature are 6.0 R_⊙ and 4000 K, respectively (Ruge et al. 2016). We set the image resolution at 35 mas or ∼5 au to match that of ALMA's observation of the HL Tau disk and use all the 50 available antennas in the 12 m array. We assume that the disk is observed for a total of 4 hr and set the integration time to 10 s per pointing. The R.A. and decl. of the center of the image are α = 04^h31^m3845 and δ = 18°13'590, J2000 (Tamayo et al. 2015). Beam deconvolution is done using CASA's CLEAN algorithm. The result is shown in Figure 5. We adopt the same nomenclature used by ALMA Partnership et al. (2015) for the dark gaps that we see in our simulations, except that we take the two gaps around the outermost planet to be the same with the planet in the middle.

Zoom In Zoom Out Reset image size

It must be noted that we are not claiming that the properties of gas-rich disks can be fully determined from simple models that do not incorporate gas and radiation forces. However, based on the ability of our simulations to reproduce the intensity profile of the HL Tau disk, we argue that reasonable matches with observations can be achieved with relatively simple particle-only models of this intrinsically much more complicated gas disk. Therefore, at least as far as understanding the dynamics of the system is involved, we make the case that simple models could be used to extract useful information about the number and properties of possible planets embedded in gas-rich disks which could be used in future, more thorough analyses of these disks.