An asymmetric dust ring around a very low mass star ZZ Tau IRS

Core accretion planet formation theories predict that only terrestrial and icy planets can form around very low mass (VLM) stars ($\lesssim$0.1–0.2 $M_{\sun}$), and the formation of gas giant planets is essentially prohibited (e.g., Liu et al., 2019, 2020; Matsumoto et al., 2020; Miguel et al., 2020). The reasons are multifold. First, a low stellar mass results in longer dynamical timescale and longer growth timescale of dust, delaying the formation of planetary cores. Secondly, the mass of dust in protoplanetary disks, the building blocks of planets, is correlated with stellar mass (e.g., Ansdell et al., 2017). A low dust mass in disks around VLM stars leads to low mass of protoplanets, suppressing pebble accretion (Ormel & Liu, 2018) and preventing cores from ever reaching the critical mass to trigger runaway gas accretion. Last but not least, observed stellar accretion rates suggest that disks around VLM stars evolve faster than those around higher mass stars, further daunting the task of planet formation (Manara et al., 2012; Liu et al., 2020). In addition, water snowline has been proposed to be the preferred site of planetesimal and planet formation due to its enhanced local dust-to-gas ratio (Ormel, 2017). Around VLM stars, snowline is close to the star at $r\sim$0.1 au. Therefore compact systems with rocky planets like the TRAPPIST-1 (Gillon et al., 2016) are expected, while planet formation at $r\gtrsim$1 au is not (Miguel et al., 2020).

Indeed, statistical studies show that at small orbital separations rocky and icy planets are common around M dwarfs (e.g., Dressing & Charbonneau, 2015), while the occurrence rate of gas giant planets dives from 14% around 2 M_☉ stars to 3% around 0.5 M_☉ stars (Johnson et al., 2010). Around VLM stars, only a handful of gas giant planets have ever been found (e.g., GJ 3512; Morales et al., 2019), and their formation mechanism is a mystery. Since their host stars are mature mean sequence stars that have finished planet formation long before, it is unclear where and when do those planets form. Protoplanetary disks around VLM stars serve as excellent laboratories for investigating planet formation in these extreme environments.

The Atacama Large Millimeter/submillimeter Array (ALMA) has revealed a variety of sub-structures in protoplanetary disks, such as rings, gaps, cavities, spirals, and crescent structures. While rings/gaps are the most common (e.g., Andrews et al., 2018; Long et al., 2018; Cieza et al., 2020), roughly 10 disks show azimuthal asymmetries (e.g., Francis & van der Marel, 2020; van der Marel et al., 2020; Tsukagoshi et al., 2019; Pérez et al., 2018; Dong et al., 2018b; Cieza et al., 2017). Disks around VLM stars also harbor rings and gaps (e.g., Kurtovic et al., 2020); however, azimuthal asymmetries have not yet been reported. Asymmetric dust crescents may be local dust traps at gas pressure maxima (e.g., Raettig et al., 2015; Ragusa et al., 2017), thus ideal locations of planet formation (e.g., Chang & Oishi, 2010). Three possible origins of asymmetries have been discussed by van der Marel et al. (2020): (a) long-lived anticyclonic vortices at gap edges, possibly opened by companions (e.g., Raettig et al., 2015), (b) gas horseshoes at the edge of eccentric cavities curved by massive companions (e.g., Ragusa et al., 2017), and (c) part of spiral arms (van der Marel et al., 2016) triggered by companions. The main difference between the former two is in the mass of the companions: a vortex can be produced at the edges of gaps opened by planets as low mass as super-Earths (Dong et al., 2018a), whereas a horseshoe needs to be triggered by a much more massive companion, i.e., a brown dwarf. Overall, all three mechanisms suggest the presence of companions (planets) inside the cavity or gap.

