Gaps and Rings in an ALMA Survey of Disks in the Taurus Star-forming Region

Feng Long; Paola Pinilla; Gregory J. Herczeg; Daniel Harsono; Giovanni Dipierro; Ilaria Pascucci; Nathan Hendler; Marco Tazzari; Enrico Ragusa; Colette Salyk; Suzan Edwards; Giuseppe Lodato; Gerrit van de Plas; Doug Johnstone; Yao Liu; Yann Boehler; Sylvie Cabrit; Carlo F. Manara; Francois Menard; Gijs D. Mulders; Brunella Nisini; William J. Fischer; Elisabetta Rigliaco; Andrea Banzatti; Henning Avenhaus; Michael Gully-Santiago

doi:10.3847/1538-4357/aae8e1

1. Introduction

Characterizing the structure of protoplanetary disks is crucial to understand the physical mechanisms responsible for disk evolution and planet formation. Given the typical size (∼100 au) of protoplanetary disks (see review by Williams & Cieza 2011), spatially resolving disks in nearby star-forming regions (≲200 pc) requires observations with subarcsecond resolution. Disk observations with the Atacama Large Millimeter/submillimeter Array (ALMA) have revealed a variety of disk substructures from thermal emission of millimeter-sized grains, dramatically changing our view of protoplanetary disks. Axisymmetric gaps and rings are the most frequently seen substructures and have been observed in disks around HL Tau, TW Hydra, AA Tau, DM Tau, AS 209, Elias 2–24, V1094 Sco, HD 169142, HD 163296, and HD 97048 (ALMA Partnership et al. 2015; Andrews et al. 2016; Isella et al. 2016; Walsh et al. 2016; Zhang et al. 2016; Cieza et al. 2017; Loomis et al. 2017; van der Plas et al. 2017; Dipierro et al. 2018; Fedele et al. 2018; van Terwisga et al. 2018), in both young and evolved systems around T Tauri and Herbig stars. Large azimuthal asymmetries also emerge in some systems (Brown et al. 2009; van der Marel et al. 2013), as well as spiral arms (Pérez et al. 2016; Tobin et al. 2016). The origin of these substructures and their role in the planet formation process are still widely debated.

In typical protoplanetary disks, millimeter-sized particles are expected to undergo fast radial drift toward the central star due to aerodynamic drag with the gas, resulting in severe depletion of millimeter-sized dust grains at large radii (Weidenschilling 1977; Birnstiel & Andrews 2014). However, this picture is contradicted by high-resolution images of millimeter-sized particles that are distributed over distances of tens or hundreds of au from the central star (see reviews by Testi et al. 2014 and Andrews 2015). Assuming that the rings revealed from ALMA are related to variations of dust density, the presence of rings indicates that inward drift of large dust grains (millimeter-sized) can be stopped or mitigated at specific radii. The physics that generates the rings therefore contributes to the persistence of millimeter-sized dust grains at large radii, even after a few megayears of disk evolution (e.g., Gonzalez et al. 2017). The accumulation of dust in these regions might trigger efficient grain growth, thereby acting as an ideal cradle for forming planets (Carrasco-Gonzalez et al. 2016). A fundamental question then is what triggers the dust accumulation into ring shapes in disks and its connection to planet formation.

The mechanisms proposed to produce ring-like substructures in disks may be categorized into those related to disk physics and chemistry and those related to planet–disk interactions. When caused by disk physics and chemistry, the presence of a gap may trace the beginning of subsequent planet formation. Some of the disk-specific mechanisms that can generate gaps and rings include zonal flows induced by magnetorotational instabilities (Johansen et al. 2009), dead zones where gas accretion is regulated by spatial variations of the ionization level (Flock et al. 2015), grain growth around condensation fronts (Zhang et al. 2015), ambipolar diffusion-assisted reconnection in magnetically coupled disk wind systems in the presence of a poloidal magnetic field (Suriano et al. 2018), disk self-organization due to nonideal MHD effects (Béthune et al. 2017), suppressed grain growth with the effect of sintering (Okuzumi et al. 2016), large-scale instabilities due to dust settling (Lorén-Aguilar & Bate 2016), and secular gravitational instabilities regulated by disk viscosity (Takahashi & Inutsuka 2016).

The disk gaps and rings could also be induced by interactions between the disk and planet(s) within the disk. On the one hand, a massive planet (≳Neptune mass) embedded in the disk forms a gap in the gas density structure around its orbit, leading to the formation of a pressure bump outside the planet orbit, trapping large dust grains into rings and forming deep dust gaps (e.g., Lin & Papaloizou 1986; Pinilla et al. 2012a; Zhu et al. 2012). On the other hand, a planet with a mass as low as 15 M_⊕ is able to slightly perturb the local radial gas velocity, inducing a "traffic jam" that forms narrower and less depleted gaps (Rosotti et al. 2016). Lower-mass planets can produce deep dust gaps without affecting the local gas structure (Fouchet et al. 2010; Dipierro et al. 2016; Dipierro & Laibe 2017). Depending on the local disk conditions (e.g., temperature and viscosity) and planet properties, a single planet can also create multiple gaps (Bae et al. 2017; Dong et al. 2017).

Connecting these rings to the known distribution of exoplanets is challenging. Statistical studies of exoplanets reveal a higher occurrence rate of giant planets around solar-type stars than M dwarfs, while this trend is not seen for smaller planets (see a recent review by Mulders 2018). For more massive stars, the formation of the cores of giant planets is expected to be more efficient (Kennedy & Kenyon 2008). The surrounding disks would then have more material to build more massive planets, as suggested by the stellar–disk mass scaling relation from recent disk surveys (e.g., Andrews et al. 2013; Pascucci et al. 2016). If gaps are carved by giant planets, then deeper and wider gaps should be more prevalent around solar-mass stars than around stars of lower mass, although this picture would be complicated by any mass dependence in disk properties (e.g., low-mass planets can more easily open gaps in inviscid disks; Dong et al. 2017). In the case of ice lines, gaps should form at certain locations determined by the disk temperature profile, which broadly scales with stellar luminosity.

The analysis of gap and ring properties with stellar/disk properties should help us to discriminate between these different mechanisms. However, the small number of systems observed at high spatial resolution (∼0 farcs 1) to date limits our knowledge about the origins of disk substructures. Moreover, the set of disks imaged at high resolution is biased to brighter disks, many with near/mid-IR signatures of dust evolution, and collected from different star-forming regions and thus environments. These biases frustrate attempts to determine the frequency of different types of substructures, how these substructures depend on properties of the star and disk, and any evolution of substructures with time.

In this paper, we investigate properties of substructures in 12 disks, selected from a sample of 32 disks in the Taurus star-forming region that were recently observed at high resolution with ALMA. The paper is organized as follows. In Section 2, we describe our ALMA Cycle 4 observations and sample selection for the 12 disks. In Section 3, we present the modeling approach for disk substructures in the visibility plane and the corresponding model results. We then discuss in detail the stellar and disk properties for the 12 disks and the possible origins for dust substructures from analysis of the gap and ring properties in Section 4. Finally, the conclusions of this work are summarized in Section 5.

2. Observations and Sample Selection

2.1. Observations and Data Reduction

Our ALMA Cycle 4 program (ID: 2016.1.01164.S; PI: Herczeg) observed 32 disks in the Taurus star-forming region in Band 6 (1.33 mm) with high spatial resolution (∼0 farcs 12, corresponding to ∼16 au for the typical distance to Taurus). Targets were selected for disks around stars with spectral type earlier than M3, excluding duplication at high resolution in archival data, close binaries (0 farcs 1–0 farcs 5), and stars with high extinction (A_V > 3 mag). Further details of the sample will be described in a forthcoming paper.

The 32 disks were split into four different observing groups based on their sky coordinates. All observations were obtained from late August to early September of 2017 using 45–47 12 m antennas on baselines of 21–3697 m (15–2780 kλ), with slight differences in each group (see Table 1). The ALMA correlators were configured identically into four separate basebands for each observation. Two basebands were set up for continuum observations, centered at 218 and 233 GHz with bandwidths of 1.875 GHz. The average observing frequency is 225.5 GHz (wavelength of 1.33 mm). The other two windows cover the two CO isotopologue lines and will not be discussed in this paper. On-source integration times were ∼4 minutes per target for one group with relatively bright disks and ∼10 minutes per target for the other three groups. Table 1 summarizes the details of observation setups in each group.

Table 1. ALMA Cycle 4 Observations

