High-resolution ALMA Observations of HD 100546: Asymmetric Circumstellar Ring and Circumplanetary Disk Upper Limits

The most up-to-date and accurate distance estimate (110.02 pc from GAIA DR2) to the star is larger than the previously derived estimate (97 pc from Hipparcos), which was used to estimate the stellar parameters, therefore we refine the stellar parameters based on d = 110.02 pc. We adopt a PHOENIX model of the stellar photosphere (Hauschildt et al. 1999) with T_eff = 9800 K (Fairlamb et al. 2015) and $\mathrm{log}(g)\,=-4.0$, then it is scaled to the GAIA DR2 distance and to the de-reddened (A_V = 0.1 mag) V-band photometry. The integrated luminosity L_* is calculated from the model, which, combined with the aforementioned T_eff, are compared to the pre-main sequence stellar tracks by Siess et al. (2000). We employed the tracks with a depleted abundance of Z, because the source is depleted in refractory elements in its atmosphere (Folsom et al. 2012). This procedure yields a stellar mass and age of M_* = 2.2 ± 0.2 M_⊙ and $t={4.8}_{-1.1}^{+2.0}$ Myr, respectively. The reported uncertainties are obtained by propagating the uncertainties on the distance, A_V (±0.1), and T_eff (a conservative ±400 K).

The maps shown in Figure 1 reveal a bright ring between 20 and 40 au with a significant flux asymmetry, and an additional inner disk coincident with the position of the star. The inner disk is unresolved (<2 au in radius), with a peak flux of 2.60 ± 0.85 mJy beam⁻¹. Between the inner and outer disk there is a dark annulus with an average brightness of ∼1 mJy beam⁻¹, which is about 8× fainter than the (faintest section of the) central annulus of the ring emission.

In Figure 3 we compare the azimuthally averaged brightness temperature of the disk emission, for which we have calculated the deprojected radius using the position angle and disk inclination parameters obtained by Pineda et al. (2014). The same figure includes the parametric disk temperature profile from Panić et al. (2014) and the temperature profile of the millimeter-sized grains from the radiative transfer model from Pineda et al. (2014) (similar values of the midplane dust temperature at 50 au (≈60 K) were found by Bruderer et al. 2012). The figure shows that at every position in the disk the parametric disk temperature from Panić et al. (2014) is much higher than the observed values. However, the more detailed radiative transfer model reveals lower temperatures for the millimeter-sized particles, with an average temperature of T_d,mean = 53 K between 20 and 50 au. We use this average disk dust temperature, T_d,mean, as a best estimate of the disk emission in the following sections. Therefore, the disk dust continuum emission is optically thin at the positions of the young planet candidates, while the central part of the ring might be optically thick.

Zoom In Zoom Out Reset image size

We model the emission with a simple parametric model that includes two Gaussian rings, a central compact source, and a vortex (to account for flux asymmetries), as follows:

$$\begin{eqnarray}\begin{array}{rcl}F(r,{r}_{1},\theta ) & = & {F}_{r1}{e}^{-{\left(r-{r}_{r1}\right)}^{2}/2{\sigma }_{r1}^{2}}\\ & & +\ {F}_{r2}{e}^{-{\left(r-{r}_{r2}\right)}^{2}/2{\sigma }_{r2}^{2}}\\ & & +\ {F}_{g}\,{e}^{-{r}_{p}^{2}/2{\sigma }_{g}^{2}}/(2\pi \,{\sigma }_{g}^{2})\\ & & +\ {F}_{V}\,{e}^{-{\left(r-{r}_{v}\right)}^{2}/(2{\sigma }_{v}^{2})}{e}^{-{\left(\theta -{\theta }_{v}\right)}^{2}/(2* {\sigma }_{\theta }^{2})}\end{array}\end{eqnarray} \tag{ 1 }$$

where r and r_r are the radii calculated at the center of the ring and point source, respectively, and taking into account the inclination angle with respect to the sky (assumed to be the same for both coordinate systems). The first two elements in the model attempt to reproduce the main disk ring-like emission (which is not well reproduced by a single Gaussian profile) and are concentric, the third one describes the central unresolved source, and the fourth element describes a possible vortex.

We use GALARIO (Tazzari et al. 2018) to sample the model image on the same visibilities as the observations. The χ² is then calculated using the sampled visibilities, and then minimized in Python to find the optimal model. Also, we use the built-in options in GALARIO to perform 2D translations in the plane of the sky (δ R.A., δ decl.) and rotation (δ PA) of the parametric model. The best-fit model in deprojected coordinates, the observed model, and the residuals are shown in Figure 4, while the best parameter values are listed in Table 2. The observed model and residuals are produced by sampling the same visibilities as the data, and then performing the imaging in CASA.

