ALMA Observations of the Asymmetric Dust Disk around DM Tau

Planets are believed to form in protoplanetary disks (e.g., Hayashi et al. 1985; Pollack et al. 1996). The early stages of planet formation can be identified by disk structures such as gaps and asymmetric structures via planet–disk interactions (e.g., Kley & Nelson 2012). Such structures have been reported in tens of dust disks with the Atacama Large Millimeter/submillimeter Array (ALMA; e.g., Andrews et al. 2018; Long et al. 2018; van der Marel et al. 2018), and roughly 10 disks show asymmetries (e.g., Dong et al. 2018b; Pérez et al. 2018; Tsukagoshi et al. 2019; Francis & van der Marel 2020; van der Marel et al. 2021). In particular, asymmetric structures such as the blob and crescent features possibly due to azimuthal gas pressure maxima (e.g., Raettig et al. 2015; Ragusa et al. 2017) could be ongoing planet-forming sites because gas pressure maxima efficiently trap dust grains, potentially leading to planetesimal formation. Two major possible origins of these asymmetries are discussed by van der Marel et al. (2021): long-lived anticyclonic vortices at gap edges possibly curved by companions (e.g., Raettig et al. 2015) and gas horseshoes due to eccentric cavities curved by massive companions (e.g., Ragusa et al. 2017). The main difference between the two is in the mass of companions: a vortex can be produced at the edges of gaps possibly opened by planets, whereas a horseshoe structure needs to be triggered by a much more massive companion, i.e., a brown dwarf. Therefore, in the vortex scenario, asymmetries could be the signposts of planets, while they are not in the horseshoe scenario. Though the origins of individual asymmetric disks have not been determined by current observations (van der Marel et al. 2021), investigating asymmetric disks could help the understanding of planet formation.

The object DM Tau (spectral type: M1, Kenyon & Hartmann 1995; T_eff: 3705 K, Andrews et al. 2011; M_*: 0.53 M_⊙, Piétu et al. 2007; distance: 145 pc, Gaia Collaboration et al. 2018) is a single-star system (Nguyen et al. 2012; Willson et al. 2016), and its protoplanetary disk has a weak asymmetry in the outer disk at r ∼ 20 au (Kudo et al. 2018). The spectral energy distribution (SED) of DM Tau shows a deficit at λ ∼ 1–10 μm, which was interpreted as the presence of a deep cavity at r ∼ 3 au around DM Tau by analyzing its SED (e.g., Calvet et al. 2002). Subsequent submillimeter interferometric observations with a beam size of 03 detected a 20 au dust ring (Andrews et al. 2011). These discrepancies were explained by recent ALMA long-baseline observations (Kudo et al. 2018); DM Tau has multiple rings at r = 3 and 20 au and low-contrast rings at r ≳ 60 au. The ¹²CO (2–1) gas disk around DM Tau has no cavity/ring structures (Kudo et al. 2018), possibly due to a high optical depth of ¹²CO, while other molecular species such as C₂H show a ring structure at r ∼ 80 au (Bergin et al. 2016). A candidate giant planet was reported at r ∼ 6 au by near-infrared (NIR) sparse aperture masking interferometry (Willson et al. 2016), but this detection has not yet been confirmed. The mass accretion rate of DM Tau, $\dot{M}$ ∼ 6 × 10⁻⁹ M_⊙ yr⁻¹ (Manara et al. 2014), is comparable with that of typical T Tauri stars (Najita et al. 2015). As small (submicron size) dust grains coupled with the gas flow into the central star, significant infrared excess at λ ∼ 1–10 μm should appear in the SED. Hence, the origin of DM Tau's 3 au cavity with both high mass accretion rate and strongly depleted dust grains in the cavity is still under debate (e.g., Manara et al. 2014; Kudo et al. 2018).

In this paper, we report follow-up observations of DM Tau with ALMA in cycle 6. Our original aim for these new observations was to confirm a weak asymmetry with a contrast of 20% in the inner ring at r ∼ 3 au reported by Kudo et al. (2018). Though we have not confirmed this asymmetry with our new cycle 6 observations, two blobs at the same radial location were newly identified in the outer ring at r ∼ 20 au. We also performed SED fitting to test the existence of dust grains inside the 3 au cavity around DM Tau.

The ALMA observations of DM Tau summarized in Table 1 were carried out with band 6 in the C43-9/10 configuration on 2019 June 5 UT under project 2018.1.01755.S, using 44 antennas with a baseline length extending from 83.1 m to 15.2 km. Since the short-baseline data are available in the ALMA archive (ID: 2013.1.00498.S; PI: L. Perez), we did not request these observations. The long-baseline data were taken with four spectral windows (SPWs): three with 128 channels spanning 1.875 GHz (31.25 MHz channel^–1) centered on 213.5, 216.3, and 228.0 GHz and one with 3840 channels spanning 1.875 GHz (448.3 kHz channel^–1, 0.64 km s⁻¹ velocity resolution) centered on the ¹²CO J = 2 → 1 rest frequency of 230.538 GHz. The bandpass and flux calibrator were the quasar J0423−0120, while the phase calibrator was J0431+1731. The mean precipitable water vapor was 0.8–0.9 mm during observations. The total on-source integration time was 134.3 minutes. The data were calibrated by the Common Astronomy Software Applications (CASA) package (McMullin et al. 2007) version 5.4.0–70, following the calibration scripts provided by ALMA. We separately conducted self-calibration of the visibilities in long- and short-baseline data. For long-baseline data, the phases were iteratively self-calibrated on solution intervals of 360 s combining all SPWs. However, as the self-calibration degraded the signal-to-noise ratio (S/N) due to the large amount of flagged data, we decided not to use self-calibrated long-baseline data in this paper. For short-baseline data, the phases were self-calibrated once with fairly long solution intervals (solint="inf") that combined all SPWs.

We combined our long-baseline data with short-baseline data to recover the missing emission at larger angular scales. We compared the visibility amplitudes at less than 200 kλ between the two data sets and confirmed their consistency. The images of both data sets were aligned in two manners as follows. ⁹

