PLANETARY SYSTEM FORMATION IN THE PROTOPLANETARY DISK AROUND HL TAURI

In this section, we describe a disk model that is utilized for estimating planet masses.

We adopt a simplified version of a disk model used in Kwon et al. (2011), who consider the volume density of gas (${\rho }_{g}(r,z)$) and the corresponding disk temperature ($T(r,z)$), both of which are functions of r and z. On the other hand, we integrate ${\rho }_{g}(r,z)$ for the vertical direction and assume that the disk temperature can be modeled by the two characteristic temperatures (T_s(r) and T_i(r)), where T_s(r) represents the temperature of disk surfaces and T_i(r) denotes the temperature at the mid-plane. Thus, the surface density distribution of gas (${{\rm{\Sigma }}}_{g}(r)$) becomes essentially comparable to the famous similarity solution (Pringle 1981), and can be written as

$$\begin{eqnarray}{{\rm{\Sigma }}}_{g}(r) & = & {\rho }_{0}\sqrt{\pi }H(r){\left(\displaystyle \frac{r}{{r}_{c}}\right)}^{-p}\\ & & \times \mathrm{exp}\left[-{\left(\displaystyle \frac{r}{{r}_{c}}\right)}^{7/2-p-q/2}\right],\end{eqnarray} \tag{ 2 }$$

where ${\rho }_{0}$ is a constant, H(r) is the pressure scale height, p is a fitting parameter, and q = 0.43 is the exponent that involves the radial temperature profile (see Table 1). Note that ${\rho }_{0}$ is computed for a given value of the disk mass (M_disk). For both T_s and T_i, a simple power-law distribution is assumed, that is, ${T}_{s}={T}_{s0}{(r/{r}_{c})}^{-q}$ and ${T}_{i}={T}_{i0}{(r/{r}_{c})}^{-q}$, following Kwon et al. (2011). Since we are currently examining sub-mm observations, T_i plays an important role in reproducing the observations.

Table 1.  Summary of Model Parameters^a

${T}_{i0}$ b	p	q	β	${\kappa }_{0,@230\mathrm{GHz}}$	rin	rc	gtd	Mdisk	i
Free parameter	1.064	0.43	0.729	0.01	2.4 AU	78.9 AU	100	0.1349 ${M}_{\odot }$	40.3°

Notes.

^aSee Kwon et al. (2011) and references therein. ^bWe treat ${T}_{i0}$ as a free parameter because the temperature cannot be well constrained at mm/sub-mm wavelengths (see Figure 3 in Kwon et al. 2011).

Download table as:  ASCIITypeset image

Based on the above disk model, Kwon et al. (2011) perform the Markov Chain Monte Carlo analysis, and derive a best set of model parameters with which the observed quantities, such as images, are well reproduced. We simply use the results of their best fit case (see Table 1). One may then wonder how reasonable it is to interpret new ALMA observations using previous results derived from observations with moderate spatial resolutions. We examine this issue by calculating the surface brightness at all of the three band to identify how well they match with the ALMA observations. Figure 4 shows the resultant surface brightness derived from Kwon's model. Note that for given values of ${{\rm{\Sigma }}}_{g}(r)$, T_i(r), the disk inner radius (T_in), and the disk inclination (i), the surface brightness was computed, following Dullemond et al. (2001) (see their Equation (21)). The results indicate that the Kwon et al. model can accurately reproduce the ALMA observations in all bands except for the gap regions. Thus, it can be concluded that disk models that can well reproduce the previous observations with moderate spatial resolutions are still valuable to characterize the global structure of disks.

Zoom In Zoom Out Reset image size

We now estimate the mass of hypothetical bodies by which the observed gaps can be opened. To proceed, we make use of the surface brightness profile obtained in the above section. Note that it is assumed that the central star is located at the observed continuum peak. We estimate their masses by following two approaches: the Hill radius analysis and a more detailed method developed from the balance of angular momentum transfer. The advantage of the first method is simple and independent of disk quantities such as surface density, temperature, or turbulent viscosity, which are sometimes poorly constrained. The advantage of the second method is a good capture of disk physics, wherein angular momentum transfer is triggered both by a planet and by gas disks due to viscosity and pressure.