UTC Date	No.	Baseline Range	pwv	Calibrators			On-source Time	Targets
	Antennas	(m)	(mm)	Flux	Bandpass	Phase	(minutes)
2017 Aug 27	47	21–3638	0.5	J0510+1800	J0510+1800	J0512+2927	4	MWC 480
						J0435+2532^a	4	CI Tau, DL Tau, DN Tau, RY Tau
						J0440+2728	4	GO Tau
						J0426+2327	1.5	IQ Tau

2017 Aug 31	45	21–3697	1.3	J1107–4449	J1427–4206	J1058–8003	9–10	CIDA 9, DS Tau

2017 Aug 31–Sep 2	45	21–3697	1.5	J0510+1800	J0423−0120	J0426+2327	8.5	FT Tau, UZ Tau E
						J0435+2532	10	IP Tau

Notes. The 12 disks discussed in this paper come from three observing groups; thus, the observation setup for the remaining group is not shown here.

^aThe scheduled phase calibrator (J0426+2327) for these disks was observed at different spectral windows from the science targets; thus, phase calibration cannot be applied from the phase calibrator to our targets. We used the weaker check source (J0435+2532) instead to transfer phase solutions.

Download table as: ASCII Typeset image

The ALMA data were calibrated using the Common Astronomy Software Applications (CASA) package (McMullin et al. 2007), version 5.1.1. Following the data reduction scripts provided by ALMA, the atmospheric phase noise was first reduced using water vapor radiometer measurements. The standard bandpass, flux, and gain calibrations were then applied (see Table 1). Based on the phase and amplitude variations on calibrators, we estimate an absolute flux calibration uncertainty of ∼10%. Continuum images were then created from the calibrated visibilities with CASA task tclean. For targets with initial signal-to-noise ratio (S/N) ≳ 100 in the image, we applied three rounds of phase (down to the integration time) and one round of amplitude self-calibration. For targets with initial S/N < 100, we applied only one round of phase and one round of amplitude self-calibration. For two disks with S/N < 30, self-calibration was not applied. After each round of self-calibration, we checked the image S/N and would cease the procedure when no significant improvement was measured in the S/N. A few disks had only two rounds of phase and one round of amplitude self-calibration. Self-calibration led to 20%–30% improvements in S/N for most disks and a factor of 2 improvement in S/N for the brightest disks. The data visibilities were extracted from the self-calibrated measurement sets for further modeling. The final continuum images were produced with Briggs weighting and a robust parameter of +0.5 in tclean, resulting in a typical beam size of 0 farcs 14 × 0 farcs 11 and a median continuum rms of 0.05 mJy beam⁻¹. These observations are not sensitive to emission larger than ∼1 farcs 3 (corresponding to ∼180 au for the typical distance of the Taurus region), which is set by the maximum recoverable scale of the chosen antenna configuration.

2.2. Sample Selection

In this paper, we analyze the subsample of disks within our program that show prominent substructures in their dust thermal emission (see dust continuum images and radial profiles for the subsample in Figures 1 and 2). Results for the full sample will be presented in a forthcoming paper.

**Figure 1.** Synthesized images of the 1.33 mm continuum with a Briggs weighting of robust = 0.5. The images are displayed in order of decreasing millimeter flux, from the top left panel to the bottom right panel, and are scaled to highlight the weaker outer emission. The beam for each disk is shown in the left corner of each panel.
Download figure:
Standard image High-resolution image

**Figure 2.** Radial intensity profiles (black lines) along disk major axis for the 12 selected disks with dust substructures, in the same order as in Figure 1. The fitted Gaussian profile is shown in red to highlight the disk substructures, except for CIDA 9 and IP Tau, which have deep inner cavities. The 1σ noise level is shown by the dashed line.
Download figure:
Standard image High-resolution image

Our sample selection of disk substructures is mainly guided by inspection of the disk radial intensity profiles. We first determine the disk major axis by using CASA task imfit to fit an elliptical Gaussian profile to the continuum emission in the image plane. The radial intensity profile along the major axis is then used for an initial classification of disk substructures, including (1) inner cavities, (2) extended emission at large radii, and (3) resolved rings or emission bumps. Twelve of our sample of 32 disks show substructures, with dust emission that cannot be fit with a single smooth central component.

This selection of disks with substructures is confirmed by quantifying the reduced χ² of fits of Gaussian profiles to the radial intensity profile along the disk major axis, within the central 1 farcs 5 of the centroid (refer to red lines in Figure 2). The disks selected for this paper have the largest χ² values. The choice to focus on 12 sources is somewhat arbitrary, but disks with even slightly lower χ² values would include those with subtle deviations from a Gaussian profile that could be well fit with a single tapered power law (see F. Long et al. 2018, in preparation). The source properties for the 12 selected disks are summarized in Table 2.

Table 2. Source Properties and Observation Results

Name	SpTy	T_eff	M_*	log(L_*)	Refs	Distance	rms	Beam Size
		(K)	(M_⊙)	(L_⊙)		(pc)	(mJy beam⁻¹)	(arcsec × arcsec)
CIDA 9	M1.8	3589	0.43	−0.7	HH14	171	0.05	0.13 × 0.10
IP Tau	M0.6	3763	0.52	−0.47	HH14	130	0.047	0.14 × 0.11
RY Tau	F7	6220	2.04	1.09	HH14	128	0.051	0.15 × 0.11
UZ Tau E	M1.9	3574	0.39	−0.46	HH14	131	0.049	0.13 × 0.11
DS Tau	M0.4	3792	0.58	−0.61	HH14	159	0.05	0.14 × 0.10
FT Tau	M2.8	3444	0.34	−0.83	HH14	127	0.047	0.13 × 0.11
MWC 480	A4.5	8460	1.91	1.24	YLiu18	161	0.07	0.17 × 0.11
DN Tau	M0.3	3806	0.52	−0.16	HH14	128	0.05	0.14 × 0.11
GO Tau	M2.3	3516	0.36	−0.67	HH14	144	0.049	0.14 × 0.11
IQ Tau	M1.1	3690	0.50	−0.67	HH14	131	0.076	0.16 × 0.11
DL Tau	K5.5	4277	0.98	−0.19	HH14	159	0.048	0.14 × 0.11
CI Tau	K5.5	4277	0.89	−0.09	HH14	158	0.05	0.13 × 0.11

Note. The distance for each target is adopted from the Gaia DR2 parallax (Gaia Collaboration et al. 2016, 2018). Spectral type and stellar luminosity are adopted from the listed references (Herczeg & Hillenbrand 2014; Liu et al. 2018) and are updated to the new Gaia distance. Stellar masses are recalculated with the stellar luminosity and effective temperature listed here using the same method as in Pascucci et al. (2016). Further details will be described in a forthcoming paper of the full sample.

Download table as: ASCII Typeset image

3. Modeling Disk Substructures

The 1.33 mm continuum images for our 12 disks (in Figure 1) reveal substructures with a wide variety of properties. Resolved rings are the most common type of substructures, characterizing half of our sample. Several disks have two or more rings. Emission bumps are detected from several disks and would likely be resolved into clear rings with higher spatial resolution. Four disks have inner disk cavities (surrounded by one or multiple rings), with different degrees of dust depletion and subtle azimuthal asymmetries. Spirals and high-contrast (with an intensity ratio higher than 2) azimuthal asymmetries are not seen in our sample. These general results are consistent with expectations based on previous results of biased samples, which showed that rings are common while large azimuthal asymmetries (azimuthal dust traps, such as vortices) and spirals are rare.

In this section, we describe the general procedure in modeling the dust substructures performed in the visibility plane and present the results of best-fit models. Disk millimeter fluxes and disk dust sizes are then measured from the best-fit intensity profiles, which will be used in later analysis.

3.1. Modeling Procedure

In order to precisely quantify the observed morphology of dust continuum emission, our analysis is performed in the Fourier plane by comparing the observed visibilities to synthetic visibilities computed from a model intensity profile. Axisymmetry is assumed, since high-contrast asymmetries in the dust emission are not seen (Figure 1; low-contrast asymmetries will be discussed briefly in Section 3.3). Each disk is initially approximated by combining a central Gaussian profile with additional radial Gaussian rings, with the model intensity profile expressed as

$\begin{eqnarray}&&I(r)=A\exp \left(-\displaystyle \frac{{r}^{2}}{2{\sigma }_{0}^{2}}\right)+\displaystyle \sum _{i}{B}_{i}\exp \left[-\displaystyle \frac{{(r-{R}_{i})}^{2}}{2{\sigma }_{i}^{2}}\right],\end{eqnarray} \tag{ 1 }$