Zoom In Zoom Out Reset image size

The combined vortex and double-ring model allows for a good fit of the image, although the residuals still show some structure, in particular close to the ring inner edge. However, none of these two rings or vortex correspond to the outer ring found by Walsh et al. (2014). We also note that the best fit confirms what is seen by eye: a significant offset between the central Gaussian source and the central position of the ring.

We generate two cuts, one through the disk major axis and one through the vortex maximum emission to investigate in more detail the ring morphology. Figure 5 shows the average flux along beam-wide strips along both directions. The profiles can clearly not be fitted with a single Gaussian flux distribution and they show significant asymmetries in the peak flux on both sides (≈15%–25%). Fitting the profiles with a superposition of five Gaussians provides a good fit, however. The best-fit parameters are summarized in Table 3.

Zoom In Zoom Out Reset image size

Table 3.  Multiple-Gaussian Fit^a

Component	Center	Peak Flux	σ
	(au)	(mJy beam−1)	(au)
NW–SE
#0	−31.1 ±  0.7	9.2 ±  0.4	8.9 ±  0.3
#1	−21.7 ±  0.2	7.6 ±  0.8	3.3 ±  0.4
#2	0.5 ±  0.3	3.3 ±  0.2	4.5 ±  0.4
#3	22.1 ±  0.2	5.8 ±  0.7	3.2 ±  0.5
#4	31.3 ±  0.6	9.8 ±  0.3	8.4 ±  0.2
NE–SW
#0	−54.3 ±  8.7	0.88 ± 0.09	14.7 ± 18.5
#1	−29.4 ±  0.3	14.11 ± 0.96	8.2 ±  0.3
#2	−0.2 ±  0.3	4.08 ± 0.08	3.4 ±  0.3
#3	22.3 ±  0.2	4.85 ± 0.67	4.2 ±  0.5
#4	30.3 ±  0.7	8.32 ± 0.39	10.0 ±  0.2

Note.

^aEach Gaussian is described as:

$$\begin{eqnarray*}&&f(x)={F}_{\mathrm{peak}}\,{e}^{-{\left(x-{x}_{\mathrm{center}}\right)}^{2}/2\left({\sigma }^{2}+{\sigma }_{\mathrm{beam}}^{2}\right)}.\end{eqnarray*}$$

Download table as:  ASCIITypeset image

We also attempted to fit the profiles with asymmetric Gaussians (see, e.g., Pinilla et al. 2017), but the results are rather poor and therefore not reported here.

Moreover, we decompose the deprojected map into its polar coordinates (radius and angle) to better understand the radial structure of the emission (see Figure 6). The plot confirms the radial asymmetry in the main ring, with a "slow" flux drop as the radius increases beyond the ≈30 au radius. On the other hand, the ring emission has a steeper inner edge and clear azimuthal asymmetry. However, all these structures are unrelated to a previous outer ring claimed by Walsh et al. (2014) at ≈190 au.

Zoom In Zoom Out Reset image size

Numerical simulations predict the presence of CPDs around young forming gas-giant planets (e.g., Szulágyi et al. 2014, 2017; Zhu et al. 2016). The expected disk sizes are supposed to be much smaller than the beam size of the observations presented here (∼3 au). Therefore, we expect any CPD emission to appear point-like in our data. However, we do not find any evidence for point-like emission close to or around the claimed protoplanet positions and place a strong 3σ detection limit of 198 μJy for an unresolved source.

The CO (3–2) emission (Figure 7) extends out to ≈27 (300 au), which is less extended than the emission detected with the ALMA Cycle0 data (Pineda et al. 2014; Walsh et al. 2014) because of the missing short baselines in our observations. Clearly the CO emission is much more extended than the continuum emission, which was already identified in the Cycle0 analysis. The first-moment (intensity weighted velocity) map is presented in Figure 7, overlaid with the continuum emission. The position–velocity (PV) diagram along the disk's major axis is presented in Figure 8. The Keplerian velocity profile for the HD 100546 system, with M_* = 2.2 M_⊙ and 42° inclination angle, reproduces the velocities at a distance >2'' from the star (red curve in Figure 8). For separations <2'', the velocities are better reproduced with an inclination angle of 32° (orange curve in Figure 8).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Given the non-detection of CPD emission toward HD 100546 b, we place upper limits on the mass or size of the CPD, depending on the assumption of optically thin or thick emission.

In the case of the CPD emission being optically thin, we estimate the total CPD mass via

