Measuring Turbulent Motion in Planet-forming Disks with ALMA: A Detection around DM Tau and Nondetections around MWC 480 and V4046 Sgr

Understanding the connection between disks around young stars and the planets formed within these environments requires knowledge about a wide range of physical processes, as well as how these processes vary between different stellar systems. One of the key disk properties is the nonthermal, non-Keplerian turbulent motion, which theoretically influences processes ranging from the growth of grains from micron sizes up to planets (e.g., Ormel & Cuzzi 2007; Birnstiel et al. 2010) and the orbital evolution of fully formed planets (e.g., Nelson & Papaloizou 2004; Paardekooper et al. 2011; Fung et al. 2014) to the chemical (Semenov & Wiebe 2011; Furuya & Aikawa 2014; Xu et al. 2017) and angular momentum (e.g., Lynden-Bell & Pringle 1974) evolution of the disk.

Efforts to broadly characterize the turbulence among a large sample of young stellar objects have relied on indirect techniques. Assuming that a planet-forming disk evolves through viscous processes, in which angular momentum is transported outward while mass is transported inward, leads to a predicted connection between accretion rate, disk mass, disk radius, and time (Hartmann et al. 1998) that can be used to constrain turbulence. In the context of an α-disk model of viscosity (Shakura & Sunyaev 1973), in which the viscosity (ν) is equal to the product of α, the local sound speed (c_s), and the gas pressure scale height (H), a number of studies have examined various predictions of viscous evolution using samples of ∼30–50 sources (Stepinski 1998; Lodato et al. 2017; Mulders et al. 2017; Rafikov 2017; Ansdell et al. 2018; Najita & Bergin 2018) and have found a wide range of α values (10⁻⁴ to 0.1). Mulders & Dominik (2012) rely on the connection between turbulence and the settling of dust grains (e.g., Youdin & Lithwick 2007) and use the average spectral energy distributions (SEDs) of Herbig stars, T Tauri stars, and brown dwarfs to find that α ∼ 10⁻⁴ best reproduces the average emission profile. High-resolution continuum imaging of dust structures also provides constraints on nonthermal motion through the sharpness and morphology of both dark gaps and bright rings (e.g., Pinte et al. 2016; Dullemond et al. 2018; Huang et al. 2018).

Gas motion provides a more direct measure of the turbulence, but observational constraints are limited to small samples. Hughes et al. (2011) studied the disks around HD 163296 and TW Hya, while Guilloteau et al. (2012) studied the disk around DM Tau. In the ALMA era, TW Hya and HD 163296 have been revisited (Flaherty et al. 2015; Teague et al. 2016; Flaherty et al. 2017, 2018; Teague et al. 2018a), finding upper limits on the turbulence of ≲5%–10% of the local sound speed for gas at distances of ∼30 au and larger from the central star. Near-infrared observations, more sensitive to gas within a few au of the central star, have found evidence of nonthermal broadening components comparable to the local sound speed among samples of similar sizes (Carr et al. 2004; Najita et al. 2009; Doppmann et al. 2011; Ilee et al. 2014).

Here we build on our previous efforts to measure turbulence directly from ALMA observations of molecular gas emission lines to include the disks around DM Tau, MWC 480, and V4046 Sgr. While still representing a modest sample, these sources were chosen to cover a range of accretion rates, X-ray luminosities, and far-ultraviolet (FUV) luminosities. This parameter space is particularly relevant for the magnetorotational instability (MRI; Balbus & Hawley 1998). The MRI relies on the coupling between gas and magnetic field to destabilize a rotating disk to generate the correlated turbulent fluctuations that can transport angular momentum and drive accretion through the disk. The outcome of the MRI, especially the level of turbulence, is sensitive to the level of disk ionization (e.g., Simon et al. 2018). Therefore, observationally establishing a relation between turbulence and other disk parameters offers important physical insight into the microphysical processes governing disk evolution.

In Section 2 we present the ALMA observations of these three systems, and in Section 3 we discuss our framework for modeling the CO emission. We present our finding of weak turbulence around MWC 480 and V4046 Sgr and nonzero turbulence around DM Tau in Section 4, while in Section 5 we examine the evidence that weak turbulence may be common among planet-forming disks, as well as the physical conditions that influence the difference in the measured level of turbulence between DM Tau and the other systems.

ALMA Band 6 observations of the disks around DM Tau, MWC 480, and V4046 Sgr were taken as part of project 2016.1.00724.S. These observations were split into short- and long-baseline components in order to capture the full spatial extent of the CO emission (out to ∼6'') at high spatial resolution (∼03). Observations were taken between 2016 December and 2017 July, with baselines ranging from 15 m to 2.6 km (Table 1).

The spectral windows were split between a continuum window centered at 232 GHz and three windows centered on CO (2–1), ¹³CO (2–1), and C¹⁸O (2–1), respectively. This paper focuses on CO (2–1), while future work will consider the isotopologues. For the spectral lines, the correlator was configured for 960 channels, each 61.035 kHz wide (∼80 m s⁻¹). The standard ALMA reduction scripts were used to produce calibrated visibilities. For each set of observations, two rounds of phase self-calibration and one round of amplitude self-calibration were performed, based on the high signal-to-noise ratio (S/N) continuum data, with the exception of the short-baseline DM Tau observations, where the amplitude self-calibration was excluded. The spectra were binned down by a factor of 2 (0.16 km s⁻¹).

Uncertainties in the flux calibration will result in small variations in the flux between the two sets of observations. To remove this effect, we compare the continuum visibilities with baselines <200 kλ among the two observational epochs. We find that the long-baseline observations are 1.06, 1.04, and 0.93 times brighter than the short baselines observations for DM Tau, MWC 480, and V4046 Sgr, respectively; the size of the amplitude differences is consistent with that expected from the uncertainty in the ALMA calibration (Butler 2012). These scale factors are applied to the model images before fitting to the long-baseline data.

The positional offset of the disk between the two epochs can vary based on, e.g., errors in the proper-motion correction or errors in the astrometric accuracy of ALMA. We estimate the positional offset between the short- and long-baseline observations by fitting an elliptical Gaussian to the continuum observations, and we apply a phase offset to each epoch to correct the disk center to the phase center. For MWC 480 the positional offsets are 73 mas, −86 mas in R.A. and decl. for the short-baseline data and 56 mas, −3 mas for the long-baseline data. For the DM Tau disk, the offsets in R.A. and decl. are 230 mas, 3.6 mas for the short-baseline data and 230 mas, −2.2 mas for the long-baseline data. No significant offset between the short- and long-baseline data was found for V4046 Sgr, although for both data sets the center of the disk was offset by 270 mas in R.A. and 17 mas in decl. from the phase center, which is included in the models. The observed locations of the disks are consistent with the stellar positions now available in the GAIA DR2 database (Gaia Collaboration et al. 2018). The amplitude differences between the short and long baselines are unaffected by these positional offsets.