where the first term represents the central emission and the second term represents a series of peaks in the radial intensity profile, and R_i and σ_i are the locations and widths of the emission components. In some cases, the central Gaussian profile is replaced with an exponentially tapered power law, which reproduces the oscillation pattern in the visibility profile (Andrews et al. 2012; Hogerheijde et al. 2016; Zhang et al. 2016) and better fits the data (with two more free parameters). The revised model is then described as

$\begin{eqnarray}&&I(r)=A{\left(\displaystyle \frac{r}{{r}_{t}}\right)}^{-\gamma }\exp \left[-{\left(\displaystyle \frac{r}{{r}_{t}}\right)}^{\beta }\right]+\displaystyle \sum _{i}{B}_{i}\exp \left[-\displaystyle \frac{{(r-{R}_{i})}^{2}}{2{\sigma }_{i}^{2}}\right],\end{eqnarray} \tag{ 2 }$

where r_t is the transition radius, γ is the surface brightness gradient index, and β is the exponentially tapered index. The model visibilities are then created by Fourier-transforming the disk model intensity profile using the publicly available code Galario (Tazzari et al. 2018). Fitting the model visibilities to the data visibilities is later performed with the emcee²⁵ package (Foreman-Mackey et al. 2013), in which a Markov chain Monte Carlo (MCMC) method is used to explore the optimal value of free parameters.

Our choice of component type and number in the model intensity profile for each disk is guided by the observed radial profile along the disk major axis (Figure 2). A resolved ring or emission bump is modeled as a Gaussian ring component. The initial guesses for the amplitude, location, and width of each component are also inferred from the radial profiles. The disk inclination angle (i), disk position angle (PA), and position offsets from the phase center (Δα and Δδ) are all free parameters in our fit. The starting point for the four parameters is estimated by fitting an elliptical Gaussian component to the continuum image with CASA task imfit. Prior ranges are set as ±20° for i and PA and ±0 farcs 5 for the position offsets. A uniform prior probability distribution is adopted for each of these parameters.

The radial grid in our model is linearly distributed within 0 farcs 0001–4'' in steps of 0 farcs 001, which is much smaller than our synthesized beam (∼0 farcs 1). We start the MCMC fit by exploring all free parameters (four disk geometric parameters, plus Gaussian profile and Gaussian ring(s)) with 100 walkers and 5000 steps for each walker. The burn-in phase for convergence is typically ∼2000 steps. A second run with parameter ranges confined from the initial run is conducted with another 5000 steps. For the second run, the autocorrelation time is typically 100 steps. The posterior distributions are then sampled using the chains of the last 1000 steps, as well as the optimal value (median value) and its associated uncertainty for each parameter. The statistical uncertainty for each parameter is estimated as the interval from the 16th to the 84th percentile.

In the next step, we perform multiple comparisons between data and model to check the goodness of our best-fit model, including visibility profiles, synthesized images, and radial cuts from the images. If significant symmetrical residuals (≳5–10σ) are present, we either include an additional Gaussian ring component or replace the central Gaussian profile with a tapered power law. These procedures are repeated until a reasonable best fit is found.

3.2. Modeling Results

Detailed results of the best-fit models are presented here, as well as the approach to derive total disk flux and disk size based on the best-fit models. Our final choice of the best-fit model for each disk is guided by using the fewest number of parameters to reproduce the axisymmetric structures with residuals less than ∼5σ. Figure 12 in the Appendix compares the best-fit model with the observed visibility profiles, synthesized images, and radial profiles for each disk. In general, our models fit the disk total flux and disk substructures reasonably well, as indicated from the consistency of data and model visibilities at the shortest baseline and the match of visibility structures at the longer baselines, respectively. For disks with azimuthal asymmetries, our model fails to accurately reproduce the amplitude of the substructure component, but captures the location and width of the ring(s) well, which are the main focus of our analysis below.

3.2.1. Best-fit Models

The best-fit model intensity profiles for the 12 disks are shown in Figure 3, with the substructure component types and numbers of each disk summarized in Table 3. The detailed information (e.g., gap and ring location/width) is provided in Table 4.

Table 3. Disk Model Parameters

Name	F_ν	R_eff	Incl	PA	Δα	Δδ	Morphology/Model Description
	(mJy)	(au)	(deg)	(deg)	(arcsec)	(arcsec)
CIDA 9	${37.1}_{-0.09}^{+0.09}$	${59.0}_{-0.17}^{+0.17}$	${45.56}_{-0.18}^{+0.19}$	${102.65}_{-0.26}^{+0.25}$	−0.51	−0.73	Inner cavity
IP Tau	${14.53}_{-0.08}^{+0.07}$	${34.58}_{-0.26}^{+0.13}$	${45.24}_{-0.33}^{+0.32}$	${173.0}_{-0.42}^{+0.43}$	0.05	0.17	Inner cavity
RY Tau	${210.39}_{-0.1}^{+0.09}$	${60.8}_{-0.13}^{+0.0}$	${65.0}_{-0.02}^{+0.02}$	${23.06}_{-0.02}^{+0.02}$	−0.05	−0.09	Inner cavity + one emission bump
UZ Tau E	${129.52}_{-0.16}^{+0.14}$	${81.61}_{-0.13}^{+0.13}$	${56.15}_{-0.07}^{+0.07}$	${90.39}_{-0.08}^{+0.08}$	0.77	−0.27	Inner cavity + two emission bumps
DS Tau	${22.24}_{-0.11}^{+0.07}$	${67.58}_{-0.32}^{+0.32}$	${65.19}_{-0.13}^{+0.13}$	${159.62}_{-0.14}^{+0.14}$	−0.13	0.22	Inner Gaussian profile + one ring
FT Tau	${89.77}_{-0.1}^{+0.09}$	${42.04}_{-0.13}^{+0.0}$	${35.55}_{-0.16}^{+0.14}$	${121.8}_{-0.27}^{+0.26}$	−0.1	0.13	Inner power-law profile + one emission bump
MWC 480	${267.76}_{-0.21}^{+0.18}$	${104.97}_{-0.16}^{+0.16}$	${36.48}_{-0.05}^{+0.05}$	${147.5}_{-0.08}^{+0.09}$	−0.01	0.0	Inner power-law profile + one ring
DN Tau	${88.61}_{-0.2}^{+0.09}$	${56.06}_{-0.13}^{+0.13}$	${35.18}_{-0.22}^{+0.2}$	${79.19}_{-0.38}^{+0.36}$	0.08	0.0	Inner Gaussian profile + two emission bumps
GO Tau	${54.76}_{-0.2}^{+0.33}$	${144.14}_{-2.45}^{+1.15}$	${53.91}_{-0.2}^{+0.2}$	${20.89}_{-0.24}^{+0.24}$	−0.17	−0.41	Inner power-law profile + two rings
IQ Tau	${64.11}_{-0.34}^{+0.25}$	${95.89}_{-1.05}^{+0.66}$	${62.12}_{-0.2}^{+0.19}$	${42.38}_{-0.23}^{+0.22}$	−0.09	0.07	Inner Gaussian profile + two emission bumps
DL Tau	${170.72}_{-0.16}^{+0.37}$	${147.39}_{-0.16}^{+0.48}$	${44.95}_{-0.09}^{+0.09}$	${52.14}_{-0.14}^{+0.15}$	0.24	−0.06	Inner power-law profile + one emission bump + two rings
CI Tau	${142.4}_{-0.24}^{+0.15}$	${173.8}_{-0.32}^{+0.47}$	${49.99}_{-0.12}^{+0.11}$	${11.22}_{-0.13}^{+0.13}$	0.33	−0.08	Inner Gaussian profile + three emission bumps + one ring

Note. The inclination, PA, and phase center offsets (Δα and Δδ) are parameters fitted with MCMC. Total flux (F_ν) and effective radius (R_eff, with 90% flux encircled) are derived from the best-fit intensity profile for each disk. The quoted uncertainties are the interval from the 16th to the 84th percentile of the model chains. The typical uncertainties for Δα and Δδ are <0 farcs 001, thus not listed. An emission bump or an resolved ring is modeled by a Gaussian ring. The faint outer ring for DL Tau and GO Tau is included to describe the tenuous outer disk, and the faint 3σ ring for MWC 480 is indicated from the fitting residual map. The three faint outer rings are not included in the description column and will not be used in the analysis in Section 4.

Download table as: ASCII Typeset image