Method A—We separately synthesized the dust continuum images of short- and long-baseline data by CASA with the CLEAN task using a multifrequency deconvolution algorithm (Rau & Cornwell 2011). We then conducted ellipse isophote fitting at 30σ in the images of both types of data. The phase centers for both data sets were corrected to the centers of ellipse isophote fitting by fixvis in the CASA tools. To test whether the new phase center is the center of the disk, we subtracted the 180° rotated image in the visibility domain. This procedure corresponds to producing a synthesized image with only the imaginary part of the visibilities. Because the visibility is complex conjugate, the subtraction of the 180° rotated image is mathematically equal to setting the real part as zero and doubling the value of the imaginary part. In other words, the real part contains information on both symmetric and asymmetric structures of objects, whereas the imaginary part contains only information on asymmetries. Therefore, by synthesizing the image with only the imaginary part, we selectively remove only symmetric structures, and only asymmetric structures can be efficiently detected. ¹⁰ We searched the minimum rms in the central region of the images with shifting images relative to the center of ellipse isophote fitting in the visibility domain by the phase shift defined as $\exp \left[2\pi i\left(u\,{\rm{\Delta }}{\rm{R}}.{\rm{A}}.+v\,{\rm{\Delta }}\mathrm{decl}.\right)\right]$, where u and v are the spatial frequencies and ΔR.A. and Δdecl. are shift values. Figures 7 and 8 in Appendix A show dust continuum images synthesized with only the imaginary part, including the image with the minimum rms We found that the shift values with the minimum rms are (ΔR.A., Δdecl.) of (0, 0 mas) and (+4, −2 mas) relative to the center of ellipse isophote fitting in the long- and short-baseline data, respectively. Finally, the pointing tables for both data sets were corrected toward the images with the minimum rms by fixplanets in the CASA tools. The new phase centers ¹¹ of the long- and short-baseline data in ICRS coordinates are (4^h33^m4874961, 18^d10^m96177) and (4^h33^m4874792, 18^d10^m96819), respectively.

Method B—We also check the shift value with the minimum χ² of the imaginary part. The value of χ² is defined as $\sum \left({W}_{j}\,\mathrm{Im}{V}_{j}^{2}\right)$, where the subscript j represents the jth data. Here $\mathrm{Im}{V}_{j}$ and W_j are the visibilities in the imaginary part and weights, respectively. We found that the shift values with the minimum χ² are (ΔR.A., Δdecl.) of (+1, +1 mas) and (+5, −2 mas) relative to the center of ellipse isophote fitting in the long- and short-baseline data, respectively. Method A measures the rms in the region we are interested in, while method B measures the rms of entire visibilities in the imaginary part; thus, we rely on the results of method A in this paper.

The final synthesized dust continuum image of the combined data is shown in Figure 1. In the CLEAN task, we set the uv-taper to obtain a nearly circular beam pattern (Table 1), and we do not use the "multiscale" option. The rms noise in the region far from the object is 8.4 μJy beam^–1 with a beam size of 35.0 × 34.2 mas at a position angle (PA) of 675.

Zoom In Zoom Out Reset image size

The ¹²CO J = 2 → 1 line data in both the long- and short-baseline data were extracted by subtracting the continuum in the uv-plane with the uvcontsub task in the CASA tools. The combined line image cube with channel widths of 0.7 km s⁻¹ was produced by the CLEAN task. We also set the uv-taper to obtain the nearly circular beam pattern (Table 1). The integrated line flux map (moment zero) and the intensity-weighted velocity map (moment 1) are shown in Figure 2, while channel maps at −1.0 to +12.3 km s⁻¹ are shown in Figure 10 in Appendix C. The rms noise in the moment zero map is 3.0 mJy beam^–1 km s⁻¹ with a beam size of 45.6 × 45.4 mas at a PA of −779, while that in the moment 1 map at the 0.7 km s⁻¹ bin is 589 mJy beam^–1. The peak S/N is 15.9 in the channel map of +2.5 km s⁻¹.

Zoom In Zoom Out Reset image size

Figure 1 shows the dust continuum images of the DM Tau disk combining long- and short-baseline data at band 6. The entire disk is shown in Figure 1(a). As reported by Kudo et al. (2018), we confirmed three components in the disk: the inner ring at r ≲ 10 au, the outer ring at r ∼ 20 au, and the extended structure at r ≳ 60 au, as noted in Figure 1(d). We discovered that the extended structure consists of two faint rings at r ∼ 90 and ∼110 au with S/Ns of ∼6 and ∼5, respectively. This W-shaped double gap structure, found in a number of other disks, could be produced by a super-Earth-mass planet in a low-viscosity environment (Dong et al. 2017, 2018a; Pérez et al. 2019; Facchini et al. 2020).

The structure of the outer ring at an S/N of ∼70 is consistent with that in Kudo et al. (2018). In the image subtracting the 180° rotated image in Figure 7, we found two prominent asymmetries at PAs of ∼270° and ∼180° (hereafter blobs A and B) in the outer ring. These two blobs can be seen in the dust continuum image (Figure 1(b)), and the azimuthal profile of the outer ring at r ∼ 20 au is shown in Figure 1(e); i.e., these two locate at the same radial location. We also synthesized the dust continuum images with different imaging parameters (e.g., using the multiscale option) in the CLEAN task in Appendix B and Figure 9 and confirmed the presence of the two blobs at roughly the same location in the images with different parameters. The contrasts of these two relative to the opposite side of the ring are ∼1.1× to ∼1.2×. Note that blob A has already been reported by Kudo et al. (2018). The total flux density derived by visibility fitting in Section 4.1 is 94.74 ± 9.47 mJy, assuming a 10% uncertainty in absolute flux calibration, consistent with previous single-dish (109 ± 13 mJy; Beckwith et al. 1990) and ALMA (93.3 ± 0.5 mJy by visibility fitting; Kudo et al. 2018) observations. The peak brightness temperature in the outer ring, except for the blobs (i.e., at the inner edge of the outer ring, which is calculated from the best-fit modeled image in visibility fitting; Section 4.1), is 17.1 ± 1.2 K, assuming a 10% uncertainty in absolute flux calibration. The optical depth τ_ν is calculated with the relationship

$$\begin{eqnarray}&&{I}_{\nu }={B}_{\nu }\,{T}_{\mathrm{mid}}\{1-\exp (-{\tau }_{\nu })\},\end{eqnarray} \tag{ 1 }$$

where I_ν , B_ν , and T_mid are the intensity, Planck function, and midplane temperature, respectively. We use the midplane temperature profile with the simplified expression for a passively heated, flared disk in radiative equilibrium (e.g., Dullemond et al. 2001),