First, we adopt the Hill radius analysis and examine the resultant mass of the bodies. Assuming that the gap width (${\rm{\Delta }}{r}_{p}$) is comparable to the Hill radius of the bodies, ${\rm{\Delta }}{r}_{p}$ is related to the mass of the bodies (${M}_{p,\mathrm{hill}}$) as

$$\begin{eqnarray}&&{\rm{\Delta }}{r}_{p}\simeq {\left(\displaystyle \frac{{M}_{p,\mathrm{hill}}}{3{M}_{*}}\right)}^{1/3}{r}_{p},\end{eqnarray} \tag{ 3 }$$

where ${M}_{*}=0.55{M}_{\odot }$ is the stellar mass (Kwon et al. 2011), and r_p is the distance from the continuum peak to the bodies. Since the radial profile of the continuum image at $\lambda =0.87\;{\rm{mm}}$ (band 7) shows the most detailed gap structure due to the highest spatial resolutions, we use the data at band 7 (see Figure 3(a)). Table 2 summarizes the resultant value of ${M}_{p,\mathrm{hill}}$ at each gap. We find that the object at the innermost orbital radius is more than 80 M_J, which is so massive that the body can be viewed as a substellar object rather than a planetary companion. All of the other bodies have masses approximately comparable to Jovian planets. Note that while we have assumed so far that the third and fourth gaps are individual gaps, we have also examined a possibility that these two are actually classified as one large gap. This may be implied by the radial surface brightness profile (see Figure 4); it may show that the third and fourth gaps are just a sub-substructure of one large gap. If this is the case, ${\rm{\Delta }}{r}_{p}={16.7}_{-0.5}^{+0.4}$ AU and ${r}_{p}={71.8}_{-1.0}^{+1.4}$ AU, which ends up with ${{\rm{M}}}_{p,\mathrm{hill}}$ that becomes $22\pm 3$M_J. As expected, the value is so large that the body is more likely to be a substellar object.

Table 2.  Summary of the Gap Properties and the Mass of Hypothetical Bodies

Gap	rp	${\rm{\Delta }}{r}_{p}$	${M}_{p,\mathrm{hill}}$ a	${M}_{p}$ b
	[AU]	[AU]	[MJ]	[MJ]
1	${13.5}_{-0.4}^{+0.4}$	${5.0}_{-0.1}^{+0.1}$	${88.8}_{-5}^{+5}$	0.85
2	${32.4}_{-0.4}^{+0.6}$	${4.1}_{-0.2}^{+0.2}$	${3.6}_{-0.6}^{+0.7}$	0.61
3	${65.2}_{-0.9}^{+1.3}$	${6.2}_{-0.8}^{+0.6}$	${1.5}_{-0.5}^{+0.5}$	0.62
4	${77.2}_{-0.7}^{+0.8}$	${4.5}_{-0.4}^{+0.4}$	${0.3}_{-0.1}^{+0.1}$	0.51

Notes.

^aThe planet mass derived from Hill radius (see Equation (3)). ^bThe planet mass derived from the method of Kanagawa et al. (2015).

Download table as:  ASCIITypeset image

Based on the above discussion, the Hill radius analysis results in the value of masses that are quite large, especially for the innermost gap. Additionally, the mass of $22\pm 3$ M_J obtained from the assumption of one large gap formed by the third and fourth gaps is unlikely since the recent mid-infrared ${L}^{\prime }$ band observation with sufficiently high sensitivity shows no signature of such a large body at the location (Testi et al. 2015). These results may suggest that the method itself is over-simplified, and hence leads to a significant overestimate of planetary masses. Accordingly, we now adopt a more detailed method and re-estimate the resultant masses.

There are a number of dedicated studies in which the gap structure is investigated by taking account of the balance of angular momentum transfer in gas disks under the action of tidal torque arising from planets (e.g., Lin & Papaloizou 1979; Goodman & Rafikov 2001). Here, we follow a method derived by Kanagawa et al. (2015). In the method, the gap depth (rather than the gap width) is required to estimate the planet mass, which can be written as (see their Equation (7))

$$\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}\\ & & \times {\left(\displaystyle \frac{{h}_{p}}{0.1}\right)}^{5/2}{\left(\displaystyle \frac{\alpha }{{10}^{-3}}\right)}^{-1/2},\end{eqnarray} \tag{ 4 }$$