$$\begin{eqnarray}&&{M}_{\mathrm{total}}=\displaystyle \frac{{d}^{2}\,{F}_{\nu }}{{B}_{\nu }({T}_{d}){\kappa }_{\nu }\,{f}_{d}},\end{eqnarray} \tag{ 2 }$$

where κ_ν is the dust opacity per dust mass, d is the distance to the source, F_ν is the observed flux, B_ν(T_d) is the blackbody function, and f_d is the dust-to-gas ratio. For the opacity we assume ${\kappa }_{\nu }=0.2(7\,\mathrm{mm}\,{\lambda }^{-1})$ cm² g⁻¹, which is consistent with the value used by Isella et al. (2014). This opacity assumes a dust composition and grain size distribution as in Isella et al. (2012). We note that this κ_ν is a factor ≈2× lower than that used by Beckwith et al. (1990) and Andrews et al. (2011), and therefore our mass upper limits are conservative. Finally, f_d is assumed to be 0.01. Given the emission upper limit determined in Section 4.4, we determine the CPD (dust and gas) mass upper limit in the optically thin case to be 1.44 M_⊕.

In the case of the CPD emission being optically thick, the disk radius is calculated from ${F}_{\nu }={B}_{\nu }({T}_{d}){\rm{\Omega }},$ where T_d is the dust temperature of the CPD and Ω is the area subtended on the sky (${\rm{\Omega }}=\pi {R}_{\mathrm{CPD}}^{2}/{d}^{2}$). Therefore, the radius is derived as

$$\begin{eqnarray}&&{R}_{\mathrm{CPD}}=\sqrt{\displaystyle \frac{{F}_{\nu }}{\pi {B}_{\nu }({T}_{d})}}\,d.\end{eqnarray} \tag{ 3 }$$

An upper limit for the CPD radius of 0.44 au is obtained using a CPD temperature equal to the mean dust temperature for the millimeter-sized particles in the radiative transfer model at that radius (T_d,mean = 53 K, see Section 4.1), while the radius is only 0.09 au for a temperature of 932 K, which is the estimated temperature from high-contrast imaging at the L- and M-bands (Quanz et al. 2015). Both numbers are much smaller than the 2.8 au radius of the Hill sphere expected for a 1 M_J planet at 53 au (HD 100546 b). Several studies have determined the CPD radius to be between 0.3 and 0.5 of the Hill radius (Quillen & Trilling 1998; Ayliffe & Bate 2012; Shabram & Boley 2013) A conservative CPD radius' upper limit of 0.44 au yields an upper limit for the planet mass of 47 M_⊕ (0.15 M_J).

Both cases, optically thin and thick limits, provide important constraints for gas-giant planet formation processes by constraining fundamental properties of CPDs. In addition, Zhu et al. (2016) provided predictions for the SEDs of CPDs, including fluxes up to the submillimeter regime based on the predicted flux at 870 μm, ∼800 μJy, which is almost a factor of 10× the noise level in the image. On the other hand, detailed numerical simulations and synthetic observations of CPDs carried out by Szulágyi et al. (2018) showed that at large separations from the central star, a small fraction of the CPDs (<R_Hill/3) are warmer than the CSD. Furthermore, based on the nominal CSD setup used by Szulágyi et al. (2018) the expected flux for a CPD around a 1 M_J planet at 52 au is ≈250 μJy, which is comparable to the upper limits reported here, and therefore is still consistent with the ALMA observations.

According to de Val-Borro et al. (2007), even a Neptune-mass planet can generate a vortex of Rossby-Wave Instability, so this is consistent with our planetary mass limit. How strong this vortex is depends on many factors apart from the planetary mass: dust-to-gas ratio, viscosity, magnetic field of the disk, etc. A detailed parameter study of various numerical simulations is needed for this system in order to constrain the planetary mass based on the vortex we observe, as has been done for IRS48 (Huang et al. 2018), which is beyond the scope of this work.

Figure 9 compares the results of a few studies that have provided upper limits for CPD masses (see also, Ricci et al. 2017). The CPD mass upper limit obtained for HD 100546 b in Section 5.1 is comparable to that reported by Ricci et al. (2017); however, our assumed dust opacity is smaller, therefore we re-scale the CPD mass estimate to the one used by Ricci et al. (2017) and plotted it using the dashed line in Figure 9. This sample includes systems covering a wide range of stellar (host) mass and environments. However, it consistently shows that the CPDs, in case they do exist, carry only a small amount of mass. This is at odds with several models that generate substantial CPDs to feed protoplanets (Gressel et al. 2013; Shabram & Boley 2013; Stamatellos & Herczeg 2015; Zhu et al. 2016). On the other hand, the current CPD mass limits are consistent with the "gas-starved" disk scenario proposed by Canup & Ward (2002, 2006), as well as the numerical simulations by Szulágyi (2017) that show a correlation between CPD mass and CSD mass.