In this paper, we report a ring and an asymmetric structure around the VLM star ZZ Tau IRS (spectral type: M4.5 to M5, $T_{\rm eff}$: 3015 to 3125 K, $L_{\rm star}$: 0.02 to 0.13 $L_{\sun}$, $M_{\rm*}$: 0.09 to 0.16 $M_{\odot}$, White & Hillenbrand, 2004; Andrews et al., 2013; Herczeg & Hillenbrand, 2014) in the Taurus star forming region at an assumed distance of 130.7 pc¹¹1We do not use the GAIA distance of 105.7 pc because the astrometric solution is poor with a value of the renormalized unit weight error (RUWE) of 2.490 (GAIA EDR3; Gaia Collaboration et al., 2020). taken from Akeson et al. (2019). The stellar luminosity is 0.017 to 0.113 $L_{\sun}$ (White & Hillenbrand, 2004; Andrews et al., 2013; Herczeg & Hillenbrand, 2014, scaled by the updated distance), with the large uncertainty mainly due to factors including flux calibration and extinction (see the discussions in Herczeg & Hillenbrand, 2014). With these $T_{\rm eff}$ and $L_{\rm star}$, the Baraffe et al. (2015) isochrone yields a stellar mass around 0.1 $M_{\odot}$ at an age of 1 to 12 Myr consistently. ZZ Tau IRS is a single star (i.e., not identified as a spectroscopic binary; Kounkel et al., 2019) while its companion ZZ Tau is 35^′′ away. Its spectral energy distribution (SED; Figure 5 in Appendix A) does not show any deficits at IR wavelengths $\lambda\lesssim$10 $\mu$m, and thus ZZ Tau IRS was not classified as having a cavity/gap structure in its disk. The large IR excess relative to the stellar flux in the SED (Figure 5) was interpreted as an edge-on disk (Furlan et al., 2011), and high-resolution optical spectroscopy of ZZ Tau IRS also suggested the existence of a close to edge-on disk based on its narrow emission lines from the optical jet (White & Hillenbrand, 2004). The total disk mass (gas $+$ dust) was estimated to be 0.014 $M_{\odot}$ by observations at millimeter wavelengths using the Submillimeter Array (Andrews et al., 2013; Akeson et al., 2019).

We used the ALMA archive data at band 7 (program ID: 2016.1.01511.S, PI: J. Patience), summarized in Table 1. The data were calibrated by the Common Astronomy Software Applications (CASA) package (McMullin et al., 2007), following the calibration scripts provided by ALMA. We performed a self-calibration of the visibilities. The phases were self-calibrated once, with fairly long solution intervals (solint = ‘inf’) that combined all spectral windows. To facilitate analyses later (§ 4), we shifted the observed dust continuum images to minimize the asymmetry relative to the phase center using the procedure described in Appendix B.

The final synthesized dust continuum image of the combined data is shown in Figure 1. In the CLEAN task, we set the $uv$-taper (0.4 $\times$ 500.0 M$\lambda$ at PA of 48°) with a robust parameter of $-$2.0 to obtain a nearly circular beam, and we did not use the ‘multi-scale’ option. The r.m.s. noise in the region far from the source is 268 $\mu$Jy/beam with a beam size of 254 $\times$ 248 mas at a position angle (PA) of $-$61.7°.

The ¹²CO $J=3\rightarrow 2$ line data (Table 1) were extracted by subtracting the continuum in the $uv$ plane using the uvcontsub task in the CASA tools. A line image cube with a channel width of 0.5 km/s was produced by the CLEAN task. The integrated line flux map (moment 0) and the intensity-weighted velocity map (moment 1) are shown in Figure 2. The channel maps at $+$1.0 to $+$13.0 km/s are shown in Figure 7 in the Appendix C. The r.m.s. noise in the moment 0 map is 61.4 mJy/beam$\cdot$km/s with a beam size of 283 $\times$ 156 mas at a PA of $-$50.2$\arcdeg$, and the r.m.s. noise in the moment 1 map is 9.56 mJy/beam at 0.5 km/s bin. The peak signal-to-noise ratio (SNR) is 24.1 in the channel map of $+$4.0 km/s.