Table 4. Gap and Ring Properties

Name	Ring Number	Ring Location	Ring Width	Gap Location	Gap Width	Contrast
		(au)	(au)	(au)	(au)
CIDA 9	1	39.52 ± 0.00	25.31 ± 0.00	⋯	49.28 ± 0.17	⋯
IP Tau	1	27.05 ± 0.07	10.40 ± 0.13	⋯	42.15 ± 0.13	⋯
RY Tau	1	18.19 ± 0.00	25.60 ± 0.13	⋯	13.85 ± 0.13	⋯
	2	49.04 ± 0.14	19.46 ± 0.13	43.41 ± 0.13	4.86 ± 0.20	1.03 ± 0.01
UZ Tau E	1	11.02 ± 0.13	12.05 ± 0.26	⋯	9.46 ± 0.13	⋯
	2	17.42 ± 1.31	28.16 ± 0.66	⋯	⋯	⋯
	3	77.04 ± 0.47	15.72 ± 0.26	69.05 ± 0.20	7.34 ± 0.46	1.07 ± 0.02
DS Tau	1	56.78 ± 0.16	17.17 ± 0.16	32.93 ± 0.32	27.03 ± 0.24	24.07 ± 1.87
FT Tau	1	32.14 ± 0.13	16.51 ± 0.13	24.78 ± 0.19	4.83 ± 0.06	1.37 ± 0.01
MWC 480	1	97.58 ± 0.08	12.56 ± 0.16	73.43 ± 0.16	33.33 ± 0.16	73.78 ± 13.00
DN Tau	1	15.36 ± 0.77	21.12 ± 0.38	⋯	⋯	⋯
	2	53.39 ± 0.95	7.68 ± 0.64	49.29 ± 0.44	3.84 ± 1.21	1.06 ± 0.09
GO Tau	1	73.02 ± 0.16	9.79 ± 0.43	58.91 ± 0.66	22.90 ± 0.86	17.83 ± 5.21
	2	109.45 ± 0.36	22.18 ± 0.86	86.99 ± 0.88	16.13 ± 0.58	4.54 ± 0.63
IQ Tau	1	48.22 ± 1.09	11.79 ± 0.92	41.15 ± 0.63	6.94 ± 1.29	1.10 ± 0.08
	2	82.79 ± 2.88	24.50 ± 1.96	⋯	⋯	⋯
DL Tau	1	46.44 ± 0.48	14.63 ± 0.48	39.29 ± 0.32	6.68 ± 0.48	1.09 ± 0.03
	2	78.08 ± 0.24	8.59 ± 0.64	66.95 ± 0.87	13.83 ± 0.72	6.36 ± 1.32
	3	112.27 ± 0.32	29.57 ± 1.27	88.90 ± 1.11	25.92 ± 0.56	2.11 ± 0.12
CI Tau	1	27.67 ± 0.24	19.28 ± 0.47	13.92 ± 0.32	8.85 ± 0.16	2.20 ± 0.12
	2	61.95 ± 0.47	29.39 ± 1.11	48.36 ± 0.41	10.90 ± 0.40	1.22 ± 0.04
	3	99.22 ± 1.58	8.69 ± 0.95	⋯	⋯	⋯
	4	152.80 ± 0.47	59.72 ± 0.47	118.99 ± 0.65	22.12 ± 0.55	1.67 ± 0.05

Note. Ring and gap properties for each disk. Gap properties and intensity contrasts are blank for the gaps that are not resolved in model intensity profiles.

Download table as: ASCII Typeset image

The inner regions of four disks are described with a Gaussian profile (Equation (1)), four disks are described by a revised power-law model (Equation (2)), and the four other disks lack millimeter emission from their inner disks (see Table 3 for details). For the inner disks described by a power law, the taper index β (>4) corresponds to a sharp outer edge of millimeter-sized dust for the emission of the inner blob, consistent with the prediction of fast radial drift of dust particles (Birnstiel & Andrews 2014).

The four disks with inner cavities are fit with (a sum of) Gaussian ring(s). In three disks (MWC 480, GO Tau, and DL Tau), an additional Gaussian ring in the outermost disk is included to account for the tenuous outer disk edge, which is detected at ∼3σ significance. The inclusion of one more component for GO Tau and DL Tau is needed to avoid generating an outer ring with a width that is much broader than observed.²⁶ The additional ring for MWC 480 at 1'' radius is needed to reproduce the 3σ ring in the residual map, which is found in the fitting of the visibilities, but is too faint to be visible in the observed image. The modeling of the MWC 480 disk by Liu et al. (2018) does not include this component, since they start the modeling in the image plane and focus on reproducing the primary structures. The detailed analysis of the substructure components will be presented in Section 4.

The disk geometry parameters are summarized in Table 3, in which the best-fit inclination and position angles are generally consistent with the values estimated from imfit within 2°–3°. The largest difference of PA is seen in DN Tau, in which our best-fit PA is 6° larger (to the east) than the initial imfit estimation, hinting for some difference in disk orientation between the emission of the inner blob and the outer ring (see also the 3σ residual in Figure 12). The differences between the inclinations and position angles from simplistic models in the image plane with imfit and those from our visibility fitting suggest that the formal errors listed in Table 3 are likely underestimated.

3.2.2. Millimeter Flux and Dust Disk Size

The disk flux densities at 1.33 mm and dust disk sizes are inferred from the model intensity profile, as described in this subsection, and are not model parameters that are directly fit in MCMC. Disk millimeter fluxes and disk dust sizes are summarized in Table 3.

Given an intensity profile, the cumulative distribution could be described as

$\begin{eqnarray}&&{f}_{\nu }(r)=2\pi {\int }_{0}^{r}{I}_{\nu }(r^{\prime} )r^{\prime} {dr}^{\prime} ;\end{eqnarray} \tag{ 3 }$

thus, the total flux is ${F}_{\nu }={f}_{\nu }(\infty )$ by definition. The millimeter flux for each disk is measured by integrating over the best-fit intensity profile. We then randomly choose 100 models in the last 1000 steps (×100 walkers) of our MCMC chain to estimate flux uncertainty as the central interval from 16th to 84th percentile. For most disks, our flux measurements at 1.33 mm are consistent with pre-ALMA interferometry measurements²⁷ within uncertainties (Andrews et al. 2013), assuming 10% and 15% absolute flux uncertainty for ALMA and pre-ALMA results, respectively. Our flux densities for CI Tau, FT Tau, and IP Tau are more than 30% brighter than those reported in Andrews et al. (2013). However, the measured flux density for CI Tau is highly consistent with a recent ALMA measurement (Konishi et al. 2018), and the FT Tau flux density is similar to a past CARMA measurement (Kwon et al. 2015). For IP Tau, the flux difference is reconciled if the SMA measurement at 0.88 mm is extrapolated to 1.33 mm with a spectral index of 2.4 (Andrews et al. 2013; Tripathi et al. 2017). These modest inconsistencies are likely related to unknown systematic flux calibration uncertainty, self-calibration, and different methods in estimating fluxes. These differences in fluxes will not affect the results in our following analysis.

The effective disk radius, R_eff, is defined here as the radius where 90% of the total flux is encircled (see, e.g., Tripathi et al. 2017). The uncertainty for disk size is estimated in the same way as the flux uncertainty. We do not compare the disk sizes with results in Tripathi et al. (2017) for the few overlapping disks, since the two works probe different wavelengths and use different size metrics.

3.3. Residuals and Azimuthal Asymmetries

The best fits to the observed visibilities yield significant residuals (>10σ) for a few disks. These residuals indicate azimuthal asymmetries for the innermost rings of CIDA 9, RY Tau, and UZ Tau E. One characteristic feature of this set of disks is that their inner regions are depleted of dust (including marginal depletions). High-contrast asymmetries have been observed in some transition disks (e.g., IRS 48 by van der Marel et al. 2013) and interpreted as vortices that could be triggered by the presence of planets. An eccentric cavity, induced by companions in the inner disk, could also create azimuthal asymmetries, with contrast levels depending on the mass of the companion (Ataiee et al. 2013; Ragusa et al. 2017). The azimuthal asymmetry of the AA Tau disk has alternatively been attributed to a misalignment between the inner and outer disks (Loomis et al. 2017). Additional azimuthal structures on top of an underlying axisymmetric disk model would be needed to better describe the emission pattern, which is beyond the scope of this paper.