$$\begin{eqnarray}&&{T}_{\mathrm{mid}}(r)={\left(\displaystyle \frac{\phi {L}_{* }}{8\pi {r}^{2}{\sigma }_{\mathrm{SB}}}\right)}^{0.25},\end{eqnarray} \tag{ 2 }$$

where L_* is the stellar luminosity (taken as 0.36 L_⊙ from Manara et al. 2014), ϕ is the flaring angle (taken as 0.02), and σ_SB is the Stefan–Boltzmann constant. At the outer ring at r = 20 au, T_mid is estimated to be 21.6 K; thus, the optical depth τ_ν is calculated as ${1.3}_{-0.2}^{+0.3}$.

Since measuring the flux density of the inner ring in the synthesized image is not straightforward, potentially due to a contamination of the bright outer ring, we derive it in the best-fit modeled image in the visibility fitting (Section 4.1). The peak brightness temperature and the flux density of the inner ring within the cavity of the outer ring calculated with the best-fit modeled image are ${11.0}_{-0.5}^{+0.6}$ K and 1.74 ± 0.17 mJy, respectively, assuming a 10% uncertainty in absolute flux calibration. The flux density is similar to that in previous studies (1.33 mJy in Kudo et al. 2018) and converted to a total mass (gas + dust) of 0.04 M_Jup, assuming a distance of 145 pc, an opacity per unit dust mass κ_ν = 2.3 cm² g⁻¹ at 230 GHz (Beckwith & Sargent 1991), a temperature of 100 K, and a gas-to-dust mass ratio of 100. Note that Francis & van der Marel (2020) suggested a gas-to-dust mass ratio of 10⁴–10⁵ for the inner ring of DM Tau, assuming a viscosity α of 10⁻³; thus, the inner ring may be 2 or 3 orders of magnitude more massive.

Kudo et al. (2018) reported a possible asymmetry in the inner ring, i.e., 20% brighter emission in the northwest. The inner ring in our new data shows that the southeastern region is ∼20% (∼3σ) brighter than the northwestern region (Figure 1(b)). However, no such asymmetries in the inner ring can be seen in the image subtracting the 180° rotated image in Figure 7, where only asymmetric signals are contained (see explanations of method A in Section 2). ¹² Therefore, as these asymmetries could be the result of image reconstruction artifacts, more data are necessary to confirm both the asymmetries and morphological variability in the inner ring.

The inner ring shown in Figure 1(c) is likely to have a different PA than that of 1578 in the DM Tau system (Kudo et al. 2018). The bright part in the southern region in the inner ring is located at a PA of ∼180°. Since the beam shape is close to circular, the shape of the inner ring is unlikely to be affected by the beam elongation. We performed ellipse isophote fitting of the inner ring at the 8σ level and found that the PA of the inner ring is 1721 ± 42. The difference is significant at 3.5σ. Furthermore, Francis & van der Marel (2020) found that the PA of the inner ring in Kudo et al. (2018) is 141° ± 7° by Gaussian fitting in the image domain, which suggests that the PA of the inner ring in our new data varies compared with previous observations in Kudo et al. (2018). These results motivated us to perform visibility analyses to test whether or not the PA of the inner ring is different from that of the system and previous observations, because there is a possibility of image reconstruction artifacts in the inner ring (see Section 4.1).

The integrated line flux map of ¹²CO J = 2 → 1 in Figure 2(a) shows a single-peak symmetric structure with a peak flux of 39.2 mJy beam^–1 km s⁻¹ at 13σ, while Kudo et al. (2018) noted that the peak emission at 9σ is shifted with ∼20 mas toward north. ¹³ Assuming the positional errors are the values of the beam size divided by the S/N, the positional error of the peak ¹²CO emissions in Kudo et al. (2018) is ∼8 mas, i.e., an ∼2.5σ deviation. These positional shifts between two epochs could also be artifacts. More data are needed to confirm or reject the time variability in the inner ring. The intensity-weighted velocity map is also shown in Figure 2(b) and is consistent with that in Kudo et al. (2018).

4.1. Visibility Fitting

As the spatial scale of the inner ring is a few times the beam size, the structure in the inner ring corresponds to high spatial frequency components in the visibilities. In general, visibility data at a high spatial frequency are more sparse even in ALMA observations, potentially resulting in image reconstruction artifacts. To confirm the different PAs between the inner and outer rings inferred in Section 3, we performed forward modeling, in which observed visibilities are reproduced with a parametric model of the disk by utilizing all spatial frequency information.

In the literature, the parametric disk model is often described with Gaussian rings (e.g., Zhang et al. 2016; Pinilla et al. 2018). On the other hand, since the radial profile of the outer ring at r = 20 au around DM Tau shows an exponential profile (Figure 1(d)), we describe the surface brightness distributions of the disk in our model with a simple power-law radial profile with an exponential taper at the outside,

$$\begin{eqnarray}&&I(r)\propto \displaystyle \sum _{i=1}^{2}{\alpha }_{i}{\left(\displaystyle \frac{r}{{r}_{{c}_{i}}}\right)}^{-{\gamma }_{{\rm{p}},i}}\exp \left[-{\left(\displaystyle \frac{r}{{r}_{{c}_{i}}}\right)}^{{\gamma }_{{\rm{e}},i}}\right],\end{eqnarray} \tag{ 3 }$$

where α_i , ${r}_{{c}_{i}}$, γ_p,i , and γ_e,i are a scaling factor, a characteristic scaling radius, exponents of the power law, and the exponential taper, respectively. We divided the disk into two global components (i; Figure 3(a)) because the inner and outer rings (component 1) are roughly 1 order of magnitude brighter than the extended outer structure (component 2). In the radial direction, we have the following scaling factors:

$$\begin{eqnarray*}\begin{array}{rcl}{\alpha }_{1}(\mathrm{component}\ 1) & = & \left\{\begin{array}{lll}{\delta }_{2} & \mathrm{for} & {r}_{\mathrm{cav}}\lt r\lt {r}_{{\mathrm{gap}}_{1}}\\ 1 & \mathrm{for} & {r}_{{\mathrm{gap}}_{1}}\lt r\lt {r}_{{\mathrm{gap}}_{2}}\\ 0 & \mathrm{for} & {r}_{{\mathrm{gap}}_{2}}\lt r\end{array}\right.\\ {\alpha }_{2}(\mathrm{component}\ 2) & = & \left\{\begin{array}{lll}0 & \mathrm{for} & r\lt {r}_{{\mathrm{gap}}_{2}}\\ 1 & \mathrm{for} & {r}_{{\mathrm{gap}}_{2}}\lt r\lt {r}_{{\mathrm{gap}}_{3}}\\ {\delta }_{3} & \mathrm{for} & {r}_{{\mathrm{gap}}_{3}}\lt r\lt {r}_{{\mathrm{gap}}_{4}}\\ {\delta }_{4} & \mathrm{for} & {r}_{{\mathrm{gap}}_{4}}\lt r.\end{array}\right.\end{array}\end{eqnarray*}$$

At r < r_cav, we set a constant value with a depletion factor (δ₁) relative to the brightness at r = r_cav. We note that, since the radial profile at ∼60 au < r < ∼80 au is likely to be flat (Figure 1(d)), we added a pseudoring in this region (i.e., at ${r}_{{\mathrm{gap}}_{2}}\lt r\lt {r}_{{\mathrm{gap}}_{3}}$ in Figure 3(a)) to reproduce the nearly flat radial profile. Two components are normalized at $r={r}_{{\mathrm{gap}}_{1}}$ (Figure 3(a)). The total flux density (F_total) is also set as a free parameter. The disk inclination (i) and PA in the inner ring and the system (meaning the outer ring + the extended structure) are set as free parameters, i.e., i_{inner ring}, PA_{inner ring}, i_system, and PA_system. We fix the phase center.

Zoom In Zoom Out Reset image size

In addition to the above disk, we also add the model of blob A at a PA of ∼270° (Figure 7), defined as the elliptical Gaussian function in the polar coordinate as follows:

$$\begin{eqnarray*}\begin{array}{rcl}{\rm{A}} & = & \displaystyle \frac{{\cos }^{2}\,{\mathrm{PA}}_{\mathrm{blob}}}{{\sigma }_{\theta }^{2}}+\displaystyle \frac{{\sin }^{2}\,{\mathrm{PA}}_{\mathrm{blob}}}{{\sigma }_{r}^{2}},\\ {\rm{B}} & = & 2\left(\displaystyle \frac{1}{{\sigma }_{\theta }^{2}}-\displaystyle \frac{1}{{\sigma }_{r}^{2}}\right)\sin \,{\mathrm{PA}}_{\mathrm{blob}}\ \cos \,{\mathrm{PA}}_{\mathrm{blob}},\\ {\rm{C}} & = & \displaystyle \frac{{\sin }^{2}\,{\mathrm{PA}}_{\mathrm{blob}}}{{\sigma }_{\theta }^{2}}+\displaystyle \frac{{\cos }^{2}\,{\mathrm{PA}}_{\mathrm{blob}}}{{\sigma }_{r}^{2}},\\ {\rm{Z}} & = & {\rm{A}}{\left(\theta -{\theta }_{\mathrm{blob}}\right)}^{2}+{\rm{B}}(\theta -{\theta }_{\mathrm{blob}})(r-{r}_{\mathrm{blob}})\\ & & +{\rm{C}}{\left(r-{r}_{\mathrm{blob}}\right)}^{2},\end{array}\end{eqnarray*}$$

$$\begin{eqnarray}&&{I}_{\mathrm{blob}}(r,\theta )\propto \exp \left(-\displaystyle \frac{1}{2}{\rm{Z}}\right),\end{eqnarray} \tag{ 4 }$$

where PA_blob is the PA of the major axis of the elliptical Gaussian function in the polar coordinate, r_blob and θ_blob are the radial and azimuthal distances of blob A in the polar coordinate, and σ_θ and σ_r are standard deviations along the azimuthal and radial directions in the elliptical Gaussian function, respectively. The value of PA_blob is set to zero. The model image is finally rotated and magnified with PA_system and i_system, respectively. The total flux of blob A is normalized to F_blob. Note that the values of FWHM_r,blob and FWHM_θ,blob are equal to 2.355σ_r and 2.355σ_θ , respectively. Figure 11 in Appendix D shows the model image of blob A. Note that since blob B has a lower brightness, we do not include it in our model. In total, there are 25 free parameters in our model (F_total, F_blob, r_cav, ${r}_{{\mathrm{gap}}_{1}}$, ${r}_{{\mathrm{gap}}_{2}}$, ${r}_{{\mathrm{gap}}_{3}}$, ${r}_{{\mathrm{gap}}_{4}}$, δ₁, δ₂, δ₃, δ₄, γ_p,1, γ_e,1, γ_p,2, γ_e,2, ${r}_{{c}_{1}}$, ${r}_{{c}_{2}}$, i_{inner ring}, PA_{inner ring}, i_system, PA_system, r_blob, θ_blob, FWHM_r,blob, and FWHM_θ,blob).

The modeled disk image was converted to complex visibilities with the public python code vis_sample (Loomis et al. 2017), in which modeled visibilities are samples with the same (u, v) grid points with observations. The modeled visibilities are deprojected ¹⁴ with the system PA and i as free parameters. The fitting is performed with a Markov Chain Monte Carlo (MCMC) method in the emcee package (Foreman-Mackey et al. 2013). The log-likelihood function lnL in the MCMC fitting is

$$\begin{eqnarray*}\begin{array}{rcl}\mathrm{ln}L & = & -0.5\displaystyle \sum \left[f\,{W}_{j}\{{\left(\mathrm{Re}{V}_{j}^{\mathrm{obs}}-\mathrm{Re}{V}_{j}^{\mathrm{model}}\right)}^{2}\right.\\ & & +\left.{\left(\mathrm{Im}{V}_{j}^{\mathrm{obs}}-\mathrm{Im}{V}_{j}^{\mathrm{model}}\right)}^{2}\}\right],\end{array}\end{eqnarray*}$$

where the subscript j represents the jth data, and ${V}_{j}^{\mathrm{obs}}$, ${V}_{j}^{\mathrm{model}}$, and W_j are observed and modeled visibilities and weights, respectively. The value of f is a factor between the weights and the standard deviations in the visibilities. To estimate the value of f, we calculate the standard deviations of the real and imaginary parts in the 3 kλ bin along the azimuthal direction in the visibility domain. The visibilities were deprojected with i_system = 361 and PA_system = 1563. Figure 12 in Appendix E shows the comparison between weights and standard deviations, and we found that the typical value of f is 0.24. The weights are overestimated, or the noise is underestimated. Our calculations used flat priors with the parameter ranges summarized in Table 2. We ran 5000 steps with 100 walkers and discarded the initial 500 steps as the burn-in phase based on the trace plot in Figure 13 in Appendix F.

Table 2. Results of MCMC Fitting and Its Parameter Ranges

rcav	${r}_{{\mathrm{gap}}_{1}}$	${r}_{{\mathrm{gap}}_{2}}$	${r}_{{\mathrm{gap}}_{3}}$	${r}_{{\mathrm{gap}}_{4}}$	δ1	δ2	δ3	δ4
(au)	(au)	(au)	(au)	(au)
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)
3.22${}_{-0.12}^{+0.31}$	${21.25}_{-0.03}^{+0.09}$	${69.24}_{-1.42}^{+0.71}$	${83.07}_{-0.45}^{+0.51}$	${111.51}_{-1.40}^{+1.76}$	${0.24}_{-0.05}^{+0.03}$	${0.012}_{-0.000}^{+0.000}$	${2.08}_{-0.13}^{+0.16}$	${3.40}_{-0.23}^{+0.28}$
{0.0 .. 5.8}	{14.5 .. 29.0}	{55.1 .. 72.5}	{72.5 .. 92.8}	{98.6 .. 118.9}	{0.0 .. 1.0}	{0.0 .. 1.0}	{1.0 .. 10.0}	{1.0 .. 10.0}

γp,1	γe,1	γp,2	γe,2	rc1	rc2	iinner ring	PAinner ring	isystem	PAsystem
				(au)	(au)	(deg)	(deg)	(deg)	(deg)
(10)	(11)	(12)	(13)	(14)	(15)	(16)	(17)	(18)	(19)
${1.87}_{-0.11}^{+0.12}$	${1.63}_{-0.19}^{+0.13}$	${2.30}_{-0.11}^{+0.07}$	${1.04}_{-0.03}^{+0.03}$	${45.88}_{-5.65}^{+4.62}$	${118.72}_{-21.26}^{+7.23}$	${26.67}_{-1.41}^{+7.62}$	${157.26}_{-2.60}^{+3.53}$	${36.11}_{-0.23}^{+0.44}$	${156.31}_{-0.48}^{+0.32}$
{−1 .. 3}	{−1 .. 3}	{0 .. 3}	{0 .. 3}	{1.45 .. 145.0}	{1.45 .. 145.0}	{25.0-75.0}	{140.0–180.0}	{25.0–45.0}	{140.0–180.0}

Fluxtotal	Fluxblob	rblob	θblob	FWHMr,blob	FWHMθ,blob
(mJy)	(μJy)	(au)	(deg)	(au)	(deg)
(20)	(21)	(22)	(23)	(24)	(25)
${94.74}_{-1.74}^{+0.82}$	${222.0}_{-10.6}^{+28.4}$	${23.65}_{-1.64}^{+0.83}$	${9.37}_{-0.45}^{+0.51}$	${1.71}_{-0.49}^{+1.03}$	${45.38}_{-2.81}^{+3.23}$
{90.0 .. 110.0}	{100.0 .. 500.0}	{14.5-29.0}	{0.0–15.0}	{0.0–5.8}	{0.0–100.0}

Note. Brackets describe parameter ranges in our MCMC calculations.

Download table as:  ASCIITypeset image

The fitting results with errors computed from the 16th and 84th percentiles, the radial profile of the best-fit surface brightness, the best-fit modeled image, and the probability distributions for the MCMC posteriors are shown in Table 2, Figures 3(a) and (d), and Figure 14 in Appendix G, respectively. Though some parameters, such as r_blob, show double peaks in the probability distributions (Figure 14 in Appendix G), since the differences in the double peaks are small, we only show the results for the best-fit model in Figure 3. We subtracted modeled visibilities from observed ones and made a CLEANed image (Figures 3(e) and f). The reduced χ² is 1.7. We confirmed that the size of the inner cavity is 3 au in radius, which is consistent with the result in Kudo et al. (2018). Though Figure 1(c) implies that the inner ring is misaligned with the outer ring, the values of PA and i in the inner ring and the system are statistically the same within 3σ in our visibility analyses.

In the residual image in Figure 3(f), we found two additional significant residual signals, labeled as blobs C and D. Figure 4(a) shows the dust continuum image subtracting the 180° rotated image, overlaying the contours of blobs B–D in Figure 3(f). Though the counterparts of blobs B and C can be seen in Figure 4(a), blob C disappears in the image with different image shifts in Figure 7, e.g., the image with ΔR.A. =1 mas and Δdecl. = 1 mas. Hence, we consider blob B to be a real structure, while blob C might be an artifact. Furthermore, blob D has no counterpart in Figure 4(a); thus, we also consider blob D to be an artifact. By ellipse Gaussian fitting in the image domain, blob B is spatially resolved with a size of 119 × 78 mas (17.3 × 11.3 au).

Zoom In Zoom Out Reset image size

The residual image also suggests a large-scale asymmetry in the eastern part of the disk (Figure 3(e)). The residuals of the real and imaginary parts in Figures 3(b) and (c) also suggest deviations in the shorter baseline at a uv-distance of less than ∼500 kλ (corresponding to the scale of ∼04). To check this large-scale asymmetry, we compare the image subtracting the modeled disk (Figure 3(e)) and the image subtracting the 180° rotated image (Figure 4(a)) in Figure 15 in Appendix H. We found that both images in Figure 15 show large-scale asymmetry in the eastern part of the disk. Such large-scale asymmetries have been reported in other disk systems, potentially due to the shadow effect (e.g., Figure 4 in Facchini et al. 2020). The large-scale asymmetry around DM Tau (Figure 3(e); to be investigated elsewhere) has clumpy structures at the ∼3σ–4σ level and could potentially induce artificial clumps in the ring of blobs A and B. The spatial distribution of low-S/N clumps caused by thermal noise is expected to be random. If blobs A and B are indeed such artificial clumps, it is unlikely that they will reoccur in different observations. We reimaged dust continuum data in cycle 5 (Kudo et al. 2018) and found that both blobs can be seen at roughly the same locations (Figure 9; Appendix B). This suggests that blobs A and B are unlikely artificial clumps caused by noise.

We also found that the central region at r ≲ 3 au is not empty, because the value of δ₁ = 0.24${}_{-0.05}^{+0.03}$ is not zero at 4.8σ. The total flux at r ≲ 3 au is roughly 50 μJy. This result indicates the existence of an unresolved ring structure at r ≲ 3 au because the SED of DM Tau suggests the (nearly) empty cavity around the central star.

4.2. SED Fitting

The visibility analyses in Section 4.1 suggest significant millimeter emissions (∼50 μJy) in the central disk region within r ≲ 3 au. As DM Tau has no or very little NIR excess in the SED (e.g., Calvet et al. 2005), it has been believed that the dust grains inside the cavity are heavily depleted. The NIR excess mainly comes from small (submicron size) dust grains. To test the contributions of large (millimeter) dust grains in both the NIR excess in the SED and the millimeter flux (∼50 μJy) of the central emission, we conducted radiative transfer modeling using a Monte Carlo radiative transfer (MCRT) code (HO-CHUNK3D; Whitney et al. 2013). For this purpose, we put additional large and small dust grains inside the inner cavity at r = 3 au. The new cavity radius at r < 3 au is referred to as r_newcav hereafter. The fiducial surface density model and other models are shown in Figure 5(a).

Zoom In Zoom Out Reset image size

The MCRT code follows a two-layer disk model with small (up to micron size) dust grains in the upper disk atmosphere and large (up to millimeter size) dust grains in the disk midplane (e.g., D'Alessio et al. 2006). The modeled disk structure and dust properties in the MCRT code are described in our previous studies (Hashimoto et al. 2015): small dust grains from the standard interstellar medium dust model (a composition of silicates and graphites; a size distribution of n(s) ∝ s^−3.5 from ${s}_{\min }=0.0025$ to ${s}_{\max }=0.2$ μm) in Kim et al. (1994) and large dust grains (a composition of carbons and silicates; a size distribution of n(s) ∝ s^−3.5 from ${s}_{\min }=0.01$ to ${s}_{\max }=1000$ μm) from model 2 in Wood et al. (2002). The radial surface density is assumed to be a simple power-law radial profile similar to Equation (3),

$$\begin{eqnarray*}\begin{array}{rcl}{\rm{\Sigma }}(r) & = & \alpha {{\rm{\Sigma }}}_{0}{\left(\displaystyle \frac{r}{{r}_{c}}\right)}^{-q}\exp \left[-{\left(\displaystyle \frac{r}{{r}_{c}}\right)}^{2-q}\right],\\ \alpha & = & \left\{\begin{array}{lll}0 & \mathrm{for} & r\lt {r}_{\mathrm{newcav}}\\ \delta & \mathrm{for} & {r}_{\mathrm{newcav}}\lt r\lt 3\,\mathrm{au}\\ 1 & \mathrm{for} & 3\,\mathrm{au}\lt r\lt 20\,\mathrm{au}\end{array}\right.,\end{array}\end{eqnarray*}$$

where Σ₀ is the normalized surface density determined from the total (gas + dust) disk mass (M_disk) assuming a gas-to-dust mass ratio of 100, r_c is the characteristic radius of 50 au, q is the radial gradient parameter, and α is the scaling factor for the surface density. As the main purpose of our MCRT modeling effort is to reproduce the DM Tau SED at λ ≲ 10 μm, we set a grid size of r = 20 au in the code, i.e., inside the 20 au cavity. We set M_disk = 0.1 M_Jup to reproduce a flux of 1.74 mJy inside the 20 au cavity at 1.3 mm and fix q = 1. The scale heights (h) of large and small grains are assumed to vary as a power law with a radius, i.e., h ∝ r^p . To simplify, we assume p = 1.25 with a typical midplane temperature profile of T ∝ r^−0.5. We fix the scale heights of 1 and 4 au at r = 100 au for large (${h}_{r=100\mathrm{au}}^{\mathrm{large}}$) and small (${h}_{r=100\mathrm{au}}^{\mathrm{small}}$) dust disks, respectively, taken from Andrews et al. (2011) for the small dust disk. The mass fraction (f) of large dust grains in the total mass of dust grains is set to 0.9. The disk inclination is set to 30°. The HO-CHUNK3D code calculates the accretion luminosity at the star based on the mass accretion rate. Half of the flux is emitted as X-rays (which heat the disk), and half is emitted as stellar flux at a higher temperature. We set a mass accretion rate of $\dot{M}\,=$ 6 ×10⁻⁹ M_⊙ yr⁻¹ (Manara et al. 2014). In the code, we vary three parameters: r_newcav (same values for large and small dust disks), δ^large, and δ^small (where the superscript represents large or small dust grains), as shown in Table 3.

Table 3. Parameters in Our MCRT Modeling

Model	Rnewcav	δlarge	δsmall	rc	p	q	${h}_{r=100\,\mathrm{au}}^{\mathrm{large}}$	${h}_{r=100\,\mathrm{au}}^{\mathrm{small}}$	Mdisk	$\dot{M}$
	(au)			(au)			(au)	(au)	(MJup)	(M☉ yr−1)
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)	(10)	(11)
Fiducial	1.0	3 × 10−1	1 × 10−3	50	1.25	1	1	4	0.1	6 × 10−9
A	1.0	3 × 10−1	3 × 10−1	50	1.25	1	1	4	0.1	6 × 10−9
B	1.0	3 × 10−1	1 × 10−2	50	1.25	1	1	4	0.1	6 × 10−9
C	1.0	3 × 10−1	1 × 10−4	50	1.25	1	1	4	0.1	6 × 10−9
D	0.1	3 × 10−1	0	50	1.25	1	1	4	0.1	6 × 10−9

Download table as:  ASCIITypeset image

Figure 5(c) shows the SEDs for each model. Our fitting procedure includes three steps, as follows.

• 1.  
  We first set r_newcav = 1 au with varying values of the depletion factor δ with the same values in large and small dust disks to reproduce the flux of the central unresolved ring structure (∼50 μJy) at band 6, and we found that δ = 3 × 10⁻¹ is a reasonable parameter (model A).
• 2.  
  As the SED for model A largely emits at ∼10 μm (Figure 5(c)), we only vary the δ of small dust grains at 10⁻²–10⁻⁴ (models B and C and fiducial model). We found that the fiducial model well reproduces the DM Tau SED at λ ≲ 10 μm.
• 3.  
  We also set r_newcav = 0.1 au (model D). However, even for δ = 0 for small dust grains, this model largely emits in the NIR wavelengths.

In summary, the fiducial model could reasonably account for both the DM Tau SED and the flux of the central unresolved ring structure. This could be because the midplane temperature at r > 1 au is too low to emit at NIR, as shown in Figure 16 in Appendix I. Our modeling suggests that small dust grains inside the 3 au cavity are depleted, while large ones remain present.

In our two modeling efforts for visibility analysis (Section 4.1) and SED (Section 4.2) fitting, we assume independent modeled disks. To check the consistency of the surface brightness for the two models, we plot them at r < 20 au in Figure 17 in Appendix J. We confirm that the two radial profiles in our modeling efforts are consistent with each other.

5.1. Multiple Blobs

Asymmetric structures referred to as blobs in this paper are interpreted as dust trapped in gas vortices (e.g., Raettig et al. 2015) or gas horseshoes (e.g., Ragusa et al. 2017). For the latter, only one blob is expected at the edge of a cavity; thus, DM Tau's multiple blobs would not be this case. Furthermore, a gas horseshoe is expected at the edge of a cavity opened by massive companions (i.e., brown dwarfs). Previous NIR sparse aperture masking interferometry (Willson et al. 2016) and radial velocity measurement (Nguyen et al. 2012) have ruled out the presence of such companions in the cavity, disfavoring the gas horseshoe origin of DM Tau's blobs. Note that though van der Marel et al. (2021) mentioned the spiral structures as a third origin of asymmetries, since scattered-light images of DM Tau are not available, it is unclear whether the spiral structures are responsible for DM Tau's blobs.

Theoretical works of vortices (e.g., Ono et al. 2018) predict that multiple small vortices at similar radial locations tend to merge into one large vortex within hundreds of orbits. The existence of multiple blobs would therefore indicate the youth of vortices. Another interpretation of multiple blobs is the destruction of a large-scale vortex due to the heavy-core instability (Chang & Oishi 2010), triggered by a close-to-unity dust-to-gas mass ratio in the core of a vortex. Recent numerical simulations by Li et al. (2020) show multiple small blobs in the ring after the destruction of a large vortex. As the orbital number (system age divided by Keplerian orbital period at the radial location of the blob) of DM Tau's blob A is more than 10⁴, DM Tau's blobs may be the outcome of vortex destruction. In this case, the dust-to-gas mass ratio is expected to be close to unity. For DM Tau, however, the azimuthally averaged dust-to-gas mass ratio is estimated to be ∼0.01 in the outer ring at r ∼ 20 au (L. Francis et al. 2021, in preparation), disfavoring the scenario. A third possibility is that the azimuthal position of trapped dust depends on the dust size (Baruteau & Zhu 2016); centimeter-sized dust grains are trapped ahead of the gas vortex center in the azimuthal direction, while millimeter-sized dust grains concentrate closer to the vortex center. To examine this possibility, multiple wavelength observations to measure the spectral index sensitive to the grain size are necessary.

We also compare the shape of DM Tau's blob A with other asymmetries. Asymmetric structures have been reported in roughly 10 protoplanetary disks (e.g., Kraus et al. 2017; Pérez et al. 2018; Tsukagoshi et al. 2019; Francis & van der Marel 2020; van der Marel et al. 2021). Table 4 summarizes the physical quantities of blobs, mainly relevant to their morphology. ¹⁵ We found that DM Tau's blob A is located at the smallest radial location in the sample. Furthermore, its aspect ratio (the azimuthal width divided by the radial width) is the largest. Figure 6 shows these observational results, which place DM Tau's blob A in a novel parameter space of asymmetries.

Zoom In Zoom Out Reset image size

Table 4. Physical Quantities of Blobs and Other Properties

Object	Radial Location	Radial Width	Azimuthal Width	Aspect Ratio	References
	(au)	(au)	(deg; au)
(1)	(2)	(3)	(4)	(5)	(6)
DM Tau	24	1.7	45; 18.7	11	1
AB Aur	170	96	122; 361.8	3.8	2
CQ Tau	50	19	59; 51.5	2.7	2
	50	19	59; 51.5	2.7	2
HD 34282	137	110	52; 124.3	1.1	3
HD 34700	155	72	64; 173	2.4	4
HD 142527	180	81	155; 486.7	6.0	2
HD 143006	74.2	11	38.4; 49.7	4.5	5
IRS 48	70	29	58; 70.8	2.4	2
MWC 758	50	7.5	49; 42.7	5.7	2
	90	15	47; 73.8	4.9	2
SAO 206462	79	20	96; 132.3	6.6	2
SR 21	55	19	82; 78.7	4.1	2
	58	19	165; 166.9	8.8	2

Note. The radial width of IRS 48 is calculated with 2.17σ_r because the radial profile of IRS 48 was found to be best fit with a fourth power in the 2D Gaussian in van der Marel et al. (2021). The aspect ratio is defined as the azimuthal width divided by the radial width. References for blob information on radial locations and shapes: (1) this work, (2) van der Marel et al. (2021), (3) van der Plas et al. (2017), (4) Benac et al. (2020), (5) Pérez et al. (2018).

Download table as:  ASCIITypeset image

5.2. Central Emission

We present new ALMA observations of DM Tau, including dust continuum images at 1.3 mm and ¹²CO J = 2 → 1 emission maps. The dust continuum data with better sensitivity than previous studies in Kudo et al. (2018) reveal (a) multiple ring structures at r = 3, 20, 90, and 110 au and (b) two blobs at PAs of ∼180° and ∼270° in the outer ring at r ∼ 20 au. To characterize the inner ring regions at r ≲ 10 au, we analyze the dust continuum emission in the visibility domain and conduct modeling efforts using the MCRT code (Whitney et al. 2013). Consequently, we find the unresolved 1 au ring inside the 3 au inner ring. Furthermore, model disks with different disk inclinations and PAs between the inner ring at r ∼ 3 au and the whole system suggest that the inner ring is statistically aligned to the whole system within 3σ.

The two blobs of DM Tau have the low contrast of ∼1.1× to ∼1.2× relative to the outer ring. These two blobs are located at the smallest radial location among 11 asymmetric disks. Furthermore, the aspect ratio of blob A is the largest. These observational results place DM Tau's blob A in a novel parameter space of asymmetries. The origin of two blobs is not determined by our observations. The early phase of vortex formation (e.g., Ono et al. 2018), the destruction of the large-scale vortex (e.g., Li et al. 2020), or double continuum emission peaks with different dust sizes (Baruteau & Zhu 2016) could account for the multiple blobs. Future high spatial resolution observations in multiple wavelengths with ALMA and JVLA would help identify the origin.

We also found significant emission with a lower-limit mass of 0.4 M_Earth inside the 3 au cavity. As the DM Tau SED shows negligible NIR excess, the inside of the 3 au cavity around DM Tau has been believed to be a dust-free region. By fitting both the DM Tau SED and the flux of the central emission using the MCRT code, our modeling shows that inside the 3 au cavity, there is an additional 1 au dust ring where large (millimeter size) dust grains are less depleted than small (submicron size) dust grains. This would be due to efficient grain growth in the 1 au dust ring. Furthermore, our modeling indicates that the disk midplane temperature in the 1 au ring is less than the typical water sublimation temperature of 150 K (e.g., Notsu et al. 2016), which suggests that only small icy planets (i.e., mini-Neptunes) could form even in the terrestrial planet formation regions around DM Tau.

We thank the anonymous referee for a helpful review of the manuscript. This paper makes use of the following ALMA data: ADS/JAO.ALMA#2018.1.01755.S, ADS/JAO.ALMA#2017.1.01460.S, and ADS/JAO.ALMA#2013.1.00498.S. ALMA is a partnership of the ESO (representing its member states), NSF (USA), and NINS (Japan), together with the NRC (Canada), NSC (Taiwan), ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by the ESO, AUI/NRAO, and NAOJ. This work is based in part on archival data obtained with the Spitzer Space Telescope, which was operated by the Jet Propulsion Laboratory, California Institute of Technology, under a contract with NASA. Support for this work was provided by an award issued by JPL/Caltech. This work was supported by JSPS KAKENHI grant Nos. 19H00703, 19H05089, and 19K03932. Y.H. is supported by the Jet Propulsion Laboratory, California Institute of Technology, under a contract with NASA. R.D. is supported by the Natural Sciences and Engineering Research Council of Canada through a Discovery Grant and the Alfred P. Sloan Foundation through a Sloan Research Fellowship.

Software: vis_sample (Loomis et al. 2017), HO-CHUNK3D (Whitney et al. 2013), CASA (McMullin et al. 2007), emcee (Foreman-Mackey et al. 2013).

Figures 7 and 8 show the dust continuum images synthesized with only the imaginary part of long- and short-baseline data, respectively, by shifting 1 mas in R.A. and decl. The image is shifted relative to the center of ellipse isophote fitting (see Section 2) in the visibility domain by the phase shift defined as $\exp \left[2\pi \left(u\,{\rm{\Delta }}{\rm{R}}.{\rm{A}}.+v\,{\rm{\Delta }}\mathrm{decl}.\right)\right]$, where u and v are the spatial frequencies and ΔR.A. and Δdecl. are shift values. The rms values are measured inside the black dotted circles with radii of 03 and 05 in Figures 7 and 8, respectively. The 1σ noise in Figures 7 and 8 is 11.6 and 89.0 μJy beam^–1, respectively, measured in the region far from the object.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We synthesized dust continuum images with different imaging parameters in the CLEAN task to examine the robustness of blobs A and B. We also reimaged the dust continuum data obtained in cycle 5 (Kudo et al. 2018) to check whether the two blobs are present in different data sets. To minimize the effect of the beam elongation, we adjusted the uv-taper parameters to obtain a nearly circular beam. Figure 9(a) shows the dust continuum image using the same parameters as Figure 1 but using the multiscale option with scales of [0, 5, 10, 15, 25] in the CLEAN task. We confirmed that the image is consistent with Figure 1. Figures 9(b) and (c) are the same as Figure 9(a) but with robust = 2.0 and −2.0. Figure 9(d) shows the dust continuum images obtained in cycle 5, with robust = 2.0. Since we do not see significant differences in the images with or without the multiscale option in Figure 9(a) and Figure 1, the multiscale option was not used in Figures 9(b)–(d). We confirmed the presence of blobs A and B at roughly the same locations in all cases.

Zoom In Zoom Out Reset image size

Figure 10 shows the ¹²CO J = 2 → 1 channel maps at −1.0 to +12.3 km s⁻¹.

Zoom In Zoom Out Reset image size

Figure 11 shows the model image of blob A used in MCMC calculations in Section 4.1.

Zoom In Zoom Out Reset image size

Figure 12 shows the comparison of the values of weight and the standard deviations in real and imaginary parts. The standard deviations (stddev) are calculated in each 3 kλ bin in the deprojected visibilities of real and imaginary parts with an i of 36° and a PA of 156°. We found that the values of weight are typically 3.85× higher than 1 stddev^{– 2} of the real and imaginary parts.

Zoom In Zoom Out Reset image size

Figure 13 shows the trace plot of 100 walkers of the parameter r_gap1 in our MCMC calculations (Section 4.1). The burn-in phase is set as the initial 500 steps.

Zoom In Zoom Out Reset image size

Figure 14 shows the histograms of the marginal distributions of the MCMC posteriors for 25 free parameters calculated in our modeling in Section 4.1.

Zoom In Zoom Out Reset image size

Figure 3(e) in our visibility analyses in Section 4.1 suggests that residual image is asymmetric in the eastern part. To check this asymmetry, the residual image is compared with the image subtracting the 180° rotated image in Figure 15. Both images show that the eastern part of the disk is brighter. Thus, DM Tau could have a disk with large-scale asymmetry.

Zoom In Zoom Out Reset image size

Figure 16 shows a profile of the midplane temperature of large and small dust grains in the fiducial model calculated in the MCRT modeling in Section 4.2.

Zoom In Zoom Out Reset image size

Figure 17 shows the azimuthally averaged radial profile at r < 20 au generated by the MCMC model fitting and MCRT modeling in Sections 4.1 and 4.2, respectively, to test the consistency of the surface brightness in the two modeling efforts.

Zoom In Zoom Out Reset image size