where M_p is a planet mass, ${{\rm{\Sigma }}}_{p}$ is the surface density at the planet orbital radius, ${{\rm{\Sigma }}}_{0}$ is the surface density where no planet exists, h_p is a disk aspect ratio, and α is viscosity.

Adopting the disk model given by Kwon et al. (2011) coupled with some disk parameters used by Kanagawa et al. (2015), we compute the mass of hypothetical bodies that can curve the observed gaps. Table 2 tabulates the outcome. The results show that the computed value of the mass at all of the gaps becomes comparative to or less than a Jupiter mass, which is consistent with the estimate done in other studies (Dipierro et al. 2015; Tamayo et al. 2015). As an additional confirmation we have also computed the mass of bodies, following a method proposed by Crida et al. (2006) (the results are not shown in this paper). For this case, the gap width is needed to estimate the mass as done in the Hill radius analysis. We have found that, within the uncertainty of disk viscosity and stellar mass, the resultant values are all consistent with those summarized in Table 2. This indicates that, while some disk quantities may be uncertain, the second approach would provide a more accurate value in the mass estimate. We therefore will use the value of masses derived from a method introduced in Kanagawa et al. (2015) in the following section.

We now discuss what kind of formation mechanism(s) of the bodies can be plausible to account for the observed gap structure. Currently, there are two planet formation mechanisms that are being actively investigated: the core accretion (CA), and the GI scenarios. In the CA, which is the widely accepted theory of planet formation, the subsequent formation of planets is divided into two phases (e.g., Pollack et al. 1996; Ida & Lin 2004; Hasegawa & Pudritz 2014). In the first phase, planetary cores form via the so-called runaway and oligarchic growth (e.g., Wetherill & Stewart 1989; Kokubo & Ida 1998), while in the second phase, the surrounding gas in disks accretes onto the formed cores (Mizuno 1980; Stevenson 1982; Ikoma et al. 2000). In the GI, contrastingly, planets form in gravitationally unstable disks as stars form in gravitationally unstable molecular clouds (e.g., Boss 1997; Gammie 2001; Rafikov 2005; Cai et al. 2006; Kratter et al. 2010; Meru & Bate 2010). In this case, massive disks are needed to trigger the GIs there, which results in fragmentation of the disks and the subsequent formation of massive clumps. In the following discussion, we attempt to examine the possibility that the GI might actually work in this system. This is naturally motivated by the observational results that HL Tau is currently viewed as the intermediate phase between Classes I and II, and therefore its disk may be massive enough to trigger the GI.

We first compute the minimum value of the surface density of gas (${{\rm{\Sigma }}}_{g,\mathrm{min}}$) that is required to trigger GIs. Based on the so-called Toomre Q condition (Toomre 1964), ${{\rm{\Sigma }}}_{g,\mathrm{min}}$ can be written as (Rafikov 2005)⁴ in the optically thick limit,

$$\begin{eqnarray}{{\rm{\Sigma }}}_{g,\mathrm{min}} & \simeq & 2.6\times {10}^{6}\;{\text{g cm}}^{-2}\\ & & \times {\left(\displaystyle \frac{{M}_{*}}{{M}_{\odot }}\right)}^{7/8}{\left(\displaystyle \frac{{r}_{p}}{1{\rm{AU}}}\right)}^{-21/8}{\left(\displaystyle \frac{\kappa }{{\kappa }_{0}}\right)}^{1/4},\end{eqnarray} \tag{ 5 }$$

where $\kappa ={\kappa }_{0}{(\nu /{\nu }_{0})}^{\beta }$ is the opacity. We set ${\kappa }_{0}=0.01$ cm² g⁻¹, ${\nu }_{0}=230\;{\rm{GHz}}$ and $\beta =0.729$, following the disk model of Kwon et al. (2011). The resultant value of ${{\rm{\Sigma }}}_{g,\mathrm{min}}$ is summarized in Table 3. One immediately finds that ${{\rm{\Sigma }}}_{g}$ at $r\;=\;$ 65.2 AU and r = 77.2 AU exceeds the value of ${{\rm{\Sigma }}}_{g,\mathrm{min}},$ indicating that the minimum condition of the GI is satisfied at these two positions. To further visualized this trend, Figure 5(a) shows the region in which the GIs can operate (see the hatched region).