Since the contrast of the asymmetry in the ZZ Tau IRS ring is of order unity, an imaging artifact such as a sidelobe of the bright ring could generate an apparent asymmetry. To characterize the disk structure, we performed forward modelling in which the observed visibilities are reproduced with a parametric model of the ring in the visibility domain.

The ZZ Tau IRS disk is assumed to consist of a ring and an asymmetry. The model ring consists of a narrow Gaussian ring (ring A) superposed on a wide Gaussian ring (ring B) as shown in Figure 3(a), and is described as follows:

where $r_{\rm peak,ringA}$, $\sigma_{r,\rm{ringA}}$, $r_{\rm peak,ringB}$, $\sigma_{r,\rm{ringB}}$, and $\delta$ are the radial peak positions and standard deviations of two Gaussian rings and a scaling factor between the two, respectively. The width of the ring (FWHM) is calculated as 2.355$\sigma_{\rm ring}$.

The asymmetry is defined as an elliptical Gaussian function in polar coordinates as follows:

where $r_{\rm peak,asym}$, $\theta_{\rm peak,asym}$, $\sigma_{r\rm{,asym}}$, and $\sigma_{\theta{\rm,asym}}$ are the radial peak position, azimuthal peak position, radial standard deviation, and azimuthal standard deviation, respectively. The flux for the asymmetry is normalized to $F_{\rm asym}$. The combined model image is magnified and rotated with an inclination ($i$) and a PA. The total flux for the combined model image is normalized to $F_{\rm total}$. The radial and azimuthal widths (FWHM) of the asymmetry are 2.355$\sigma_{r,\rm{asym}}$ and 2.355$\sigma_{\theta,\rm{asym}}$, respectively. In total, there are 13 free parameters in our model ($F_{\rm total}$, $F_{\rm asym}$, $r_{\rm peak,ringA}$, $\sigma_{r,\rm{ringA}}$, $r_{\rm peak,ringB}$, $\sigma_{r\rm{,ringB}}$, $\delta$, $r_{\rm peak,asym}$, $\theta_{\rm peak,asym}$, $\sigma_{r,\rm{asym}}$, $\sigma_{\theta{\rm,asym}}$, $i$, PA).

The model image was converted to complex visibilities with the public Python code vis_sample (Loomis et al., 2017), in which model visibilities are samples with the same ($u$, $v$) grid points as the observations. The model visibilities are deprojected²²2Visibilities are deprojected in the $uv$-plane using the following equations (e.g., Zhang et al., 2016): $u^{\prime}=(u\,{\rm cos}\,{\rm PA}-v\,{\rm sin}\,{\rm PA})\times{\rm cos}\,i,v^{\prime}=(u\,{\rm sin}\,{\rm PA}-v\,{\rm cos}\,{\rm PA}),$ where $i$ and PA are free parameters in our visibility analyses in § 4. 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 ln$L$ in visibility analyses is