Cleaned images were generated by combining the short- and long-baseline visibilities. Using natural weighting results in a beam that is ∼03–05 wide with an rms of 2–4 mJy beam^–1 per 0.16 km s⁻¹ wide channel.

The modeling code used to derive the turbulence has been described previously in Flaherty et al. (2015, 2017, 2018) and is based on earlier work by Rosenfeld et al. (2013a) and Dartois et al. (2003). We summarize the relevant equations below, including the key free parameters and the assumptions made for each source.

The surface density is assumed to follow a power law, with an exponential tail, as expected for a viscously evolving disk (Lynden-Bell & Pringle 1974; Hartmann et al. 1998):

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{\mathrm{gas}}(r)=\displaystyle \frac{{M}_{\mathrm{gas}}(2-\gamma )}{2\pi {R}_{c}^{2}}{\left(\displaystyle \frac{r}{{R}_{c}}\right)}^{-\gamma }\exp \left[-{\left(\displaystyle \frac{r}{{R}_{c}}\right)}^{2-\gamma }\right],\end{eqnarray} \tag{ 1 }$$

where M_gas, R_c, and γ are the gas mass (in M_⊙), critical radius (in au), and power-law index, respectively. The disk extends from R_in to 1000 au.

The temperature structure is defined as a power law with radius, with a vertical gradient connecting the cold midplane with the warm atmosphere:

$$\begin{eqnarray}&&{T}_{\mathrm{mid}}={T}_{\mathrm{mid}0}{\left(\displaystyle \frac{r}{150\mathrm{au}}\right)}^{q}\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&{T}_{\mathrm{atm}}={T}_{\mathrm{atm}0}{\left(\displaystyle \frac{r}{150\mathrm{au}}\right)}^{q}\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&\begin{array}{l}{T}_{\mathrm{gas}}(r,z)\\ \ =\ \left\{\begin{array}{ll}{T}_{\mathrm{atm}}+({T}_{\mathrm{mid}}-{T}_{\mathrm{atm}}){\left(\cos \displaystyle \frac{\pi z}{2{Z}_{q}}\right)}^{2} & {\rm{if}}\ z\lt {Z}_{q}\\ {T}_{\mathrm{atm}} & {\rm{if}}\ z\geqslant {Z}_{q}\end{array}\right.\end{array}\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&{Z}_{q}={Z}_{q0}{\left(r/150\mathrm{au}\right)}^{1.3}.\end{eqnarray} \tag{ 5 }$$

The parameter Z_q is the height above the midplane at which the gas temperature reaches its maximum value. Once the temperature structure and surface density profile have been specified, the hydrostatic equilibrium calculation is performed at each radius to derive the gas volume density at each height above the midplane. The velocity field is assumed to be Keplerian, with corrections for the height above the midplane and the pressure support of the gas, as in Rosenfeld et al. (2013a).

The line profile is a Gaussian whose width is set by the thermal and nonthermal motions. We associate the nonthermal term with turbulence (δv_turb), although we consider other sources of non-Keplerian, nonthermal motion in Section 4.1.1. We focus on cases in which the nonthermal term is proportional to the local isothermal sound speed:

$$\begin{eqnarray}&&{\rm{\Delta }}V=\sqrt{\left(2{k}_{B}T(r,z)/{m}_{{CO}}\right)(1+\delta {v}_{\mathrm{turb}}^{2})}.\end{eqnarray} \tag{ 6 }$$

This parameterization is physically motivated by numerical simulations that predict that the turbulence scales with the local sound speed (e.g., Shi & Chiang 2014; Simon et al. 2015; Flock et al. 2017).

One important component in fitting the CO observations is the midplane temperature. As demonstrated in Flaherty et al. (2018), an overestimate of the midplane temperature can lead to an underestimate of the turbulence. The data themselves do not provide strong constraints on the midplane temperature; C¹⁸O (2–1) has a strong degeneracy between T_mid and the CO abundance, while CO (2–1), unlike other systems (e.g., Rosenfeld et al. 2013a; Pinte et al. 2018a; Dullemond et al. 2019), does not have enough sight lines at low optical depth that reach to the midplane.

For DM Tau, the value of T_mid0 is chosen such that CO freeze-out at the midplane (T_mid = 19 K) occurs at 70 au, as predicted by the chemical models of Loomis et al. (2015). To preserve the location of CO freeze-out while the temperature structure, through the parameter q, is varied during the model fitting, the exact value of the T_mid0 parameter is defined based on the value of q (i.e., T_mid0 = 19(70/150)^−q). For MWC 480, we use the temperature structure derived by Piétu et al. (2007) based on observations of ¹³CO, which corresponds to T_mid0 = 18.8 K. For the disk around V4046 Sgr, we use a midplane normalization based on a D'Alessio model fitting of the SED, which results in T_mid0 = 12 K (C. Qi 2020, private communication). In addition to these fiducial values, we also consider lower T_mid0 values to examine the effect of a different assumed value of the midplane temperature on the turbulence estimates.

In our initial attempts at fitting the CO emission, we found that the models were unable to simultaneously reproduce the emission at radii ≲200 au and ≳200 au. This likely reflects complexities in the CO freeze-out/desorption process that are not included in our initial model prescription; we implicitly assume that thermal desorption and photodissociation set the CO molecular layer, while recent chemical models have shown that nonthermal desorption effects are important around IM Lup (Öberg et al. 2015) and TW Hya (Teague et al. 2017). Similarly, extended cold CO emission has been observed around AS 209 (Huang et al. 2016) and TW Hya (Schwarz et al. 2016). To account for nonthermal desorption of CO, we reintroduce CO to the disk in the region between ${N}_{{{\rm{H}}}_{2}}$ = 1.3 × 10²¹ cm⁻² and 4.8 × 10²¹ cm⁻² when T_gas < 20 K, as described in Appendix A.

In our fiducial models we allow T_atm0, q, R_c, i, δv_turb, and R_in to vary, along with the systemic velocity of the disk (v_sys). Table 2 lists the fixed parameters for each source. All distances are based on GAIA DR2 parallax measurements (Gaia Collaboration et al. 2018; Lindegren et al. 2018). The position angle of the disk is determined based on initial fitting of the system (PA = 1548, 1478, 7565 for DM Tau, MWC 480, and V4046 Sgr, respectively) and was kept fixed owing to its independence from the other model parameters, while γ is always fixed at 1. For the disk around V4046 Sgr, we use the sum of the masses of the two stars to set the Keplerian motion of the disk. The sign on the inclination is adjusted so as to minimize asymmetric residuals (Appendix B). Using GAIA DR2 distances, Simon et al. (2019) derive stellar masses for MWC 480 and DM Tau of 2.11 ± 0.01 M_⊙ and 0.55 ± 0.02 M_⊙, respectively. Utilizing these updated masses will most strongly affect the derived inclination (e.g., Czekala et al. 2016), while their affect on turbulence is likely to be small, especially for DM Tau, and will be investigated in more detail in future work.

Table 2.  Fixed Model Parameters

Star	Stellar Mass (M⊙)	Mgas (M⊙)	XCO	Distance (pc)	PA	Zq0 (au)	γ
DM Tau	0.54 (1)	0.04 (2,3)	2 × 10−5 (3)	144.5 ± 1.1 (4)	1548	70	1
MWC 480	1.85 (5)	0.046 (6)	10−4	161.1 ± 2.0 (4)	1478	50	1
V4046 Sgr	1.75 (7)	0.09 (8)	3 × 10−6 (7)	72.26 ± 0.34 (4)	757	35	1

References.  (1) Simon et al. 2000; (2) Andrews et al. 2011; (3) McClure et al. 2016; (4) Lindegren et al. 2018; Bailer-Jones et al. 2018; (5) Piétu et al. 2007; (6) this work; (7) Rosenfeld et al. 2012; (8) Rosenfeld et al. 2013b.

Download table as:  ASCIITypeset image

The posterior distributions for each parameter are estimated using the Markov Chain Monte Carlo (MCMC) routine EMCEE (Foreman-Mackey et al. 2013), based on the Affine-Invariant algorithm originally proposed in Goodman & Weare (2010). The uncertainties in the data are estimated based on the calculated dispersion among baselines of similar distances in line-free channels. Typical MCMC chains consist of 80 walkers and 800 steps, with convergence on the final solution occurring within 300 steps. The first 500 steps are removed as burn-in, after which the median of the posterior distributions varies by <1%. We use linear spacing in all parameters except R_c, where we fit log(R_c). We assume uniform priors that simply restrict parameters to physically realistic values (e.g., δv_turb > 0, T_atm0 > 0, R_in > 0).

4.1. Nonzero Turbulence around DM Tau

Unlike in the disks around HD 163296 and TW Hya, we find that a nonzero turbulence provides the best fit to the disk around DM Tau, with a 99.7% confidence interval on the turbulent motion of ${0.279}_{-0.004}^{+0.005}$c_s (Table 3), corresponding to velocities of 60–80 m s⁻¹. Model spectra, as well as residuals generated from the difference between model and data visibilities, are shown in Figure 1. Full posterior distributions, along with the progression of the walkers through the chains, are shown in Appendix C. The nonzero turbulence model spectrum is well matched to the data, and outside of a persistent unresolved (diameter <40 au) positive residual at the center, indicating that the model underestimates the data, the channel maps do not show strong residuals. As discussed below, similar central residuals are seen around MWC 480 and V4046 Sgr and may be due to deviations from our simple prescriptions for density and temperature, or a sign of changing CO abundance with radius.

Zoom In Zoom Out Reset image size

Table 3.  CO (2–1) Model Results

Model	q	$\mathrm{log}({R}_{c}(\mathrm{au}))$	δvturb (cs)	Tatm0 (K)	Tmid0 (K)	Inclination (deg)	Rin (au)	vsys (km s−1)
DM Tau
Fiducal	$-{0.371}_{-0.002}^{+0.003}$	${2.444}_{-0.003}^{+0.004}$	${0.279}_{-0.004}^{+0.005}$	${24.68}_{-0.09}^{+0.11}$	14.3a	${36.0}_{-0.09}^{+0.12}$	<9	6.017 ± 0.001
sys = 1.2	−0.344 ± 0.003	${2.452}_{-0.004}^{+0.003}$	${0.257}_{-0.005}^{+0.004}$	${28.08}_{-0.12}^{+0.15}$	14.6a	${36.21}_{-0.13}^{+0.12}$	<9	${6.018}_{-0.002}^{+0.001}$
sys = 0.8	$-{0.403}_{-0.002}^{+0.003}$	2.44 ± 0.01	0.300 ± 0.005	${21.26}_{-0.08}^{+0.09}$	14.0a	${35.81}_{-0.13}^{+0.12}$	<9	6.016 ± 0.001
highres	$-{0.352}_{-0.002}^{+0.001}$	2.440 ± 0.002	${0.252}_{-0.004}^{+0.003}$	${24.92}_{-0.07}^{+0.08}$	14.5a	${35.99}_{-0.08}^{+0.09}$	<9	6.045 ± 0.001
$\delta {v}_{\mathrm{turb}}$ = 0	$-{0.352}_{-0.002}^{+0.004}$	${2.456}_{-0.003}^{+0.004}$	0b	27.54 ± 0.09	14.5a	${37.73}_{-0.11}^{+0.16}$	<9	${6.013}_{-0.001}^{+0.002}$
low ${T}_{\mathrm{mid}0}$	$-{0.310}_{-0.003}^{+0.002}$	${2.469}_{-0.005}^{+0.004}$	0.328 ± 0.004	${25.57}_{-0.07}^{+0.09}$	11.0c	35.86 ± 0.11	<9	6.019 ± 0.001
MWC 480
Fiducial	−0.467 ± 0.002	${2.161}_{-0.001}^{+0.004}$	<0.12	${60.3}_{-0.2}^{+0.3}$	18.8b	$-{40.34}_{-0.06}^{+0.05}$	<9	5.102 ± 0.001
highres	−0.457 ± 0.002	${2.172}_{-0.001}^{+0.002}$	<0.08	62.1 ± 0.2	18.8b	$-{40.37}_{-0.04}^{+0.05}$	<9	5.137 ± 0.001
$\delta {v}_{\mathrm{turb}}$ = 0	−0.463 ± 0.002	${2.168}_{-0.003}^{+0.002}$	0b	62.3 ± 0.1	18.8b	$-{40.49}_{-0.05}^{+0.06}$	<9	5.102 ± 0.001
low ${T}_{\mathrm{mid}0}$	$-{0.413}_{-0.003}^{+0.001}$	${2.208}_{-0.003}^{+0.004}$	<0.18	${70.02}_{-0.05}^{+0.36}$	14b	−40.06 ± 0.05	<9	5.141 ± 0.001
V4046 Sgr
Fiducial	−0.565 ± 0.003	${1.939}_{-0.003}^{+0.004}$	<0.15	${27.84}_{-0.09}^{+0.13}$	12b	${34.69}_{-0.04}^{+0.03}$	<7	2.857 ± 0.001
highres	$-{0.564}_{-0.001}^{+0.003}$	${1.939}_{-0.003}^{+0.001}$	<0.12	${28.26}_{-0.05}^{+0.10}$	12b	${34.68}_{-0.03}^{+0.02}$	<7	${2.8957}_{-0.0006}^{+0.0005}$
δvturb = 0	−0.579 ± 0.002	${1.956}_{-0.003}^{+0.004}$	0b	${28.80}_{-0.07}^{+0.04}$	12b	${34.76}_{-0.03}^{+0.04}$	<7	${2.857}_{-0.001}^{+0.001}$
low Tmid0	$-{0.565}_{-0.004}^{+0.002}$	${2.078}_{-0.002}^{+0.012}$	<0.18	${29.27}_{-0.06}^{+0.17}$	8b	${34.59}_{-0.03}^{+0.05}$	<6.5	2.857 ± 0.001

Notes. sys = 1.2: a trial simulating if the data were systematically 20% brighter than measured; sys = 0.8: a trial simulating if the data were systematically 20% fainter than measured; highres: a trial fitting to the data at full spectral resolution, rather than the binned spectra used in the fiducial model; v_turb = 0: a trial in which turbulence is fixed at zero; low T_mid0: a trial in which the midplane gas temperature is lower than the fiducial model.

^aT_mid0 = 19(70/150)^−q. This forces the midplane temperature to 19 K at 70 au, matching the location of the CO condensation front in the chemical modeling of Loomis et al. (2015). ^bHeld fixed during the fitting process. ^cT_mid0 = 14(70/150)^−q.

Download table as:  ASCIITypeset image

The MCMC modeling indicates that a nonzero turbulence exists at a highly significant level. To verify this conclusion, we perform an additional MCMC run with the turbulence fixed at zero. This results in a model that systematically underestimates the flux of the disk (Figure 1). Comparing the median models from both trials using the Akaike Information Criterion (Akaike 1974), the difference is confirmed, with the zero turbulence model being worse than the nonzero turbulence model at a level that is much greater than 10σ. We also can confirm the nonzero turbulence by fitting to the unbinned spectra while allowing turbulence to vary. Again we find nonzero turbulence (${0.252}_{-0.004}^{+0.003}$c_s), confirming our fit to the binned spectrum.

We find an uncertainty on the turbulence at the ∼2% level, but this does not include the systematic uncertainty on the amplitude calibration, which may affect the derived temperature structure, and in turn the turbulent broadening. The 4%–7% variation in amplitude measured between the short- and long-baseline observations of the three systems under study here is consistent with some uncertainty in the flux calibration of these data. To test the robustness of our result in the face of the amplitude calibration uncertainty, we run additional trials in which the model is systematically scaled by ±20%, simulating a true disk flux that is 20% fainter/brighter than is actually observed. Again we find significant turbulence, with δv_turb = 0.300c_s ± 0.005c_s for a faint disk and δv_turb = ${0.257}_{-0.005}^{+0.004}$c_s for a bright disk. These results differ in the gas temperature with T_atm0 = ${21.26}_{-0.08}^{+0.09}$ K for the faint disk and T_atm0 = ${28.08}_{-0.12}^{+0.15}$ K for the bright disk, which in turn affects the exact value of the sound speed. This anticorrelation between T_atm0 and δv_turb is also seen in the fiducial model posterior distribution (Figure C2). The range in δv_turb values between the systematic uncertainty models is larger than in the fiducial model since the amplitude calibration uncertainty allows for a larger range in T_atm0. The anticorrelation between δv_turb and T_atm0 does indicate that uncertainty in the gas temperature structure influences our turbulence constraints.

As discussed earlier, the assumed midplane temperature also influences the final turbulence result. We run an additional MCMC ensemble in which T_mid0 has been reduced to 11 K, as suggested by observations of N₂H⁺ (Qi et al. 2019), and lower than employed in the chemical models of Loomis et al. (2015). Again we find nonzero turbulence, at a level of δv_turb = 0.328c_s ± 0.004c_s. The anticorrelation between midplane temperature and turbulence arises because both influence the strength of the emission at the edges of the emitting regions within a given channel map (Flaherty et al. 2018). These regions have low optical depth and probe the gas temperature closer to the midplane; increasing the midplane temperature will increase the intensity of the emission from these edges. At the same time, increasing turbulence will increase the amount of emission that "bleeds over" from nearby channels, leading to an increase in the optical depth and emission that arises from higher and warmer layers. As a result, a given intensity along the edges of the emission within a specific channel can be achieved with either a cold midplane and strong turbulence or weaker turbulence and a warmer midplane. When combined with the results from the amplitude calibration trials, the low-T_mid0 trial indicates that systematic uncertainties dominate over statistical uncertainties in our measurement of turbulence. Based on the sample of trials run here, we find δv_turb = 0.25c_s–0.33c_s, significantly larger than zero, but with a larger uncertainty than indicated by any individual trial.

Previous analysis of CS (3–2) and CS (5–4) (Guilloteau et al. 2012) found significant turbulence around DM Tau, at a level of 120 m s⁻¹, with CS being used instead of CO because of its larger molecular weight and hence smaller thermal broadening. When including the uncertainty on the turbulence velocity and the temperature structure, this translates to a 3σ range of (0.12–0.69)c_s, consistent with our measured range of (0.25–0.33)c_s when accounting for systematic uncertainties (Figure 2). The higher-S/N and higher spatial resolution ALMA data provide a tighter constraint on the turbulence around DM Tau.

Zoom In Zoom Out Reset image size

Interpreting these results in the context of theoretical models requires knowledge of the emitting region of the molecular emission. Figure 2 shows the emitting region of CO derived from our models, as well as that of CS from Guilloteau et al. (2012); the emission from both molecules arises from between two and three pressure scale heights above the midplane. The emitting regions are at low column densities, where FUV photons are expected to ionize the disk, making it more susceptible to MRI (Perez-Becker & Chiang 2011). At δv_turb = 0.25c_s–0.33c_s, the turbulence is at the low end for model predictions from MRI (Simon et al. 2015, 2018), consistent with gravito-turbulence (Forgan et al. 2012; Shi & Chiang 2014; Shi et al. 2016), and at the high end for the vertical shear instability (VSI; Richard et al. 2016; Flock et al. 2017). Distinguishing between these different models for the source of the turbulence may require a measure of the vertical structure of the turbulence (Simon et al. 2015), more detailed knowledge of the physical conditions within the outer disk (Simon et al. 2018; Lyra & Umurhan 2019), or complementary observations of other disk properties that are influenced by the turbulence (e.g., dust settling; Flock et al. 2017).

4.1.1. Anisotropic versus Isotropic Motion

Within our models we assume that the nonthermal non-Keplerian component of the motion is isotropic turbulence, but other forms of non-Keplerian motion may exist in protoplanetary disks. Deviations are predicted to be induced by planets (Bae & Zhu 2018; Pérez et al. 2018), and recent observational studies have found evidence for deviations from Keplerian motion due to the gravitational influence of unseen planets (Pinte et al. 2018b, 2019; Teague et al. 2018b, 2018c). Zonal flows, radial density variations associated with magnetic field concentrations, and the VSI, arising in part from the change in Keplerian velocity with vertical distance from the midplane, produce corrugated variations in the gas motions on velocity scales that are comparable to the nonthermal non-Keplerian motion observed around DM Tau (e.g., Simon & Armitage 2014; Flock et al. 2017). If deviations from Keplerian motion are spread through a large region of the outer disk and are of a sufficiently high velocity at a sufficiently small spatial scale, they may mimic the appearance of isotropic turbulence.

To test the sensitivity of our data to structured motion, we employ a toy model in which we add a velocity component in the orbital direction that varies as a function of radius:

$$\begin{eqnarray}&&{{dv}}_{\phi }={dv}\sin \left(\displaystyle \frac{2\pi R}{{dR}}\right).\end{eqnarray} \tag{ 7 }$$

Here dv represents the maximum deviation from Keplerian motion, while dR is the spatial scale of the perturbation. This perturbation is added to a model with parameters (T_atm0, R_c, etc.) based on the DM Tau fiducial model. Visibilities are generated from the model images and are used to create a cleaned image, using the same procedure as for the data.

In Figure 3 we show the central velocity channel of three models: (1) large spatial scale and high velocity (dR = 150 au, dv = 1.5 km s⁻¹), (2) small spatial scale and high velocity (dR = 70 au, dv = 1.5 km s⁻¹), and (3) large spatial scale and small velocity (dR = 150 au, dv = 0.3 km s⁻¹ ∼ c_s). Anisotropic motion is easily distinguishable from isotropic turbulence in the large size scale and large velocity model owing to the corrugated emission structure but becomes more difficult to detect as the size and velocity scales become smaller. When the distance between the minimum and maximum velocity deviation (dR/2) becomes smaller than the beam size (∼40 au for DM Tau), the features are no longer spatially resolved and blend together into a general broadening of the emission. For small velocity scales the azimuthal perturbations to the image are small enough that they are not resolvable as distinct features. Leveraging the full three-dimensional data cube, rather than a single channel as shown in Figure 3, provides greater sensitivity to smaller fluctuations (Teague et al. 2018b, 2018c), although it is still limited by the spatial resolution of the data.

Zoom In Zoom Out Reset image size

Theoretical models of structured non-Keplerian motion suggest that the radial and velocity scales of typical features are small enough to make them difficult to distinguish from purely isotropic motion in our data. The size scale of these features generally scales with the pressure scale height (H), and for DM Tau our model fits indicate a pressure scale height that reaches 70 au at 500 au from the central star. With a beam size of ∼40 au, features as large as (5–10)H would be easily resolved. Global nonideal MHD simulations find that when the net poloidal field is relatively strong (thus driving high accretion rates), the disk tends to form zonal flows of scales¹³ of ∼2H (Béthune et al. 2017; Suriano et al. 2018), while in shearing box simulations, zonal flows are found at a broader range of disk magnetizations, on scales of 4H or beyond, and are dependent on the box size (Bai & Stone 2014; Simon & Armitage 2014) and the magnetization of the disk (Riols & Lesur 2019). VSIs produce corrugated variations in the gas motion with velocities of 10%–20% of the local sound speed and size scales of ∼2H (Flock et al. 2017). At these scales, zonal flows and VSI are difficult to distinguish from isotropic motion in our data, although other factors, such as the efficient lifting of dust grains by VSI (Flock et al. 2017), may break this degeneracy. Planet-induced motions are unlikely to be misinterpreted as turbulence because of the highly localized nature of the perturbation (Pérez et al. 2018). In modeling the CO emission from around DM Tau we assume that a single turbulence value, relative to the local sound speed, applies uniformly throughout the disk. The residuals do not show a strong deviation from this assumption, outside of the unresolved emission at the center of the disk. This suggests that instabilities that generate features over a wide radial range within the disk, such as is possible under certain physical conditions with zonal flows and VSIs, as well as MRI assuming that the surface layers are sufficiently ionized, are consistent with our results around DM Tau.

4.2. Weak Turbulence around MWC 480 and V4046 Sgr

For both the disk around MWC 480 and the disk around V4046 Sgr we do not detect significant nonthermal motion of the gas. Figures 4 and 5 show the data as compared to the model defined by the median of the posterior distributions (posterior distributions are shown in Appendix C). For both MWC 480 and V4046 Sgr, the MCMC routine returns a posterior that nominally indicates a detection of turbulence. We run additional trials with turbulence fixed at zero and find that the zero turbulence fit is nearly indistinguishable, in both the spectra and the channel maps, from the nonzero turbulence result, although, as discussed below, there are still residuals between the best-fit models and the data. For this reason we conservatively interpret the results for MWC 480 and V4046 Sgr as upper limits rather than detections.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For MWC 480, we are able to constrain turbulence to <0.12c_s. Fitting to the unbinned data pushes the limit down to <0.08c_s, corresponding to velocities below 20–40 m s⁻¹ beyond 100 au. For V4046 Sgr a fit to the binned spectrum constrains turbulence to <0.15c_s, while a fit to the unbinned spectrum pushes the upper limit down to <0.12c_s, corresponding to velocities below 30–55 m s⁻¹ beyond 100 au. As with DM Tau, we consider models with colder midplane temperatures (T_mid0 = 14 K for MWC 480, and 8 K for V4046 Sgr), and in both cases the upper limits rise to <0.18c_s but are still consistent with a nondetection of turbulence. This suggests that systematic effects (e.g., the assumed midplane temperature) are significant compared to statistical uncertainties when constraining the turbulence level in these systems.

In the case of both MWC 480 and V4046 Sgr, the imaged residuals show structure (Figures 4 and 5), although at a level that is less than 10% of the peak flux. Both systems show an unresolved positive residual at the center of the disk, with a diameter <50 au and <20 au for MWC 480 and V4046 Sgr, respectively, which may be a result of a deviation from simple prescriptions for the temperature and surface density radial profiles, a change in CO abundance with radius (e.g., Schwarz et al. 2018), or a change in the velocity profile due to a warp (e.g., Walsh et al. 2017) or rapid inward radial flow (e.g., Rosenfeld et al. 2014). Around V4046 Sgr negative residuals extend out to ∼200 au, corresponding to the radius beyond which our simple photodesorption treatment sets the CO abundance. This suggests that a more detailed treatment of photodesorption would bring the model in closer agreement with the data. None of these residuals are indicative of an underestimate of the turbulence level. In the channel maps turbulence acts to broaden emission rather than to change the radial profile (Simon et al. 2015). Also, the similarity between the model derived when turbulence was fixed at zero and the model derived when turbulence was allowed to vary supports our interpretation of the model results as an upper limit. We can rule out turbulence levels comparable to that around DM Tau, but we cannot completely rule out modest turbulent motions within these systems.

5.1. Is Weak Turbulence Common among Planet-forming Disks?

5.2. Physical Conditions for Turbulence

Physical conditions within the disk can set the turbulence level and may explain the lack of turbulence among many of the studied sources. In the context of MRI, Simon et al. (2018) explore the turbulence levels among models with different ionization conditions and magnetic field strengths. They find that the presence of FUV emission drives strong turbulence (0.2c_s–0.7c_s), while the removal of this ionization component, especially when combined with weak magnetic fields, can lead to weak turbulence (0.02c_s–0.09c_s). Among the sources with constraints on turbulence from molecular line emission, DM Tau exhibits the weakest FUV flux (Table 4) despite having the only turbulent disk. X-ray emission can also ionize gas, but again DM Tau does not possess the strongest X-ray emission among our sample. This suggests that FUV and X-ray emission strength alone does not set the turbulence level.

Table 4.  System Properties

Star	δvturb	M*	Mdisk	Age	${\dot{M}}_{\mathrm{acc}}$	LX-ray	LFUV
		(M⊙)	(M⊙)	(Myr)	(M⊙ yr−1)	(L⊙)	(L⊙)
HD 163296	<0.05cs (1)	2.3	0.09 (2)	5	5 × 10−7 (3)	10−4 (4, 5)	3.21–5.58 (6)
TW Hya	<0.08cs (7)	0.6	0.05 (8)	10–12	2 × 10−9 (9, 10, 11)	6.8 × 10−4 (12, 13, 14, 15)	7 × 10−3 (16)
MWC 480	<0.08cs (17)	1.85 (18)	0.046 (17)	7–8	5.3 × 10−7 (19, 3, 20)	2.6 × 10−6 (21)	...
V4046 Sgr	<0.12cs (17)	0.9,0.9	0.09 (22)	12–23 (23, 24, 25)	1.3 × 10−8 (26)	3.1 × 10−4 (27, 28)	10−2 (16)
DM Tau	0.25–0.33cs (17)	0.54 (29)	0.04 (30,31)	1–5	2.9 × 10−9 (16, 32)	7.8 × 10−5 (33)	3 × 10−3 (34, 16)
Typical disk ranges	...	...	3 × 10−4–0.03 (35)	1–10	5 × 10−11–10−7 (36)	2 × 10−6–0.02 (37)	2 × 10−3–0.2 (16)

References. (1) Flaherty et al. 2017; (2) Isella et al. 2007; (3) Mendigutía et al. 2013; (4) Günther & Schmitt 2009; (5) Swartz et al. 2005; (6) Meeus et al. 2012; (7) Flaherty et al. 2018; (8) Bergin et al. 2013; (9) Alencar & Batalha 2002; (10) Herczeg et al. 2004; (11) Ingleby et al. 2013; (12) Brickhouse et al. 2010; (13) Dupree et al. 2012; (14) Kastner et al. 1999; (15) Kastner et al. 2002; (16) France et al. 2014; (17) this work; (18) Piétu et al. 2007; (19) Costigan et al. 2014; (20) Donehew & Brittain 2011; (21) Henning et al. 2010; (22) Rosenfeld et al. 2013b; (23) Mamajek & Bell 2014; (24) Torres et al. 2006; (25) Binks & Jeffries 2014; (26) Curran et al. 2011; (27) Espaillat et al. 2013; (28) Sacco et al. 2012; (29) Simon et al. 2000; (30) Andrews et al. 2011; (31) McClure et al. 2016; (32) Ingleby et al. 2011; (33) L.I. Cleeves 2020, private communication; (34) Yang et al. 2012; (35) Andrews et al. 2013; (36) Manara et al. 2016; (37) Feigelson et al. 2005.

Download table as:  ASCIITypeset image

Observations constrain the FUV and X-ray emission near the stellar surface, while the relevant parameter for MRI models is the FUV and X-ray field as seen by material in the outer disk. High-energy radiation can be severely attenuated before reaching the outer disk if intervening material, such as a wind launched by the inner disk, acts as an obscuring screen. TW Hya (Pascucci et al. 2011), MWC 480 (Fernandes et al. 2018), V4046 Sgr (Sacco et al. 2012), and HD 163296 (Klaassen et al. 2013) show evidence of a disk wind based on optical/infrared forbidden lines (TW Hya, V4046 Sgr), infrared emission (MWC 480), and extended CO emission (HD 163296). In contrast, the optical forbidden line emission from DM Tau can be explained entirely with a gas disk in Keplerian motion (Simon et al. 2016). Whether or not a wind is the defining factor between a turbulent and nonturbulent disk depends on the density of the wind, with vertical column densities of 0.01 g cm⁻² sufficient to absorb a majority of the FUV radiation (Perez-Becker & Chiang 2011).

Regardless of the ionizing source, the ionization level of the outer disk can be constrained based on observations of ionized molecular species. The presence of, e.g., DCO⁺, HCO⁺, and N₂H⁺ in the outer disk of HD 163296 (Flaherty et al. 2017), TW Hya (Qi et al. 2013), DM Tau (Dutrey et al. 1997; Öberg et al. 2011b; Teague et al. 2015), MWC 480 (Öberg et al. 2010; Huang et al. 2017), and V4046 Sgr (Öberg et al. 2011a; Huang et al. 2017) suggests that some ionizing radiation reaches the outer disk. Öberg et al. (2011b) measure an ionization fraction of ∼4 × 10⁻¹⁰ in the CO layer around DM Tau, comparable to the ionization fraction in the models from Simon et al. (2018) that assumed no FUV photons but still included X-rays and cosmic-rays. A subset of these models displayed vigorous turbulence in the upper disk layers.

Cleeves et al. (2017) found that the H¹³CO⁺ abundance around IM Lup varied in response to an X-ray flare, indicating that the ionization level may change on short timescales. In the presence of variable ionization, MRI will respond with varying turbulence levels on its growth/decay timescale. Since the MRI relies on the differential rotation of gas, its growth/decay timescale is, at its fastest, the dynamical timescale of the gas. Molecular line observations with ALMA probe regions of the disk beyond ∼30 au, corresponding to dynamical timescales of hundreds to thousands of years. In their chemical modeling, Cleeves et al. (2017) find that the HCO⁺ abundance would be enhanced for 10–20 days at 100 au following an X-ray flare, much shorter than the time needed for MRI to respond to changes in the ionization, resulting in little variation in turbulence with time. Longer-timescale fluctuations in, e.g., the wind mass-loss rate may still lead to oscillations in the turbulence level over the course of the disk lifetime.

For reasonable values of the magnetic field strength, as informed by theoretical models (e.g., Bai & Stone 2011; Bai 2016; Simon et al. 2018), stronger fields generally result in stronger turbulence, and possibly higher accretion rates. The predicted correlation between turbulence and accretion rate seems to be contradicted in our sample, with DM Tau showing the strongest turbulence at yet the weakest accretion rate. This may be due to the fact that the measured accretion rate is the instantaneous accretion rate onto the star and does not reflect the accretion rate in the outer disk where the turbulence is measured, or that the turbulence signature does not reflect MRI-generated turbulence but instead reflects other sources of nonthermal motion (e.g., unresolved spatial variations in azimuthal velocity).

Age may also play a factor in the strength of the turbulence by driving the evolution of the physical conditions that set the turbulence level. DM Tau is one of the youngest targets within our sample; as a member of the Taurus star-forming region, its age likely lies between 1 Myr (Luhman et al. 2010) and 5 Myr (Simon et al. 2000; Guilloteau et al. 2014). Larger turbulence at younger ages is consistent with the findings of Najita & Bergin (2018) of large α values among class I objects. But ages of individual objects are difficult to constrain; DM Tau and MWC 480 are both part of the Taurus region, and would be expected to have similar ages, yet MWC 480 is typically assumed to be older than DM Tau (Simon et al. 2000). Accounting for internal magnetic fields can push stellar ages to systematically higher values (Simon et al. 2019). Similarly, the lack of turbulence around the class I object HL Tau (Pinte et al. 2016) suggests that not all objects of the same age have the same turbulence properties.

In the absence of MRI, other processes may drive turbulence under the right conditions. Forgan et al. (2012) find that gravitational instabilities can drive large-scale accretion, with the driving source transitioning from large-scale spiral arms to gravito-turbulence as the disk mass decreases. The measured disk masses (Table 4) are 4%–8% of the stellar mass, resulting in minimum Toomre Q values of 5–15, inconsistent with the presence of gravito-turbulence. Powell et al. (2019) derive larger disk masses for HD 163296 and TW Hya based on the structure of the dust emission, but they find that even these large masses do not push these disks into the gravitationally unstable regime. The weak levels of turbulence may also open up the possibility of hydrodynamic instabilities, which can produce α ∼ 10⁻⁴ to 10⁻³, playing a larger role (e.g., VSI, convective overstability, and zombie vortex instability, as reviewed by Lyra & Umurhan 2019). The VSI develops owing to the change in orbital velocity with vertical height within a disk that has sufficiently rapid cooling (Simon & Armitage 2014; Flock et al. 2017), conditions that are likely in the outer disk (Lin & Youdin 2015). VSI can operate in a magnetized disk (Cui & Bai 2019) and produces highly anisotropic motion, dominated by vertical motion reaching 10%–20% of the local sound speed (Flock et al. 2017). Such motions can be ruled out for some of the sample (HD 163296, TW Hya, MWC 480). Relying on the ability of turbulence to spread dust in the radial direction, Dullemond et al. (2018) use the resolved widths of narrow dust rings around AS 209, Elias 24, HD 163296, GW Lup, and HD 143006 to rule out α ≪ 5 × 10⁻⁴ for grain sizes ≫0.1 cm, suggesting that some modest level of turbulence may exist in these systems.

Using new ALMA observations of CO (2–1) emission, we detect turbulence around DM Tau at a level of δv_turb = ${0.279}_{-0.004}^{+0.005}$c_s, although this range extends from 0.25c_s to 0.33c_s when accounting for systematic effects (e.g., assumptions about the midplane temperature). Around MWC 480 and V4046 Sgr we find no strong evidence for nonzero turbulence, restricting the turbulence levels to <0.08c_s and <0.12c_s, respectively. The finding of larger turbulence around DM Tau than around V4046 Sgr is consistent with the enhanced dust settling in the disk around V4046 Sgr relative to the disk around DM Tau (Qi et al. 2019).

The upper limits on turbulence in the disks around MWC 480 and V4046 Sgr put them at a turbulence level similar to that around HD 163296 (Boneberg et al. 2016; Flaherty et al. 2017; Dullemond et al. 2018; Liu et al. 2018) and TW Hya (van Boekel et al. 2017; Flaherty et al. 2018; Teague et al. 2018a). While this is still a modest sample, the results are consistent with studies of high-resolution images of continuum emission, sensitive to disk properties closer to the midplane and smaller radii than the CO emission, which found modest turbulence among additional planet-forming systems (α ∼ 10⁻⁴; Pinte et al. 2016; Huang et al. 2018; Pérez et al. 2019), suggesting that weak turbulence may be a common feature among planet-forming disks. At the same time, the robust detection of turbulence around DM Tau provides an important anchor point in comparison with models of instabilities within planet-forming disks.

Among this small sample we can begin to look for clues as to the explanation for these results. In the context of MRI, the strength of the ionizing radiation reaching the outer disk and the magnetic field are two defining factors (e.g., Simon et al. 2018). While the FUV and X-ray flux emitted by DM Tau is not large relative to MWC 480, V4046 Sgr, HD 163296, or TW Hya, DM Tau is the only source among this sample that does not show evidence for an inner disk wind that could potentially block the ionizing radiation from reaching the outer disk. The presence of ionized molecular species in all of these systems indicates that some modest level of ionization exists in the outer disk, although the exact ionization level has only been quantified in the case of the disk around DM Tau. Instead, the magnetic field strength may be the defining factor in setting the strength of the turbulence. Magnetic field strengths remain largely unconstrained, with observational upper limits (e.g., Vlemmings et al. 2019) still lying above the levels examined by theoretical models. Age also may play a role, given that DM Tau is the youngest source among our sample, although further measurements are needed to determine whether this trend holds up.

We thank the referee for their detailed reading and helpful comments. The work of K.F. was performed in part at the Aspen Center for Physics, which is supported by the National Science Foundation grant PHY-1607611. A.M.H. is supported by a Cottrell Scholar Award from the Research Corporation for Science Advancement. This paper makes use of the following ALMA data: ADS/JAO.ALMA#2016.1.00724.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. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program under grant agreement No. 716155 (SACCRED).

Facility: ALMA. -

Software: astropy (Astropy Collaboration et al. 2013), EMCEE (Foreman-Mackey et al. 2013), MIRIAD.

To incorporate CO in regions with photodesorption into our parametric structure, without resorting to detailed chemical models, we compare the timescales for freeze-out, photodesorption, and photodissociation. Regions in which photodesorption operates more quickly than freeze-out and photodissociation will provide additional CO emission in the outer disk. Regions in which freeze-out or photodissociation dominates will be heavily depleted in CO and will not contribute CO emission.

The freeze-out timescale is set by the frequency at which CO molecules collide with dust particles (Flower et al. 2005):

$$\begin{eqnarray}&&{\tau }_{\mathrm{fre}}={\left({n}_{g}\pi {a}_{g}{v}_{\mathrm{th}}S\right)}^{-1},\end{eqnarray} \tag{ A1 }$$

where n_g is the number density of grains, a_g is the radius of a grain particle, v_th is the thermal velocity of the gas, and S is the sticking efficiency, assumed to be 1.

CO photodesorption is influenced by the UV radiation field, its ability to penetrate the disk, and the intrinsic efficiency with which CO desorbs from a dust grain in response to an incoming photon:

$$\begin{eqnarray}&&R={I}_{\mathrm{ISRF}}\exp (-\gamma {A}_{V}){Y}_{\mathrm{pd}}{a}_{g},\end{eqnarray} \tag{ A2 }$$

where R is the number of CO molecules that are photodesorbed per second per dust grain, I_ISRF is the interstellar radiation field (=1 Habing = 10⁸ photons s⁻¹ cm⁻²; Habing 1968), γ is a measure of the UV extinction relative to the visual extinction (∼2 for small dust grains), A_V is the visual extinction, and Y_pd is the experimentally measured number of photodesorptions per UV photon. A detailed calculation of the interstellar radiation field would involve accounting for the nearby stellar population (e.g., Cleeves et al. 2016), while here we assume a radiation field of 1 Habing for simplicity. The factor Y_pd can vary by ∼3–4 depending on the exact shape of the incident spectrum (Fayolle et al. 2011; Chen et al. 2014) and the temperature at which the ice is deposited on the dust grains (Öberg et al. 2009), although we ignore these effects in our calculation and assume Y_pd = 2.7 × 10⁻³. Given the photodesorption rate (R), the timescale to double the CO density is

$$\begin{eqnarray}&&{\tau }_{\mathrm{pde}}=\displaystyle \frac{{n}_{\mathrm{CO}}}{{{Rn}}_{g}},\end{eqnarray} \tag{ A3 }$$

where n_CO is the number density of CO molecules.

Once CO is released from the grains, it must survive as a molecule, without being dissociated by the same UV radiation field that released it from the dust grain in the first place. The rate at which CO is photodissociated is (Visser et al. 2009)

$$\begin{eqnarray}&&k=\chi {k}_{0}{\rm{\Theta }}\exp (-\gamma {A}_{V}),\end{eqnarray} \tag{ A4 }$$

where χ is the radiation field in units of the Draine (1978) field (which corresponds to roughly 1.7 Habings), k₀ is the unattenuated photodissociation rate (=2.592 × 10⁻¹⁰, in units of photodissociations per second per CO molecule), and Θ is the attenuation factor accounting for shielding of the radiation field by H₂ and self-shielding by CO. The photodissociation timescale is simply τ_pdi = 1/k.

Using the disk around V4046 Sgr as a reference, with the parameters taken from Rosenfeld et al. (2013b) and assuming 0.1 μm dust grains and a dust-to-gas-ratio of 10⁻³, we find a region in the upper layers of the outer disk in which photodesorption operates more quickly than freeze-out and photodissociation (Figure B1). Near the disk surface, photodesorption operates more quickly than freeze-out owing to the strong UV field and the low densities. In the uppermost layers photodissociation operates more rapidly than photodesorption; while both processes rely on the same UV field, the decreasing number of dust grains from which CO can desorb, as well as the diminishing effects of self-shielding, leads photodissociation to operate more quickly than photodesorption in the highest layers. These trends limit effective photodesorption to a narrow layer within the disk, similar to what has been seen in more detailed chemical modeling (Öberg et al. 2015; Teague et al. 2017). The boundaries for the photodesorption region appear to follow contours in the vertically integrated column density. Given this behavior, we parameterize CO nonthermal desorption as a return to normal CO abundance in the region between ${N}_{{{\rm{H}}}_{2}}$ = 1.3 × 10²¹ cm⁻² and 4.8 × 10²¹ cm⁻². The CO emission is most strongly dependent on the boundary at lower column densities, corresponding to larger heights above the midplane, due to the highly optically thick nature of CO. Our choice of ${N}_{{{\rm{H}}}_{2}}$ = 1.3 × 10²¹ cm⁻² for photodesorption matches the boundary used in the inner disk for photodissociation of CO, ensuring that there is no discontinuity in the CO location across the freeze-out boundary. Including this parameterization for CO photodesorption eliminates most of the residuals in the outer disk (Figure B1).

For a perfectly thin disk, there is a degeneracy between positive and negative inclinations, but protoplanetary disks have a nonzero thickness, which introduces small differences in the emission based on the sign of the inclination. In the case of the disk around HD 163296, Rosenfeld et al. (2013a) demonstrated that spatially resolving the upper and lower molecular layers of the disk indicates that the inclination is positive. While we do not resolve the layers of the disks around DM Tau, MWC 480, or V4046 Sgr, the nonzero thickness still influences the emission. Lines of sight through the half of the disk closer to the observer terminate at smaller radii, and hence warmer gas temperatures, than lines of sight through the far side of the disk. This leads to a small asymmetry in the emission from the disk that can be detected when the data are of sufficiently high S/N.

Zoom In Zoom Out Reset image size

Figure C1 shows the residuals for the three systems between models with a positive and negative inclination (following the disk orientation convention outlined in Czekala et al. 2019). The disks around DM Tau and V4046 Sgr exhibit residuals that are asymmetric along the minor axis of the disk when the inclination is negative, while these residuals disappear when the inclination is positive. Conversely, the disk around MWC 480 has asymmetric residuals for a positive inclination but more symmetric residuals for negative disk inclinations. This is evidence that the disks around DM Tau and V4046 Sgr have positive inclinations (i.e., the northeast and northwest regions of the disk around DM Tau and V4046 Sgr, respectively, are closer to the observer), while the disk around MWC 480 has a negative inclination (i.e., the southwest region of the disk is closer to the observer than the northeast).

Here we show the two-dimensional and one-dimensional posterior distribution functions for the fiducial MCMC trials for DM Tau (Figure C2), MWC 480 (Figure C3), and V4046 Sgr (Figure C4). Figure C5 shows the chains for DM Tau trials, demonstrating that the walkers quickly converge toward a preferred solution.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size