Zoom In Zoom Out Reset image size

Table 3.  Summary of Results

Gap	${{\rm{\Sigma }}}_{g}({r}_{p})$ a	$T({r}_{p})$ a	${{\rm{\Sigma }}}_{g,\mathrm{min}}({r}_{p})$	${T}_{\mathrm{min}}({r}_{p})$	GI	Mgap	Migration
	[g cm ${}^{-2}$]	[K]	[${\rm{g}}\;{\mathrm{cm}}^{-2}$]	[K]		[MJ]	Hillb	Cridac
1	45	30	$1.4\times {10}^{3}$	64	No	0.1	Type II	0.14
2	49	21	$1.4\times {10}^{2}$	22	No	0.2	Type II	0.44
3	34	15	29	8	Yes	0.36	Type II	0.68
4	26	14	18	6	Yes	0.41	Type II	1.10

Notes.

^aThe value derived from the disk model of Kwon et al. (2011). ^bThe results from Hill radius analysis (see Equation (8)). ^cThe results from a method proposed by Crida et al. (2006).

Download table as:  ASCIITypeset image

However, we cannot fully conclude yet that the GI can really operate there. This is because it is important to examine the cooling condition that assures that disks fragment into self-gravitating clumps (e.g., Gammie 2001). To proceed, we calculate the minimum disk temperature (T_min(r)) that is required to satisfy the cooling condition. The actual formula can be given as (Rafikov 2005), in the optically thick limit,

$$\begin{eqnarray}{T}_{\mathrm{min}}(r) & \simeq & 1.4\times {10}^{3}{\rm{}}\;{\rm{K}}{\left(\displaystyle \frac{r}{1{\rm{}}{\rm{AU}}}\right)}^{-6/5}\\ & & \times {\left(\displaystyle \frac{{M}_{*}}{1{M}_{\odot }}\right)}^{2/5}{\left(\displaystyle \frac{\kappa }{{\kappa }_{0}}\right)}^{2/5}{\left(\displaystyle \frac{{{\rm{\Sigma }}}_{g}}{10{\text{g cm}}^{-2}}\right)}^{2/5}.\end{eqnarray} \tag{ 6 }$$

The resultant value of ${T}_{\mathrm{min}}({r}_{p})$ is summarized in Table 3. We find that ${{\rm{\Sigma }}}_{g}({r}_{p})$ and $T({r}_{p})$ at $r=65.2\;{\rm{AU}}$ and $r=77.2\;{\rm{AU}}$ exceed ${{\rm{\Sigma }}}_{g,\mathrm{min}}$ and ${T}_{\mathrm{min}}({r}_{p})$. As done in Figure 5(a), the region where GIs can be feasible in terms of ${T}_{\mathrm{min}}({r}_{p})$ is denoted by the hatched region in Figure 5(b). Our analysis therefore suggests that, when the observed gaps originate from the presence of unseen bodies, the GI is very likely to be a mechanism for the gaps at r = 65.2 and 77.2 AU. Note that our finding is consistent with the results of Kwon et al. (2011), which indicate that the HL Tau disk may be gravitationally unstable in the outer region, especially the region between 50 and 100 AU in radius.

In the case of HL Tau, therefore, the GIs can occur in the region beyond $r\sim 52$ AU, where both of the conditions are simultaneously met.

As discussed above, the gaps located in the outer part of the disk may be related to the GIs, which can form massive clumps that are possibly the progenitors of planets. How about the gaps located in the inner part of the disk (at r = 13.5 and 32.4 AU)? Is there any possibility that an origin of these gaps can be discussed in the same framework as above? To this end, we focus on planetary migration that is induced mainly by the disk gas. Here, we attempt to discuss how all of the observed gaps can be established, coupling GIs with planetary migration.

Planetary migration arises from the tidal resonant interactions between massive bodies and the surrounding gas disks, which leads to the efficient angular momentum transfer between them (e.g., Kley & Nelson 2012, for a review). There are mainly two modes in migration that are classified by the mass of migrating bodies: type I and type II (e.g., Ward 1997; Nelson et al. 2000; Tanaka et al. 2002; Paardekooper et al. 2010; Hasegawa & Ida 2013). Type I migration is effective for low-mass bodies such as cores of gas giants and terrestrial planets, and the direction of the migration is very sensitive to the disk structure such as the surface density and the temperature. Type II migration becomes important for massive bodies such as Jovian planets in our solar system, and its direction is generally regulated by the bulk gas flow in the disks, that is, spiraling toward the central star.

