MODELING DUST EMISSION OF HL TAU DISK BASED ON PLANET–DISK INTERACTIONS

We initialize our hydrodynamic simulation using the results from Kwon et al. (2011). The disk extends from 2.4 to 156 AU, and is simulated using a polar grid containing 3072 and 768 cells in the radial and azimuthal direction, respectively. We adopt a disk surface density described by

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{g}(r)={{\rm{\Sigma }}}_{0}{\left(\displaystyle \frac{r}{{r}_{0}}\right)}^{0.23}\mathrm{exp}[-{(r/{r}_{c})}^{2.23}],\end{eqnarray} \tag{ 1 }$$

where ${r}_{c}=79\;{\rm{AU}}$ is a characteristic radius (e.g., Andrews et al. 2009) and ${r}_{0}=10\;{\rm{AU}}$. The mass of the star is ${M}_{*}=0.55{M}_{\odot }$. Although Kwon et al. (2011) finds a disk mass of $0.135{M}_{\odot }$, we find that the disk with this mass quickly becomes gravitationally unstable under small perturbations. Hence we reduce the disk mass to be ${M}_{{\rm{disk}}}\approx 7.35\times {10}^{-2}$M_⊙, which leads to ${{\rm{\Sigma }}}_{0}=4.83\times {10}^{-4}$ $[{M}_{*}/{(10{\rm{AU}})}^{2}]\approx 23.61$ g cm⁻². The whole disk is assumed to have a constant Shakura–Sunyaev viscosity parameter $\alpha ={10}^{-3}$. The locally isothermal sound speed is chosen as

$$\begin{eqnarray}&&\left(\displaystyle \frac{{c}_{s}}{{v}_{\phi }}\right)(r)=0.05{\left(\displaystyle \frac{r}{{r}_{0}}\right)}^{0.285}\;\;.\end{eqnarray} \tag{ 2 }$$

We assume that the initial dust surface density ${{\rm{\Sigma }}}_{d}(r)$ follows ${{\rm{\Sigma }}}_{g}(r)$, with an initial dust-to-gas ratio of 1%. We have performed gas+dust two-fluid simulations using different dust particle sizes ranging from μm to millimeter, though only one dust size is used in a single hydro simulation. Such gas+dust distributions are used to calculate the dust temperature as well as the dust emissions, which are ultimately compared with ALMA observations. For comparison with the ALMA images, we have used the 0.15 millimeter size dust particles from gas+dust two-fluid simulations (see below for dust models). The initial Stokes number of 0.15 millimeter size dust particles is $\lt 0.01$ within ∼100 AU. For such a small Stokes number, we adopt short friction time approximation (Johansen & Klahr 2005) to circumvent problems with small time steps.

The disk self-gravity and indirect acceleration term from the central star are always included. An extended disk that includes an inner disk of $r\in [0.01,0.24]$ and an outer disk of $r\in [15.6,156]$ is used in disk self-gravity calculations. We emphasize that the effects of both disk self-gravity and an indirect acceleration term are important due to the relatively massive disk in HL Tau. In fact, simulations without including such terms could give erroneous results. A fixed boundary condition is used for gas at both inner and outer disk boundaries. For dust, we adopt an inflow (outflow) boundary condition at the inner (outer) boundary.

Partnership et al. (2015) identify seven pairs of dark and bright rings in the 1.0 mm image of the HL Tau disk, and among them four prominent dark rings could be considered as gaps, i.e., the dark rings D1, D2, and the adjoining D5 and D6. Based on a large body of previous literature on gap opening and formation as functions of disk gas temperature, planet mass, and disk viscosity (e.g., Crida et al. 2006; Li et al. 2009), we have performed tens of disk+planet hydro simulations in order to determine the likely planet mass values that could roughly match the multiple gap widths and depths as shown in the ALMA observations. Based on these simulations, we adopt a nominal model (run0), which has a set of planet mass parameters as ${M}_{{\rm{p}}}\approx 0.35,0.17,$ and $0.26{M}_{{\rm{Jup}}}$ at $13.1,33.0$, and 68.6 AU, respectively. The mass of the planet in the first gap is consistent with the mass of ∼0.3M_Jup as estimated from the gap structure (Kanagawa et al. 2015). These three planets are fixed at their orbital locations, and planetary radial migration is turned off.

To explore how our results could vary with some key parameters, we have carried out many additional runs. In Table 1, we list a few simulations where we reduce the masses for all three planets by half (run1), disable the disk self-gravity (run2), and reduce the disk viscosity (run3). Other changes such as the disk surface density and mass are not presented here. Most of the results presented in this paper use the snapshot at 4000 orbits (∼170,556 years at 10 AU) to generate simulated observations, except for run1 where the snapshot at 3400 orbits is used.