Zoom In Zoom Out Reset image size

We calculate an upper limit for the potential planet's mass using Equation (7) from Szulágyi (2017),

$$\begin{eqnarray}&&{M}_{\mathrm{CPD}}\times {10}^{4}=3.17\,{M}_{\mathrm{CSD}}\,{M}_{p}-4.33\,{M}_{\mathrm{CSD}},\end{eqnarray} \tag{ 4 }$$

which relates the CPD, planetary (M_p), and CSD mass (all in units of M_J). This assumes that the planet is still accreting from the surrounding CSD, which is supported by the previous detection of L'- and M-band thermal emission. Assuming the optically thin CPD estimate case from above, we place a planetary mass upper limit of 1.65 M_J, using our CPD mass upper limit of 1.44 M_⊕ (0.0045 M_J) and the CSD mass of 50 M_J (Pineda et al. 2014). This upper limit on the planetary mass estimate is clearly less stringent than the one derived using the optically thick approximation in Section 5.1; however, the upper limit calculated using the relation between CPD and CSD does have a less strong assumption and might be more realistic than the one reported in Section 5.1.

The central emission is compact and represents the innermost circumstellar material. We fitted a single Gaussian over a 80 mas region with a total flux of 13.6 ± 1.0 mJy, a deconvolved FWHM of the major and minor axis of 80 ± 8 mas and 56 ± 6 mas, respectively, and with a position angle of 177° ± 14°.

We use Equation (2), the same dust properties used in Section 5.1, and a disk temperature of 300 K, to derive an inner disk mass of 15 M_⊕. The stellar accretion rate of the central star is estimated to be $\dot{{M}_{* }}={10}^{-{7.04}_{-0.15}^{+0.13}}\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$ (Fairlamb et al. 2015). Thus, the central disk depletion lifetime (${M}_{\mathrm{disk}}/\dot{M}$) is only 500 yr. Therefore, the disk must be replenished with material from the outer ring/disk (e.g., Pinilla et al. 2016).

Garufi et al. (2016) presented an unsharp masked version of the HD 100546 disk based on SPHERE/ZIMPOL polarimetric differential imaging data. This image shows the disk inner rim, a spiral to the NE, and an arm-like structure to the north. In Figure 10 we show the SPHERE Q_ϕ image with our ALMA continuum map overlaid in contours, while Figure 11 similarly compares it to the CO-integrated intensity. The SPHERE data are aligned to match the star position with the center of the compact dust continuum emission.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

This comparison confirms that the disk inner rim is well-traced by the SPHERE observations and by the ALMA observations (continuum and CO). The NE-spiral feature observed in the SPHERE data coincides with the central region of the ring in the continuum emission, which indicates that the spiral-like feature in scattered light does not have a counterpart in the midplane. However, this feature location and the general orientation are coincident with a spiral-like enhancement seen in the CO-integrated intensity. This coincidence might suggest that the spiral-like feature is real and present in the disk surface. This is consistent with the fact that small dust grains and gas are well coupled in those disk regions.

Based upon the low angular resolution CO (3–2) ALMA Cycle0 observations, a warp disk was claimed by Pineda et al. (2014) by comparing the PV diagram along the major axis. Since the stellar mass is well constrained, by overplotting the expected Keplerian curves it was clear that a single-disk inclination could not reproduce the observations (see Loomis et al. 2017, for a similar claim of a warp in AA Tau). Further analysis of the same data, but using a more complex modeling tool, also suggested the presence of a change in the disk inclination (Walsh et al. 2017). Our results also show a similar pattern in the centroid velocity map (Figure 7), where the inner disk region is slightly twisted.

The possible disk warp has been suggested before to explain different observations (Quillen 2006; Panić et al. 2010). Figure 8 shows the PV diagram along the disk major axis, which shows the same behavior seen from the previous low angular resolution, with the kinematics of the outer section of the disk (>2'') being better described by an inclination angle of ≈42°, while the inner section of the disk (<05) is better described by an inclination angle closer to ≈32°. This means that the whole disk is not well described by a single inclination angle.

Also, it has been proposed that departures from the Keplerian velocity field in the disk kinematics could provide an independent way to identify the presence of a CPD in HD 100546 (Perez et al. 2015). Unfortunately, the image fidelity and sensitivity of the CO (3–2) data presented here do not allow us to identify such a feature.