The inner emission blobs of CI Tau and GO Tau also return modest residuals of ∼5σ. These inner emission regions have a narrow extent of 0 farcs 1–0 farcs 2 in radius, so subtle radial variations might be present but are not well enough resolved to interpret here. The hot-Jupiter candidate around CI Tau found by Johns-Krull et al. (2016) is in a 9-day orbit and likely does not affect the rings detected here on much larger radii.

4. Results and Discussion

Previous measurements of disk substructures have been biased to brighter disks or disks in which the presence of substructures has already been inferred from other observations. The disk substructures identified in this survey are seen for the first time²⁸ at ∼0 farcs 1 resolution in an unbiased study that covers a wide range in fluxes within a given range of stellar mass.

From our full sample of 32 disks, we have identified 12 disks with substructures in their dust continuum emission. Four disks have inner cavities in the millimeter continuum, encircled by single rings for CIDA 9 and IP Tau and multiple rings for RY Tau and UZ Tau E. Three disks (FT Tau, DS Tau, and MWC 480), have millimeter continuum emission characterized by an inner disk encircled by a single ring. Five disks (CI Tau, DL Tau, GO Tau, IQ Tau, and DN Tau) have an inner disk encircled by multiple rings. The location and shape for each of these components are modeled as symmetric Gaussian profiles and are fit in the visibility plane (Section 3.1).

Table 4 summarizes the results from our fits, including the size of inner cavities, the radial location and width of gaps and rings, and the flux contrast ratio between the rings and gaps. The properties of the substructures are disparate, with radial locations from 10 to 120 au, rings with emission that accounts for 10%–100% of the total flux from the disk, and widths that are usually ∼0.2 times the radial location of the gap but can be wider. The presence of most of these substructures does not obviously depend on any disk or stellar property.

In this section, we synthesize these disparate properties in an attempt to identify the physical mechanism(s) that produce cavities and rings. We begin by exploring the parameter space occupied by our sample to describe the star and disk properties of our substructures. We then apply our results to expectations for the properties of gaps and rings that could be introduced by condensation fronts and by planets. The bulk of gaps could be carved by planets with masses close to the minimum planet mass able to produce gas pressure bumps, while less than half of the gaps are close to volatile condensation fronts.

4.1. Source Properties for Disks with Substructures

Substructures in our sample are present in objects that cover a wide range in stellar and disk mass. The left panel of Figure 4 shows the location of our 12 disks with dust substructures in the M_*–M_dust plane, as well as the other 20 disks in our full sample, which do not show dust substructures at our current resolution. The full Taurus sample from Andrews et al. (2013) is also included in this plot to provide a broader comparison. The dust masses are estimated from the 1.33 mm continuum flux density (e.g., Beckwith et al. 1990) by

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

with a dust opacity ${\kappa }_{\nu }=2.3\,{\mathrm{cm}}^{2}\,{{\rm{g}}}^{-1}\times {(\nu /230\mathrm{GHz})}^{0.4}$ and a Planck function B_ν(T_dust) for a dust temperature of 20 K for each disk at distance D, assuming that the dust is optically thin. Stellar effective temperatures and stellar luminosities are adopted from Herczeg & Hillenbrand (2014) and then updated for individual Gaia DR2 distances (Gaia Collaboration et al. 2018) for the full sample. The stellar masses are then calculated from the Baraffe et al. (2015) and nonmagnetic Feiden (2016) evolutionary tracks, following Pascucci et al. (2016). The disk dust masses are calculated with updated Gaia DR2 distance for individual objects with millimeter fluxes adopted from our measurements for the 12 disks and adopted from Andrews et al. (2013) for the other Taurus members. The uncertainties of our estimated dust masses only consider the uncertainties of flux measurements and do not take into account the differences in dust temperature and dust optical depth.

**Figure 4.** Left: stellar mass vs. disk mass for the 12 disks with substructures (in blue open circles for the four disks with inner cavities), the 20 disks without substructures in current observations (in orange, using disk masses from Andrews et al. 2013), and the full Taurus sample (in gray, upper limits in triangles) of Andrews et al. (2013). The relationship between stellar mass and disk mass for transition disks (dotted blue line) is taken from Pinilla et al. (2018), with the shaded region showing the typical data scatter. The typical errors in $\mathrm{log}$ (M_dust) of ∼0.04 dex, including the 10% flux calibration uncertainty, and in $\mathrm{log}$ (M_*) of 0.1 dex are shown in the left corner. Right: disk effective radius vs. disk dust mass for the 12 disks with substructures, with colors and symbol shapes separating disks with single, double, and multiple rings and disks with inner cavities.
Download figure:
Standard image High-resolution image

**Figure 4.** Left: stellar mass vs. disk mass for the 12 disks with substructures (in blue open circles for the four disks with inner cavities), the 20 disks without substructures in current observations (in orange, using disk masses from Andrews et al. 2013), and the full Taurus sample (in gray, upper limits in triangles) of Andrews et al. (2013). The relationship between stellar mass and disk mass for transition disks (dotted blue line) is taken from Pinilla et al. (2018), with the shaded region showing the typical data scatter. The typical errors in $\mathrm{log}$ (M_dust) of ∼0.04 dex, including the 10% flux calibration uncertainty, and in $\mathrm{log}$ (M_*) of 0.1 dex are shown in the left corner. Right: disk effective radius vs. disk dust mass for the 12 disks with substructures, with colors and symbol shapes separating disks with single, double, and multiple rings and disks with inner cavities.
Download figure:
Standard image High-resolution image

Our sample focuses on disks around M dwarfs and solar-mass stars in Taurus, with a requirement that the spectral type of the star is earlier than M3 (corresponding to ∼0.25 M_⊙ in the Baraffe et al. 2015 evolutionary tracks and ∼0.45 M_⊙ in the magnetic Feiden 2016 tracks, for an age of 2 Myr). Disks with dust substructures cover this full stellar mass range of our whole sample. The two disks around early spectral types (the A4 star MWC 480 and the F7 star RY Tau) both have prominent dust rings. Disk dust masses for our 12 disks scatter over more than one order of magnitude, even in a narrow stellar mass bin. Dust substructures seem to be more common in brighter disks. A more complete analysis of the statistics with respect to the parent sample of 32 objects will be presented in a forthcoming paper.