We first determine which type of migration can take place in the HL Tau system. This can be done by estimating the so-called gap-opening mass (M_gap) above which planets are sufficiently massive to open up a gap in their gas disk via disk–planet interactions, and hence the planets will undergo type II migration. In principle, M_gap can be formulated either by the Hill radius analysis or by a disk viscosity analysis (e.g., Crida et al. 2006; Hasegawa & Pudritz 2014). We first adopt the Hill radius analysis to compute M_gap, which can be written as

$$\begin{eqnarray}{M}_{\mathrm{gap}} & \simeq & 2.7{M}_{\oplus }\\ & & \times {\left(\displaystyle \frac{r}{1{\rm{}}{\rm{AU}}}\right)}^{3/2}{\left(\displaystyle \frac{{M}_{*}}{1{M}_{\odot }}\right)}^{-1/2}{\left(\displaystyle \frac{T}{100{\rm{K}}}\right)}^{3/2}.\end{eqnarray} \tag{ 7 }$$

Note that the HL Tau disk would have a high value of disk accretion rates, equivalent to a high value of the disk turbulence. As a result, our estimate may become conservative in the sense that the Hill radius analysis would provide a lower value of M_gap. The resultant value of M_gap is summarized in Table 3. We find that type II migration would be the dominant mode of migration at all of the gaps.

We now can examine how valid the Hill radius analysis is by adopting a disk viscosity analysis. We rely on the criterion of type II migration derived from Crida et al. (2006), which can be written as (see their Equation (15))

$$\begin{eqnarray}&&R=\displaystyle \frac{3}{4}\displaystyle \frac{h}{{R}_{{\rm{H}}}}+\displaystyle \frac{50}{\hat{q}\mathrm{Re}}\lesssim 1,\end{eqnarray} \tag{ 8 }$$

where h is the pressure scale height of the disk, R_H is the Hill radius, $\hat{q}$ is the mass ratio of the planet to the central star M_p/M_⋆, and Re is the Reynolds number ${\hat{r}}_{p}^{2}{{\rm{\Omega }}}_{p}$/ν, where ${\hat{r}}_{p}$ is the distance to the planet from the central star, ${{\rm{\Omega }}}_{p}$ is the angular velocity of the planet, and ν is the disk viscosity, respectively.

Based on the criterion, planets can undergo type II migration when the value is less than unity, and they can experience type I migration when it becomes larger than unity. Table 3 summarizes the resultant value. We find that all of the hypothetical bodies may undergo Type II migration, except the outermost body. The resultant R for the outermost one marginally exceeds unity. In addition, it can naturally be expected that the disk viscous parameter (α) would have a relatively large uncertainty, which can result in an overestimate of R. Thus, it can be concluded that both the approaches lead to a roughly consistent consequence, which suggests that type II migration can be effective for hypothetical bodies.

If type II migration is the case, we can develop an intriguing conjecture as follows: the hypothetical bodies that might play some role in opening up the two inner gaps ($r=13.5\;{\rm{AU}}$ and r = 32.4 AU) may have formed originally in the outer part of the disk where the GIs can operate. Then, the bodies may subsequently have undergone type II migration. Since the direction of type II migration is inward for most of the cases, they would  have eventually arrived at their current locations. In other words, the outer part of the disk would act as the main site of forming all of the massive bodies, and all of the observed gaps in the HL Tau system might be established by the following inward migration of the bodies.

The unprecedented high spatial resolution images of HL Tau derived from the ALMA observations have a substantial potential to revolutionize our picture of planet formation. As discussed above, the observed multiple gap structure may be a consequence of the coupling of GIs with inward planetary migration. Nonetheless, it is obvious that further observational and theoretical studies are needed to examine our conjecture. For instance, the current image quality can be further improved by conducting long-baseline observations. Moreover, optically thin line observations such as ¹³C¹⁸O or deuterium molecules are demanded to investigate the gas distribution near the mid-plane of the disk. We also cannot fully exclude the hypothesis that the estimated Jovian bodies in the analysis might form via CA.