Table 1.  Simulation Parameters

Model	DSG	${M}_{{\rm{p}}}({M}_{{\rm{Jup}}})$	α	Orbit No.
run0	Y	0.35, 0.17, 0.26	1e-3	4000
run1	Y	0.17, 0.09, 0.13	1e-3	3400
run2	N	0.35, 0.17, 0.26	1e-3	4000
run3	Y	0.35, 0.17, 0.26	2e-4	4000

Download table as:  ASCIITypeset image

The dust opacities adopted in this work were calculated as in Isella et al. (2009) by assuming that dust grains are compact spheres made of astronomical silicates (Weingartner & Draine 2001), organic carbonates (Zubko et al. 1996), and water ice, with fractional abundances as in Pollack et al. (1994). Single grain opacities were averaged on a grain size distribution to obtain the mean opacity used in the radiative transfer model. We did that by adopting a typical MNR power-law size distribution (Mathis et al. 1977), $n(a)\propto {a}^{-3.5}$, between a minimum grain size of $5\times {10}^{-4}$ mm and a maximum grain size of 10 mm. The resulting dust opacity at the wavelength of 1 mm is 4.6 cm² g⁻¹, and it is dominated by dust grains with sizes between 0.1 and 0.2 mm. Based on this assumption, in our hydrodynamic models, we adopt dust grains with a size of 0.15 mm as best tracers of the dust emission at 1 mm.

To calculate synthetic maps of the disk continuum emission, we first need to calculate the dust temperature throughout the disk. This is controlled by micron-size dust particles that are well coupled with gas. We adopt a dust-to-gas ratio of 0.01 and use the gas surface density from simulations to interpolate the micron-size dust density distribution on a 3D spherical grid along with a scale height profile of, ${h}_{\mu {\rm{m}}-\mathrm{dust}}{(r)=1.0{\rm{AU}}\times (r/20{\rm{AU}})}^{1.25}$. We set up a thermal Monte Carlo run to calculate the dust temperature using the radiative transfer code RADMC-3D. The calculated dust temperature profile shows the typical two-layer vertical structure of passive irradiated circumstellar disks (Dullemond et al. 2002). In the surface layer, the temperature decreases from ∼280 K at 5 AU to ∼75 K at 100 AU. In the midplane, the temperature decreases from ∼110 K at 5 AU to ∼20 K at 100 AU. Next, using the obtained dust temperature, we calculate the millimeter thermal emission based on the density and opacity of our dust disk. We convert the surface density of 0.15 millimeter size dust from two-fluid simulations to a 3D spherical density profile using a scale height of, ${h}_{\mathrm{mm}-\mathrm{dust}}(r)={f}_{\mathrm{sett}}\times {h}_{\mu {\rm{m}}-\mathrm{dust}}(r)$, where ${f}_{\mathrm{sett}}=0.1$ is a parameter that accounts for the settling of 0.15 millimeter size dust toward the midplane. This leads to a scale height of ∼ 0.75 AU at 100 AU for relatively large dust particles, which is consistent with the findings by Pinte et al. (2015). Using our gaseous disk parameters, the estimated settling time of 0.15 mm dust is within 10⁵ years at the scale height of the gas disk, thus it is a reasonable choice of dust scale height, considering the age of the HL Tau disk. We then generate the synthetic images at the millimeter wavelength using the ray tracing method of RADMC-3D. As a last step, we "observe" our models and compare them with the ALMA observations in the Fourier domain, following the procedure described in Isella et al. (2009).

Figure 1 shows the surface density distributions of gas and 0.15 millimeter size dust in our nominal model. The gas and dust disks exhibit similar morphology, showing three gaps and three higher density rings/bands of gas created by the disk–planet interactions, along with visible spiral arm features. The rings/bands in the dust distribution show asymmetries associated with the spiral density features produced by the planets, especially in the inner and middle rings. The dust-to-gas surface density ratio (initially at 0.01) shows the evolution of dust+gas distributions while dust is subject to gas drag and is concentrated at the high pressure region, causing significant enhancements at the locations of rings. Dust filtration effects are expected in the gaps because the larger size dust particles are more concentrated in the rings. The deficit of dust at large radii $\gt 120\;{\rm{AU}}$ is due to the fact that dust particles drift radially inward and there is no supply of dust particles through the outer boundary.

Zoom In Zoom Out Reset image size