Four disks in our sample have resolved inner cavities (including marginal depletion), including three new discoveries and confirmation of the inner cavity of RY Tau found by Pinilla et al. (2018). None of these four disks show any signature of a cavity based on the SED (see also Figure 13 in Andrews et al. 2013), so they all have warm dust near the star, similar to some of the cavities found by Andrews et al. (2011). In an analysis of 29 disks with inner millimeter cavities observed by ALMA, Pinilla et al. (2018) found a flatter M_*–M_dust relation when compared to the correlation obtained for all disks from several different star-forming regions. Inner cavities may therefore be more common among more massive disks, regardless of stellar mass. Three of our four inner cavity disks are consistent with this correlation and are in the upper end of masses for all Taurus disks. The exception, IP Tau, has the smallest and faintest disk in our sample. Some outliers may be expected from this relationship for circumbinary disks, such as the disk around CoKu Tau 4 (D'Alessio et al. 2005; Ireland & Kraus 2008). With a cavity radius of 21 au, a companion to IP Tau would need to be located at ∼10 au, or ∼0 farcs 06 (Artymowicz & Lubow 1994). Previous binary searches would have been unable to resolve such a close companion. Unlike CoKu Tau 4, IP Tau must have an inner disk to explain the IR excess (Furlan et al. 2006) and active, though very weak, accretion (Gullbring et al. 1998). Future high angular resolution and high-sensitivity observations of faint disks will test whether the Pinilla et al. (2018) relationship is robust to the selection bias that past ALMA observations (Cycle 0 to Cycle 3) had toward brighter disks.

For the four disks with inner cavities, only one ring is detected from IP Tau and CIDA 9, while two rings are detected from RY Tau and three from UZ Tau E. The outer substructures of RY Tau and UZ Tau E are not clearly seen in the images but are detected in the uv-plane and in the cross-cut of the image along the major axis in Figure 2. In contrast, the disks of IP Tau and CIDA 9 have large inner cavities, with faint dust emission detected at low S/N (see the right panel of Figure 4) that suggests larger depletion of millimeter grains. However, the outer substructures of RY Tau and UZ Tau E would still have been detected if the S/N were scaled down to match the fainter signal of IP Tau and CIDA 9. One possibility is that narrow cavities, as in the case of RY Tau and UZ Tau E,²⁹ will evolve into larger and more depleted cavities. Nevertheless, more observations of inner cavity disks spanning different cavity sizes, ages, and dust depletion factors are required to test this hypothesis.

As shown in the right panel of Figure 4, the number of substructures (single, double, or multiple rings) seems independent of whether millimeter continuum emission is present in the inner disk. Disks with similar brightness have dust distributions that are diverse in shape and in size. However, disks with multiple rings tend to be brighter and more extended, implying that substructures that originate at larger radii may act as mechanisms (e.g., dust traps) preventing the loss of millimeter-sized dust due to inward drift, thus retaining both the high dust mass and large disk size (e.g., Pinilla et al. 2012b). This explanation should be valid if R_gas/R_dust is higher in smaller disks, for which dust inward migration is very efficient without "traps" formed at larger radii, though this is not seen in Ansdell et al. (2018). Alternatively, the initial distribution of disk sizes in young stellar objects may be bimodal, with some large and some small, depending on the alignment of the rotation vector and magnetic field (Tsukamoto et al. 2015; Wurster et al. 2016)—although this scenario may be complicated by initial angular momentum distribution, magnetohydrodynamic structure, and turbulence (see review by Li et al. 2014). Disk sizes for the other 20 disks in the full sample that do not show dust substructures are generally more compact. A detailed comparison of disk sizes between the two subsamples will be presented in our survey overview paper.

4.2. Properties of Gaps and Rings

In this section, we explore the general properties of gaps and rings revealed from our observations and their implications for disk properties and evolution. We analyze a total of 19 gap and ring pairs (e.g., a gap and the associated ring emission exterior to the gap) from our 12 disks. The very faint, outermost components of the MWC 480, GO Tau, and DL Tau disks are excluded in the analysis, since they were added to the fit to characterize the extended tenuous disk outer edge and do not necessarily represent a physical ring. Four ring components that are not well resolved (the gap interior to the ring peak is not present), as seen in the model intensity profiles (see also Table 4), are excluded from this analysis. The inner cavities are also not considered in some of this analysis, since the gap locations cannot be well determined, although the gaps exterior to the inner ring are included.

Each substructure is described by the gap location, the gap width, and the intensity contrast ratio, as measured in our model fits. The gap location is defined here as the radius where the intensity profile reaches a local minimum interior to the ring. The gap width is defined as the full width at half depth, in which the depth is the difference between the intensity at the gap location and the peak value of its outer ring. We also measure the ring–gap contrast ratio as the intensity ratio at ring peak and gap location. The uncertainty for each parameter is estimated from the 16th–84th percentile (1σ) values from the chains of the last 1000 steps of our MCMC calculations.

As shown in Figure 5, gaps are located from 10 to 120 au with no preferred distance. Most gaps are narrow and unresolved.³⁰ With higher resolution, more narrow and lower-contrast substructures would be expected to emerge (e.g., for TW Hydra, more rings were revealed in observations with 0 farcs 02 resolution observations (Andrews et al. 2016) than with 0 farcs 3 resolution (Zhang et al. 2016)). Narrow gaps around 100 au are absent in the current observations. A weak trend might be seen between gap location and gap width, in which gaps located farther out have larger width, broadly consistent with the case of planet–disk interaction as will be discussed in Section 4.3.1. Moreover, as shown in Figure 6, gap location does not depend on disk mass. We might see a desert of gaps located outside 40 au for less massive disks, which seems to be consistent with the L_mm–R_eff relation that fainter disks tend to be smaller in sizes (Tazzari et al. 2017; Tripathi et al. 2017). Alternatively, in order to retain a massive disk, substructures should be formed at larger radii, or the outer disk will be drained through fast inward drift.

farcs — **Figure 5.** Gap location vs. gap width for the 19 well-resolved gap and ring pairs (including the four inner cavities, with location set to 0 and shown as open circles; see the text for more details of which ring components are excluded). The gray dashed lines represent the typical beam size (012), adopting a typical distance of 140 pc.
Download figure:
Standard image High-resolution image

**Figure 6.** Gap location as a function of disk dust mass for 15 gaps (excluding the four inner cavities).
Download figure:
Standard image High-resolution image

Most ring–gap pairs have intensity contrast ratios lower than 3, with a few very depleted exceptions (the ring of DS Tau and MWC 480, the first ring of GO Tau) with ratios exceeding 20. Of the four inner cavities, IP Tau and CIDA 9 have nearly empty inner holes, while RY Tau and UZ Tau E only have a factor of two depletion. Figure 7 illustrates the relative flux in each ring with respect to the total disk flux, which peaks around 0.2, with a tail toward higher fraction. Except for the two rings around inner cavities for IP Tau and CIDA 9, the two rings in DS Tau and FT Tau (two single-ring disks) have more than 60% of the total disk flux. The rings in disks with multiple substructures (e.g., CI Tau, DL Tau, GO Tau) generally hold ∼20% of the total disk flux. This quantity is approximately proportional to the fraction of dust mass within the ring, though the optical depth of the dust and the temperature differences between the ring and the rest of the disk lead to substantial uncertainties. Since the back-reaction of dust on gas is strongest when dust-to-gas density ratio is of the order of unity (Youdin & Goodman 2005), a typical 20% accumulation of dust in the rings suggests that the creation of an individual ring in general may not be so relevant for the global disk dynamical evolution. However, the total effect of all rings, including any rings in the inner disk that we could not detect and rings emerging from single-ring systems, may affect the dynamical evolution of the disk.

**Figure 7.** Left: fractional flux of the disk in each ring. The four rings associated with inner cavities are shown as open circles. For disks with multiple rings, the adjacent numbers indicate the relative position of the ring component from the star. Right: histogram of fractional flux in the ring for 19 rings.
Download figure:
Standard image High-resolution image

4.3. Possible Origins for Gaps and Rings

The exciting discovery of gap- and ring-like features in the young HL Tau system (ALMA Partnership et al. 2015) suggests that dust particles get trapped in local gas pressure bumps. This and subsequent observations, including those presented here, have revealed that rings are prevalent in protoplanetary disks, with an importance that has motivated the development of many theoretical explanations of the observed substructures. Pressure bumps could be created outside the orbit of a planet (e.g., Pinilla et al. 2012b) or the outer edge of a low-ionization region (the so-called dead zone; Flock et al. 2015). Other magnetohydrodynamic effects, including zonal flows (e.g., Johansen et al. 2009), could also play a role in gas evolution and gas pressure distribution. Our focus of this section is to compare the rings and gaps in our sample to two popular hypotheses, that rings are carved by embedded planets or induced by condensation fronts, followed by future perspectives in discerning different mechanisms at play.

4.3.1. Planet–Disk Interactions

One of the widely invoked explanations for the observed gaps and rings in protoplanetary disks is related to the presence of embedded planets orbiting around the central star. The mass of the planet, the viscosity, and pressure forces of the disk combine to determine the dynamical evolution of disk–planet interactions and thus the resulting distribution of the millimeter-sized dust grains. For the purposes of this subsection, we assume that the gaps are carved by planets and then use the ring and gap properties to estimate the masses of the hidden planets. We then compare these results to statistics from exoplanet observations.

For planet–disk interactions, the gap location should occur at the orbital radius of the planet. The ring of millimeter-size dust grows at the location of the local pressure maximum in the gas, outside the planet orbit. A more massive planet will build a steeper pressure gradient, thereby forming a deeper and wider gap than would be created by a less massive planet (Fung et al. 2014; Kanagawa et al. 2015; Rosotti et al. 2016).

The minimum mass of a planet that could form a gap in the gas density structure, leading to a local pressure bump beyond the planet's orbit, may be described analytically by

$\begin{eqnarray}&&\displaystyle \frac{{M}_{{\rm{p}}}}{{M}_{\star }}\propto {\left(\displaystyle \frac{H}{{r}_{{\rm{p}}}}\right)}^{a}{\alpha }^{b},\end{eqnarray} \tag{ 5 }$

where r_p is the distance of the planet to the star and α is the turbulence parameter, with power-law indices a ∈ [2, 3] and b ∈ [0, 1] (see derivations in Lin & Papaloizou 1993; Duffell & MacFadyen 2012, 2013; Ataiee et al. 2018). If the disk is in vertically hydrostatic equilibrium (H = c_s/Ω_k), assuming uniform α and a power-law profile for the temperature (T ∝ r^−1/2), the minimum planet mass able to create a pressure bump scales as ${M}_{\star }^{c}{r}_{{\rm{p}}}^{d}$ with c ∈ [−1/2, 0] and d ∈ [1/2, 3/4]. Assuming that the gap–ring distance (the distance of gap minimum and the ring peak, an alternative measurement for the gap width) scales with the Hill radius of the planet ( ${R}_{\mathrm{Hill}}={r}_{{\rm{p}}}{({M}_{{\rm{p}}}/(3{M}_{* }))}^{1/3}$ ), and taking the planet masses given by Equation (5), the gap width normalized to the gap location is expected to scale as ${M}_{\star }^{e}{r}_{{\rm{p}}}^{f}$ with e ∈ [−1/3, −1/2] and f ∈ [1/6, 1/4], i.e., a weak dependence on both parameters. Figure 8 shows the distance between the ring peak and gap center, normalized to the gap location (presumably the location of any potential planet). The normalized gap–ring distance is typically 0.2–0.3, with only two gaps as high diagnostic outliers (DS Tau and the closest gap of CI Tau). Given the lack of a clear trend between planet mass indicator and planet location, the mass of most planets in our sample (except for the outliers of CI Tau and DS Tau) might be close to the minimum planet mass able to produce a pressure bump beyond the planet orbit.

**Figure 8.** Left: gap-to-ring distance normalized to gap location (an indicator of planet mass) as a function of gap location. Right: histogram of the gap-to-ring distance.
Download figure:
Standard image High-resolution image

If we simply assume that the gap radius corresponds to 4R_Hill (Dodson-Robinson & Salyk 2011), most of our gaps are related to planets with mass of 0.1–0.5 M_J, as shown in Figure 9. These estimated planet masses³¹ have large uncertainties and should be interpreted as upper limits, since gap radius could extend to (7–10)R_Hill (e.g., Pinilla et al. 2012b). An alternative way to estimate the mass of a planet associated with the gap is by linking the diagnostic of the gap–ring distance to hydrodynamic simulations (Rosotti et al. 2016). The planet mass derived from this diagnostic highly depends on disk viscosity. When the turbulence parameter (Shakura & Sunyaev 1973) is assumed to be α = 10⁻⁴, a value consistent with recent turbulence constraints by Flaherty et al. (2015) and Flaherty et al. (2018), the diagnostic of 0.25 corresponds to a planet mass of ∼15 M_⊕ (0.05 M_J). For a viscous accretion disk (α = 10⁻²; Hartmann et al. 1998), the related planet mass is about 0.3 M_J (Rosotti et al. 2016; S. Facchini et al. 2018, in preparation). The closest gap of CI Tau has the largest normalized gap–ring distance, corresponding to a Jupiter-mass planet in all cases. The large, spatially resolved gap in DS Tau provides us with a hint that a massive planet may also reside in this disk. These simulations were performed for a central star with 1 M_⊙ (Rosotti et al. 2016; S. Facchini et al. 2018, in preparation); thus, the inferred planet mass should be rescaled to the same planet/stellar mass ratio for stars with different masses. The estimated planet masses are subject to large uncertainties in disk properties defined in the models, including disk temperature, viscosity, and scale height. This picture is also complicated by the ability of super-Earths to open multiple gaps in inviscid disks (Bae et al. 2017; Dong et al. 2017, 2018). Dynamical interactions between a gaseous disk and an embedded planet could also produce spiral arms, which are not detected in our sample, probably because of the weak coupling of millimeter grains with the gas.

**Figure 9.** Gap-to-ring distance normalized to gap location (an indicator of planet mass) as a function of stellar mass. The dashed line represents how the planet mass indicator scales with stellar mass for a given planet mass, assuming that the gap–ring distance is 4R_Hill.
Download figure:
Standard image High-resolution image

The analysis above excludes the four inner cavities. The infrared SEDs of these four systems indicate that small dust grains are still present in the inner disks, despite the depletion of millimeter-sized grains that is observed. Low-mass planet(s) with traffic jam effects could be responsible for the cavity opening (Rosotti et al. 2016).

The planets inferred from these analyses are challenging to compare to the statistics of known exoplanets. Very few stars have massive planets at large radii, and even fewer planets are found around M dwarfs (see review by Bowler 2016); however, current sensitivities prevent the detection of the lower-mass planets that would create the gaps located at 20–100 au identified here. Systems with super-Earths/sub-Jovian mass planets at radii larger than 30 au would be unlike our own solar system (ice giants of ∼20 M_⊕ at 20–30 au), but may be prevalent. The most common types of planets found in Kepler transit and microlensing are super-Earths and Neptunes (Winn & Fabrycky 2015; Suzuki et al. 2016; Pascucci et al. 2018), although the distribution at larger radii is not well constrained (Clanton & Gaudi 2016; Meyer et al. 2018).

4.3.2. Condensation Fronts

As disk temperature decreases toward larger radii, a series of major volatile molecules freeze out onto dust grains. These phase transition regions, also referred to as condensation fronts or ice lines, are locations where dust opacities and collisional growth/fragmentation are expected to change, producing features like gaps and rings seen in dust images (e.g., HL Tau; Banzatti et al. 2015; Zhang et al. 2015; Okuzumi et al. 2016). If correct, this explanation should be universal and apply to all disks.

For an irradiated flared disk, Kenyon & Hartmann (1987) parameterize the disk midplane temperature as

$\begin{eqnarray}&&T(r)={T}_{\star }{\left(\displaystyle \frac{{R}_{\star }}{r}\right)}^{1/2}{\phi }_{\mathrm{inc}}^{1/4},\end{eqnarray} \tag{ 6 }$

where the flaring angle ϕ_inc is assumed to be 0.05 (Dullemond & Dominik 2004). The temperature profile for each disk is shown in Figure 10. Taking the stellar luminosity³² $({L}_{\star }=4\pi {R}_{\star }^{2}\sigma {{T}_{\star }}^{4})$ into account, the ice line location for a specific species scales as $r\propto {{L}_{\star }}^{1/2}{T}_{\mathrm{conden}}^{-2}$ . If gaps are formed around ice lines, the gap locations should occur in regions defined by major volatiles. Since the ice line regions associated with the abundant H₂O and NH₃ ice are unresolved for most of our observations (Figure 10), our focus here is on molecules with lower condensation temperature: N₂ (12–15 K), CO (23–28 K), and clathrate-hydrated CO and N₂ (41–46 K). The condensation temperatures, adopted from Zhang et al. (2015), correspond to gas number densities of 10¹⁰–10¹³ cm⁻³ in the disk midplane.

Figure 10 compares the location of gaps to the expected location of ice lines for each target. We also summarize the comparison in Figure 11 in the $r-\sqrt{{L}_{* }}$ plane. Five gaps lie close to the CO ice line (i.e., blue-shaded region matched with gap location), the closer-in gap of CI Tau and the gap of MWC 480 are located around the clathrate-hydrated CO+N₂ ice line, and the outer gap of GO Tau is located at a region with a temperature consistent with the N₂ ice line. The other seven gaps are unrelated to any of the ice lines we consider here. The four inner cavities are not included in this discussion, since temperature in the inner cavity could match with condensation temperatures of a series of major volatiles. For the disks with multiple, well-resolved gaps (CI Tau, DL Tau, and GO Tau), one to two gaps in each disk correspond to an ice line location, but one to two rings in each also do not match an ice line. Among the gaps that are consistent with condensation fronts, explaining the wide and deep gaps of DS Tau and MWC 480 using condensation fronts alone would be challenging (see simulated images around ice lines in Pinilla et al. 2017).

These comparisons between gap locations and snow lines depend on the disk temperature profile (Pinte & Laibe 2014) and the condensation temperatures for different species, both of which suffer from significant uncertainties. A detailed radiative transfer modeling on MWC 480 by Liu et al. (2018) yields a midplane temperature for millimeter-sized grains that is a factor of 1.5–1.9 cooler than the parameterized temperature from Equation (6). A simple experiment of radiative transfer modeling³³ on the other 11 disks yields temperature profiles that are similar to analytic solutions (Equation (6)) beyond 20 au, where our gaps are located, for all sources except for the brighter objects, RY Tau. Therefore, the comparison of ice line locations to gap locations is not significantly affected.

The temperature range for volatile condensation also varies under different conditions and depends on the dust grain size, composition, and surface area that the molecule freezes onto. Taking CO ice as an example, Öberg et al. (2011) suggest an average condensation temperature of 20 K for CO ice, lower than the value adopted from Zhang et al. (2015), due to a different assumption for the gas number density in the disk midplane. The disk region corresponding to the CO ice line would then move outward. The specific gaps that are associated with CO ice lines would change, but their total number would not increase.

The behavior of dust particles around ice lines is currently not well understood. From a physical point of view, inward drifting particles lose their surface ice when crossing the condensation fronts, causing a higher dust-to-gas ratio just outside the ice line, and the evaporated gas may diffuse outward and recondense onto dust outside the ice line, both leading to enhanced grain growth beyond the ice line (Cuzzi & Zahnle 2004; Stammler et al. 2017). Differential grain growth inside and outside the ice lines, producing spectral index variations, has been suggested around the H₂O ice line (Banzatti et al. 2015), with observational evidence from the outburst system V883 Ori (Cieza et al. 2016). Any change in spectral index of dust emission has not yet been detected for the CO ice line (Stammler et al. 2017). Suppressed grain growth for dust particles outside the ice lines due to sintering could also cause pileups of dust in the region where smaller grains have lower drift velocity (Okuzumi et al. 2016). In addition, Zhang et al. (2015) suggested rapid grain growth at radii corresponding to different ice condensation fronts, reducing the detectable millimeter flux in the gaps. The nature of ice sublimation and condensation near these ice lines is complex, because of the strong dependence on the evolution of dust aggregates in terms of fragmentation/coagulation for dust with different compositions, as well as radial drift and turbulent mixing overall (Pinilla et al. 2017).

Given the physical and chemical complexity of the problem, reflected in the model uncertainties described above, and the absence of a clear correlation and correspondence between the gap locations and ice lines, it is still very challenging to test and reconcile the observed substructures with the locations of different ice lines. Volatile condensation fronts of CO, N₂, and clathrate-hydrated CO+N₂ may not be a universal solution for all observed gaps and rings for our sample, but could play a role in shaping some disk dust structures.

4.3.3. Distinguishing Different Mechanisms of Gap Creation

Distinguishing between the different mechanisms responsible for gap and ring formation will yield a better understanding of planet formation, either by providing an indirect probe of planets or by providing a diagnostic of physics that would likely be important in the growth of planetesimals. The rings and gaps in our sample have a wide range of properties. The radial locations of the rings do not seem to prefer the expected location of snow lines, and the gaps of DS Tau and MWC 480 are likely too wide to be explained by ice lines. Around ice lines, the composition of particles is altered and hence their aerodynamical behavior with the gas (Pinilla et al. 2017), but without significant changes of the gas density profile, and therefore pressure bumps are not created.

On the other hand, the variety of ring properties might be expected for hidden exoplanets, which would occur at a wide range of locations and masses. Unfortunately, the planet mass function and their radial distributions are not yet well enough known at such low masses to be able to test whether Neptune-mass planets are prevalent at these distances.

Comparing dust morphology for different grain sizes observed at different wavelengths will allow us to distinguish particle trapping from the other mechanisms (e.g., ice lines). Additional tests of these and other theories will emerge with better maps, including the distribution of gas in the disk. These observations require deep integrations but are feasible for small numbers of disks as follow-ups from surveys such as ours. Indeed, CO cavities and gaps have been reported in a few bright disks (e.g., Isella et al. 2016; Fedele et al. 2017; Boehler et al. 2018), lending some evidence for planet–disk interaction. However, some degeneracy is introduced because the CO gaps could be interpreted as either a gap in the gas density or one in the thermal structure (Facchini et al. 2017).

5. Conclusions

This paper presents the analysis of dust substructures detected at 1.33 mm continuum emission from 12 disks, observed at ∼0 farcs 12 resolution in an ALMA Cycle 4 survey of 32 disks in the Taurus star-forming region. Rings and gaps are the most common type of substructure in our selected disks and exhibit a wide variety of properties. Disk model fitting is performed in the visibility plane to quantify the amplitude, location, and width for each substructure component. We then study the stellar and disk properties for the selected sample in the context of dust morphology and the origins of dust substructures from analysis of the gap and ring properties. Our main findings are summarized as follows:

1.
The 12 disks with detected substructures span a wide range in stellar mass and disk brightness. Disks with multiple rings tend to be more massive and more extended (larger effective radius) than those with single rings.
2.
Four disks are identified with inner dust cavities with a radius of 5–25 au and different levels of depletion. The IR SEDs reveal the presence of small dust grains where large grains are depleted, consistent with expectations for dust filtration by a low-mass planet. Three of these four disks are relatively massive, consistent with expectations from the M_*–M_dust relationship for disks with cavities. These disks may be a collection of heterogeneous sources; more observations toward the fainter end might help to understand the origins of disk inner cavities.
3.
We resolve 19 gap–ring pairs in the model intensity profiles of the 12 disks. The locations, sizes, and contrasts of these gap–ring pairs are used to investigate how these features form. Dust gaps are located from 10 to 120 au. Most gaps have narrow widths that are smaller than the beam size (∼16 au). A typical fraction of 20%–30% of millimeter fluxes are accumulated in each ring, with a few exceptions.
4.
The presence of wider gaps at larger radii hints for planet–disk interaction. The low intensity contrast in most ring and gap pairs suggests the possible link to low-mass planets. We follow the diagnostic used in planet–disk interaction simulations (the separation of ring and gap normalized to gap location) to infer planet mass and find that super-Earths and Neptunes are good candidates if disk turbulence is low (α = 10⁻⁴), in line with the most common type of planets discovered so far.
5.
We do not observe a concentration of gap radii around major ice line locations. While five gaps are located close to the expected radii of CO ice lines, and another one and two gaps are related to N₂ and clathrate-hydrated CO+N₂ ice lines, several other gaps have no relationship with the estimated radii of any ice line we considered. Forming gaps and rings around condensation fronts may not be a universal explanation for all our observed substructures, but could play some role in shaping some of the disks. If ice lines cause rings, it remains unclear why condensation fronts would affect only some of the disks.

Multiwavelength observations that probe different sizes of dust grains will help to discern particle trapping (e.g., planets) from other mechanisms (e.g., ice lines), with the expected narrower dust accumulation at longer wavelength if dust particles are trapped in rings. Characterizing the gas distribution across the gaps will also be essential to better constrain their origins. Follow-up observations in these forms are timely and will accelerate our understanding of disk evolution and planet formation.

We thank the anonymous referee for useful comments. F.L. thanks Sean Andrews for discussions on comparing ALMA fluxes to past interferometer fluxes, and Zhaohuan Zhu, Kaitlin Kratter, and Andrew Youdin for insightful discussions. F.L. and G.J.H. thank Subo Dong for discussions of planet mass functions and radial distributions. F.L. and G.J.H. are supported by general grants 11473005 and 11773002 awarded by the National Science Foundation of China. P.P. acknowledges support by NASA through Hubble Fellowship grant HST-HF2-51380.001-A awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS 5-26555. G.D. acknowledges financial support from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant agreement no. 681601). D.H. is supported by European Union A-ERC grant 291141 CHEMPLAN, NWO, and by a KNAW professor prize awarded to E. van Dishoeck. C.F.M. acknowledges support through the ESO Fellowship. M.T. has been supported by the DISCSIM project, grant agreement 341137 funded by the European Research Council under ERC-2013-ADG. F.M., G.v.d.P., and Y.B. acknowledge funding from ANR of France under contract no. ANR-16-CE31-0013 (Planet-Forming-Disks). D.J. is supported by NRC Canada and by an NSERC Discovery Grant. B.N. and G.L. thank the support by the project PRIN-INAF 2016 The Cradle of Life—GENESIS-SKA (General Conditions in Early Planetary Systems for the rise of life with SKA). Y.L. acknowledges support by the Natural Science Foundation of Jiangsu Province of China (grant no. BK20181513) and by the Natural Science Foundation of China (grant no. 11503087). This paper makes use of the following ALMA data: 2016.1.01164.S. ALMA is a partnership of ESO (representing its member states), NSF (USA), and NINS (Japan), together with NRC (Canada), MOST and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO, and NAOJ. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement.

Software: Galario (Tazzari et al. 2018), emcee (Foreman-Mackey et al. 2013), RADMC-3D (Dullemond et al. 2012), CASA (v5.1.1; McMullin et al. 2007).

Appendix: Fitting Results for Individual Disk

The comparison of data and best-fit model for each disk is shown in Figure 12, in which we check the goodness of our fit through visibility profile, synthesized image, and radial cut. In most cases, the maximum residual in the image is ∼5σ.

Gaps and Rings in an ALMA Survey of Disks in the Taurus Star-forming Region

Article metrics

Permissions

Author e-mails

Author affiliations

ORCID iDs

Dates

Abstract

1. Introduction