where the subscript $j$ represents the $j$-th data. $V_{j}^{\rm obs}$, $V_{j}^{\rm model}$, and $W_{j}$ are observed and model visibilities and weights corresponding to 1/$\sigma_{j}^{2}$ ($\sigma_{j}$ is the rms noise of a given visibility. See details in CASA Guides³³3https://casaguides.nrao.edu/index.php/DataWeightsAndCombination.), respectively. The values of $V_{j}^{\rm obs}$ and $W_{j}$ (i.e., 1/$\sigma_{j}^{2}$) are provided in the measurement set of ALMA data. Our calculations used flat priors with the parameter ranges summarized in Table 2. We ran 3,000 steps with 100 walkers, and discarded the initial 1,000 steps as the burn-in phase based on trace plots in Figure 8 in Appendix D.

The fitting results with uncertainties computed from the 16th and 84th percentiles, the radial profile of best-fit surface brightness for the ring with a raw resolution, the best-fit model image, and the probability distributions for the MCMC posteriors are shown in Table 2, Figure 3(a), Figure 3(d), and Figure 9 in the Appendix E, respectively. We subtracted modeled visibilities from observed visibilities (Figure 10 in Appendix F), and made a CLEANed image (Figure 3e) in which we did not find significant residuals. To calculate the value of the reduced-$\chi^{2}$, we derived a factor between the weights and standard deviations (referred to as stddev) in the visibilities by calculating the standard deviations of the real and imaginary parts in 2 k$\lambda$ bins along the azimuthal direction in the visibility domain. The values of weights $W_{j}$ are provided in the measurement set of ALMA data while the standard deviations (stddev) are calculated in our fitting. The visibilities were deprojected with $i=$ 60$\arcdeg$ and PA $=$ 135$\arcdeg$. Figure 11 in Appendix G shows a comparison between the weights and standard deviations (stddev^-2). We found that the weights are typically 4.0$\times$ higher than the standard deviations. Finally, we calculated a reduced-$\chi^{2}$ of 1.91. We also conducted visibility analyses using a model with one Gaussian ring $+$ one asymmetry (i.e., $\delta=$ 0), and found significant residuals. Thus, the model with two Gaussian rings $+$ one asymmetry is necessary to reproduce the observed dust continuum image.

Though the observed dust continuum image in Figure 1(a) shows two peaks along the ring major axis, our visibility analyses suggest that only the one at PA $=$ 130$\arcdeg$ is real.

The radial ring position (i.e., the peak position) and its FWHM (Ring A $+$ Ring B) in the model image were measured to be 58 and 49 au in Figure 3(a). The contrast in the asymmetry, which is the ratio between the peak brightness at the asymmetry and that at its opposite side on the ring, is measured to be 3.0 in Figure 3(d).

We analyzed ALMA band 7 archive data (dust continuum and ¹²CO $J=3\rightarrow 2$ emission) of ZZ Tau IRS. The source has been classified as a VLM star with a mass of $\sim$0.1–0.2 $M_{\sun}$ in the literature. The observed ¹²CO $J=3\rightarrow 2$ kinematics suggest a dynamical mass of $\sim$0.1–0.3 $M_{\sun}$, consistent with previous estimates. The ALMA dust continuum image at a spatial resolution of 0$\farcs$25 shows a ring and an asymmetric structures around ZZ Tau IRS. Our visibility analysis confirmed that the observed asymmetry is real rather than apparent.

Our best-fit model of the ALMA dust observations shows that the asymmetry has a radial FWHM of 11.3 au, an aspect ratio of 7.2 (ratio of azimuthal FWHM to azimuthal FWHM), a brightness contrast of 3.0, and the asymmetry contributes 5% of the total disk flux. These metrics are comparable with asymmetries seen around intermediate mass stars ($\sim$2 $M_{\sun}$), implying a similar origin, despite ZZ Tau IRS being 10 times less massive. While future high spatial resolution observations in multiple wavelengths with ALMA and JVLA are needed to determine the exact origin of the asymmetry, the three leading hypotheses all invoke companions inside the gap.

The inner disk around ZZ Tau IRS, not visible in our observations, has been inferred to be edge-on based on the SED. Our observations showed that the inclination of the outer disk at $r=$ 50 au is 60°. These results imply that the ZZ Tau IRS inner disk may be misaligned relative to the outer ring, and it is possible that the misalignment (and the gap) is produced by a giant planet with $M_{\rm p}\gtrsim M_{\rm J}$ on an inclined orbit inside the gap. Future scattered light observations may test this hypothesis by searching for shadows on the outer disk cast by the inner disk.

Overall, the structures found in the ZZ Tau IRS disk suggest the presence of massive companions including gas giant planets at $\sim 10$ au. This is surprising as the core accretion theories have generally prohibited the formation of such planets around VLM stars (Miguel et al., 2020; Liu et al., 2020, § 1). Hence, ZZ Tau IRS will be a prime target in the investigation of planet formation around VLM stars.