Since the HL Tau system is relatively young, the viscous timescale at disk radii $\gt 10\;{\rm{AU}}$ is longer than the age of the disk. So, the disk is not expected to be in a quasi-steady state. When monitoring the evolution of the ring and gap structures over time, we indeed find that the accretion rate in simulations is not steady and the depths of the gaps are still evolving (say between 4000 and 7000 orbits), though understandably the first gap is approaching a steady state faster.

Figure 2 shows the ALMA observation of HL Tau, the synthetic observations of our nominal model, and a residual map obtained by subtracting the model from the observations in the Fourier domain. Overall, our model reproduces the observations quite well. (Good agreement is obtained between our model and ALMA Band 6 and 7 observations, and we are only showing the Band 7 image due to its higher spatial resolution.) What is most remarkable is that the flux density level between the model and observations can be matched so closely, given the fact that we have used a generic smooth gaseous disk with a uniform initial dust-to-gas ratio throughout such a disk. The simulated observation of our model can explain pretty well the three major gaps at ∼13, 33, and 68 AU, but it does not reproduce the detailed structures of multiple rings at larger radii due to our choice of using only one planet at 68 AU. The large difference at the center region shown in the residual is due to our choice of 2.4 AU as the inner disk boundary. Interestingly, the spiral arms and the asymmetry in the rings that are visible in the gas and dust distributions (see Figure 1) are mostly smeared out in the simulated observation due to the high optical depth of the continuum emission at 850 μm.

Zoom In Zoom Out Reset image size

One of the key findings from the ALMA observations is that rings and gaps are somewhat eccentric. To quantify this feature in our model image, by using a face-on configuration, we find the local maxima (for the rings) or minima (for the gaps) in the azimuthal direction and fit the locations of these points with ellipses. We obtain the eccentricities of gaps as 0.246, 0.274, and 0.277, respectively, whereas we obtain the eccentricities of rings as 0.191, 0.079, and 0.157, respectively. This reproduces the trend seen in the actual ALMA observations, though quantitative comparisons are somewhat difficult. This reproduces the trend seen in the actual ALMA observations. Because we have calculated the eccentricities using our simulations which have much higher spatial resolutions than ALMA observations, simulations can yield more pronounced non-axisymmetric features on smaller scales, this makes more quantitative comparisons with observations difficult.

The left panel of Figure 3 shows the radial profile of the intensity along the major axis of the disk corresponding to a position angle of 138°. The flux densities of our model match the observation well, especially in the region around the first and second ring+gap. At the third ring and beyond, the relative difference between our model and the observation increases. Again, it is worth emphasizing that the relative heights among different regions were all obtained self-consistently through two-fluid gas+dust simulations along with the radiative transfer calculations.

Zoom In Zoom Out Reset image size

The right panel of Figure 3 shows the optical depth τ and the spectral index α of the dust emission calculated between 0.87 and 1.3 mm along the disk major axis. The dark rings are (partially) optically thin while the rest of the disk is mostly optically thick. For $\tau \ll 1$, α depends on the slope of the dust opacity β ($\alpha \sim 2+\beta$). Vice versa, for $\tau \gtrsim 1$, α approaches 2. The radial profile of α predicted by our model is also in remarkable agreement with the observations (see Figure 3 of Partnership et al. 2015).

Finally, note that since we used only one dust size in our simulations, the spectral index variation in our model is only determined by the properties of the dust disk (such as optical depth). In real observations, the radial variation of spectral index might also depend on the dust properties (such as size evolution) through the disk (Guilloteau et al. 2011; Testi et al. 2014, p. 339).

In Figure 4, we compare the other three models with the ALMA Band 7 image. Similar to Figure 3, we draw cross-cuts of the simulated observations of run1, run2, and run3 along the major axis of the disk and compare the obtained flux density profiles with the observation. Run1, where the planet masses are half of the nominal model, fails to produce the second gap and the depth of the first gap is too shallow compared to the observation. This is expected since less massive planets cannot open deep gaps. Run2, where the disk self-gravity is turned off, fails to produce the third gap. This is because, as the third planet gradually creates a gap, the gravitational force of the disk interior of the planet will try to pull the outer disk inward, further promoting the formation of a partial gap. Without such disk self-gravity, the gap is much shallower. In run3, where the disk viscosity is reduced by a factor of five, the gaps are too deep because the spiral shocks produced by planets in low viscosity disks are more effective in clearing the regions surrounding the planets.

Zoom In Zoom Out Reset image size

Note that the results shown in Figure 4 do not mean that we could not find a suitable fitting when parameters are different from our nominal model. On the contrary, it is likely that a right combination of planet masses, disk viscosity, temperature, mass, as well as the disk evolution time could be found to produce a reasonable fit to the current ALMA observation. The parameter studies in this subsection are meant to help us to gauge how sensitive the final outcome is to various parameters.