A Three-dimensional View of Turbulence: Constraints on Turbulent Motions in the HD 163296 Protoplanetary Disk Using DCO⁺

The planet formation environment in disks around young stars is highly dynamic, subject to large azimuthal and radial velocities, coupled gas and dust dynamics, as well as perturbations from magnetic fields, embedded planets, and/or passing stars. While these phenomena are challenging to fully characterize, the Atacama Large Millimeter/submillimeter Array (ALMA), with its ability to spatially resolve Doppler motions from multiple molecular emission lines at high signal-to-noise ratios (S/Ns), is advancing our understanding of the protoplanetary disk kinematic structure. ALMA observations have used the dominant Keplerian motion to measure the mass of the central star(s) (Czekala et al. 2015, 2016), while deviations from Keplerian rotation have been investigated for warps (Casassus et al. 2015), fast radial streams (Rosenfeld et al. 2014; van der Plas et al. 2017), winds (Salyk et al. 2014), gas-pressure support, self-gravity (Rosenfeld et al. 2013), and even protoplanetary candidates (Factor et al. 2017).

The amplitude and nature of the motion within the disk represent important parameters in the planet formation process. In particular, non-thermal motions can influence processes ranging from the collisional growth of small dust grains (Testi et al. 2014) to the ability of massive planets to open a gap (e.g., Fung et al. 2014). Because of these effects, constraints on the strength of turbulence, especially near the disk midplane where planet formation is likely to occur, are crucial for understanding the evolution of planets and planetary systems.

Constraints on the vertical structure of turbulence are also valuable for understanding the physical mechanism driving turbulence. The magnetorotational instability (MRI Balbus & Hawley 1991, 1998), which in protoplanetary disks relies on the coupling of magnetic fields to weakly ionized gas in a rotating disk, generates turbulence that can increase from a few percent of the local sound speed near the midplane up to nearly the sound speed in the upper layers (Fromang & Nelson 2006; Simon et al. 2013, 2015). Gravitoturbulence, which relies on gravitational instabilities near the midplane, also predicts large non-thermal motions, but with a nearly constant vertical profile (Forgan et al. 2012; Shi & Chiang 2014). The vertical shear instability may also produce turbulence at the level of a few percent of the sound speed (Nelson et al. 2013), although the vertical structure of the turbulence arising from this instability has not been fully characterized in detail.

Previous radio observations using the SMA and PdBI have found tentative evidence for turbulence (e.g., Hughes et al. 2011; Guilloteau et al. 2012). In our analysis of the ALMA Science Verification data of the disk around HD 163296 (Flaherty et al. 2015), we placed a tight constraint on the turbulence in the upper layers (v_turb < 0.04c_s) using CO(3-2) emission, but the modest S/N of the optically thinner ¹³CO(2-1) and C¹⁸O(2-1) emissions prevented us from tightly constraining the motion near the midplane. Teague et al. (2016) use CO(2-1), CN(2-1), and CS(2-1) to measure the turbulence in the disk around TW Hya, finding non-thermal motion from the CO(2-1) emission layer, although the uncertainty in the amplitude calibration severely limits the measurements.

We take advantage of the complex chemical environment within protoplanetary disks to probe the vertical structure of the turbulence. CO is the second most abundant molecule after H₂ and is highly optically thick. The emission from less abundant isotopologues, such as ¹³CO or C¹⁸O, originates deeper in the disk, but the ability of these isotopologues to trace the midplane is still limited by the freeze-out of CO gas onto dust grains (e.g., Qi et al. 2015). Other molecules may remain abundant in the gas phase at colder temperatures. DCO⁺ emission, which has been observed in a number of protoplanetary disk systems (van Dishoeck et al. 2003; Guilloteau et al. 2006; Qi et al. 2008; Öberg et al. 2010, 2011; Mathews et al. 2013; Teague et al. 2015; Qi et al. 2015; Yen et al. 2016; Huang et al. 2017), is predicted to arise from close to the cold midplane (Willacy 2007; Öberg et al. 2012, 2015; Teague et al. 2015). Recent observations have found that DCO⁺ may exist in regions of the disk where CO was previously believed to be frozen out (Öberg et al. 2015). This spatial distribution makes it a particularly useful tool for measuring gas kinematics in a way that complements the much brighter CO emission.

Here we present an analysis of ALMA observations of DCO⁺(3-2), C¹⁸O(2-1), and CO(2-1) from the disk around the Herbig star HD 163296. Flaherty et al. (2015) previously found v_turb < 0.04c_s in the upper layers of this system based on the ALMA Science Verification CO(3-2). With these new high S/N observations we can probe similar levels of non-thermal motion for a region closer to the midplane, where planet formation is expected to occur. In Section 2 we present our data; in Section 3 we lay out an initial analysis of the structure, which informs our detailed model (Section 4) and the subsequent results (Section 5). In Section 6 we discuss the implications of the weak turbulence seen in all three tracers.

Cycle 2 observations (2013.1.00366.S) with the Atacama Large Millimeter/submillimeter Array of HD 163296 were taken over the course of five nights (2014 June 4, 14, 16, 17, and 29). During each night, the on-source integration time was 52 minutes, and the nearby quasar J1733-1304 was used for phase and gain calibration. Baselines ranged from 10kλ to 480kλ, corresponding to a spatial resolution of ∼05.

Four spectral windows were observed: spectral window 0 was centered on CO(2-1), spectral window 2 was centered on DCO⁺(3-2), spectral window 3 was centered on C¹⁸O(2-1), and spectral window 1 was a continuum window centered on 232.71 GHz. The spectral windows centered on line emission had 3840 channels, each 0.015 MHz (20 m s⁻¹) wide. The continuum spectral window had a total bandwidth of 2 GHz covered by 128 channels. Data were self-calibrated using the continuum spectral window.

The spectral line data were continuum-subtracted and binned by a factor of 10, to 0.2 km s⁻¹ wide channels, to increase the S/N and decrease the time taken to model the emission. We do not anticipate that this spectral binning will substantially limit our ability to constrain turbulence; we found similar constraints on turbulence based on CO(3-2) emission when using either a 0.1 km s⁻¹ resolution spectrum (<0.04c_s) or a 0.3 km s⁻¹ spectrum (<0.06c_s; Flaherty et al. 2015). Spatial resolution and S/N play a larger role in our ability to tightly constrain turbulence in the outer disk. Uncertainties on the visibilities were derived based on the dispersion around each baseline in line-free channels. Images were generated using robust weighting for the spectral lines and uniform weighting for the continuum window, resulting in beam sizes of ∼05 and ∼04, respectively. A clean mask based on the emission pattern of a Keplerian disk was applied to the CO(2-1) data. Additionally, we extracted DCO⁺(5-4) and Band 7 continuum data from the ALMA Science Verification observations of HD 163296, with self-calibration applied based on the HD 163296 Band 7 Science Verification CASA guide.¹¹

Derived quantities for the cycle 2 observations, which reach peak S/Ns of 50–1000, are listed in Table 1. Rosenfeld et al. (2013) measure CO(2-1) and C¹⁸O(2-1) fluxes of 46 ± 5 Jy km s⁻¹ and 5.8 ± 0.6 Jy km s⁻¹, respectively, from the ALMA science verification data. Qi et al. (2011) report fluxes for these same lines of 54.17 ± 0.39 Jy km s⁻¹ and 6.40 ± 0.16 Jy km s⁻¹ based on SMA observations. Qi et al. (2015) and Boneberg et al. (2016) report C¹⁸O(2-1) integrated fluxes from the ALMA Science Verification data of 6.8 ± 0.7 Jy km s⁻¹ and 6.2 ± 0.4 J km s⁻¹. SMA observations of the 1.3 mm continuum by Isella et al. (2007) measure a flux of 705 ± 12 mJy. ALMA observations of DCO⁺(3-2) have been reported by Yen et al. (2016) and Huang et al. (2017), with Huang et al. (2017) deriving an integrated line flux of 1.29 ± 0.04 Jy km s⁻¹. Our integrated line and continuum intensities are consistent with these previous measurements when accounting for uncertainties in calibration, which are much larger than the statistical uncertainties. The flux of the amplitude calibrator source varied from 1.211 ± 0.001 Jy to 1.313 ± 0.005 Jy over the course of our observations, while the continuum and individual lines varied by 5%–10% from night to night, suggesting larger errors due to night-to-night calibration. An additional uncertainty arises from the imprecise knowledge of the flux of the calibrator source; the ALMA calibrator database lists a 15% uncertainty on the flux at 233 GHz of J1733. Throughout our analysis, we assume a 20% uncertainty on the amplitude calibration of our data.

Table 1.  Data Summary

Line	Channel Width	Beam Size	Beam PA	rms	Integrated Intensity	Peak Flux
		(FWHM)		(mJy/beam)		(Jy/beam)
CO(2-1)	0.20 km s−1a	058 × 047b	756	6.5	44.1 ± 0.3 Jy km s−1	0.53
C18O(2-1)	0.21 km s−1a	058 × 049b	−888	1.4	5.37 ± 0.06 Jy km s−1	0.13
DCO+(3-2)	0.21 km s−1a	059 × 050b	896	1.1	1.09 ± 0.01 Jy km s−1	0.06
continuum	2 GHz	042 × 035c	−891	0.1	626 ± 19 mJy	0.13

Notes.

^aBinned down by a factor of 10 from the original spectral resolution. ^bDerived with robust = 0.5. ^cDerived with uniform weighting.

Download table as:  ASCIITypeset image

While previous observations of DCO⁺ in the disk around HD 163296 found a broad ring (Mathews et al. 2013; Qi et al. 2015), the high S/N and spatial resolution of our data are able to reveal new details about this structure (Figure 1). Emission is clearly detected out to ∼250 au (Figure 2) with peak S/N ∼ 50. The channels near the line peak reveal additional structure. Rather than smooth emission, we see a peanut-shaped feature, which is suggestive of the presence of two distinct rings of DCO⁺. Recent observations have begun to reveal gas rings around other stars in both DCO⁺ (Öberg et al. 2015) and CO (Huang et al. 2016; Schwarz et al. 2016), and such behavior here is not surprising, especially in light of the multi-ringed structure recently detected in the dust continuum emission at millimeter wavelengths (Isella et al. 2016; Zhang et al. 2016).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To obtain an initial estimate of the structure, we employ an image deconvolution routine modeled after Lucy (1974) and Scoville et al. (1983). This procedure utilizes concentric rings of emission moving in Keplerian orbits to match the full three-dimensional (position–position–velocity) data set. By iteratively varying the emissivity of each ring independently until it matches the emission profile, while accounting for the finite spectral and spatial resolution, this process can fit an arbitrary radial distribution of gas. The known velocity profile, Keplerian to first order, allows for a deconvolution of structure even below the nominal spatial resolution (∼05 = 60 au).

When applied to the DCO⁺(3-2) data, this routine finds three distinct rings of emission at ∼60, ∼140, and ∼260 au (Figure 3). The two outer rings were evident based on the appearance of the channel maps, while the third ring appears as enhanced emission within the line wings. The overall distribution of DCO⁺, with emission covering ∼60 au to ∼250 au, is consistent with that found from previous studies. Qi et al. (2015) detect DCO⁺(4-3) emission from 40 au to 250 au, while Mathews et al. (2013) find DCO⁺(5-4) in a broad ring centered at roughly 90 au. Our high S/N data is able to resolve the broad features seen in these previous data into its underlying constituents.

Zoom In Zoom Out Reset image size

Deconvolution reveals details about the radial surface brightness structure but does not constrain the density or temperature within the rings. The previous measurements of DCO⁺(4-3) (Qi et al. 2015) and DCO⁺(5-4) (Mathews et al. 2013) emission can be combined with our DCO⁺(3-2) measurement to obtain an initial estimate of the excitation temperature of the DCO⁺ gas. The excitation temperature defined by the J = (u₁–l₁) and J = (u₂–l₂) transitions (u₁ > u₂) is

$$\begin{eqnarray}&&{T}_{\mathrm{ex}}=\displaystyle \frac{h{\nu }_{{u}_{1}{u}_{2}}}{k\mathrm{ln}\left(\tfrac{{S}_{{u}_{1}{l}_{1}}{g}_{{u}_{2}}{A}_{{u}_{2}{l}_{2}}{\nu }_{{u}_{2}{l}_{2}}}{{S}_{{u}_{2}{l}_{2}}{g}_{{u}_{1}}{A}_{{u}_{1}{l}_{1}}{\nu }_{{u}_{1}{l}_{1}}}\right)},\end{eqnarray} \tag{ 1 }$$

where ${\nu }_{{u}_{1}{u}_{2}}$ is the frequency of the transition between levels u₁ and u₂, ${S}_{{u}_{1}{l}_{1}}$ is the integrated flux density between the u₁ and l₁ levels, g is the statistical weight, and A is the Einstein A value. Including a 20% uncertainty on the amplitude calibration for each data set, we find a temperature of 16 ± 4 K and 11 ± 4 K when comparing DCO⁺(3-2) to DCO⁺(5-4) and DCO⁺(4-3), respectively. These low temperatures suggest that DCO⁺ arises from the cold temperature midplane rather than from the warm surface layers. This calculation assumes LTE, which is likely applicable for DCO⁺. The DCO⁺(3-2) and DCO⁺(5-4) lines have critical densities of 1.9 × 10⁶ cm⁻³ and 1.4 × 10⁷ cm⁻³, respectively. According to the density and temperature structure derived in Flaherty et al. (2015), typical densities of ∼10⁹ cm⁻³ are found near the midplane. Densities below 10⁶ cm⁻³ are only seen in the upper reaches of the outer disk (z > 80 au for r > 200 au), an area that does not contribute significantly to the total emission. The excitation temperature derived from the line intensity ratios is an intensity-weighted disk average temperature. In a realistic disk, the temperature decreases with distance from the star, and based on our previous analysis (Flaherty et al. 2015), we expect the temperature to vary by ∼30% between the inner and outer rings, which is comparable to the uncertainty in our temperature estimate.

In order to constrain the turbulence from the line emission, we employ a parametric disk structure and ray-tracing code to generate model visibilities that can be compared directly to the data. The modeling code comes from Flaherty et al. (2015), which is based on earlier works by Rosenfeld et al. (2013) and Dartois et al. (2003). We provide a brief summary here with more detail available in Flaherty et al. (2015).

The underlying surface density structure is assumed to follow the functional form predicted 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{ 2 }$$

Here M_gas is the total gas mass, R_c is the critical radius, and γ controls the power-law shape of the surface density for r≪R_c. The temperature is assumed to follow a power law with radius, with a vertical gradient connecting the cold midplane (T_mid) with the warm atmosphere (T_atm):

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

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

$$\begin{eqnarray}&&{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{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&{Z}_{q}=70\ \mathrm{au}{(r/150\mathrm{au})}^{1.3}.\end{eqnarray} \tag{ 6 }$$

The hydrostatic equilibrium calculation is performed to derive the volume density given the temperature and surface density structure. The velocity field is Keplerian, with corrections for the height above the midplane and the pressure support of the gas, as in Rosenfeld et al. (2013). The line is assumed to be a Gaussian with width

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

We included a non-thermal line-broadening term, v_turb, which is assumed to be proportional to the local isothermal sound speed ($\sqrt{2{k}_{B}T(r,z)/\mu {m}_{H}}$). We interpret this additional broadening as due to turbulence, although other deviations from Keplerian and thermal motion may contribute to this term.

Our main goal is to use the individual lines to constrain the turbulent structure at the various vertical emission layers of these molecules within the disk. As such, we start with similar underlying structures, but allow for variations between the molecules to best fit the data. Although there is certainly a single temperature/density/kinematic structure that applies to all of the emission from this system, the assumptions needed to fit an individual line (e.g., assuming a temperature structure to break its degeneracy with the surface density for an optically thin line), and our incomplete knowledge of the complex physics/chemistry in protoplanetary disks (e.g., selective photodissociation, CO non-thermal desorption), can lead to differences in the model parameters derived from different lines. The constraints on the model parameters simply reflect the values needed, in the context of our assumed parametric model, to reproduce the flux from the emitting region of a particular line.

With this framework in mind, we consider the detailed structure for each molecule in turn. CO is confined to a molecular layer, whose upper boundary is set by photodissociation and whose lower boundary is set by freeze-out. The photodissociation boundary is defined by the height at which the vertical column density reaches Σ₂₁ = 0.79, where Σ₂₁ is the surface density in units of 1.59 × 10²¹ cm⁻². Freeze-out occurs below a temperature of 19 K, and for CO, which is less sensitive to midplane abundance, we assume that freeze-out leads to a drop in CO abundance by eight orders of magnitude. For C¹⁸O(2-1), which is more sensitive to the abundance at the cold midplane, we assume an abundance decrement of a factor of five, consistent with recent studies (Qi et al. 2015). We assume isotope abundances of ¹⁸O/¹⁶O = 557 and ¹³C/¹²C = 69 (Wilson 1999). The abundance of CO relative to H₂, X_CO, is allowed to vary when modeling the optically thin C¹⁸O, and is fixed at 10⁻⁴ when modeling optically thick CO(2-1). We assume LTE for all models which, as discussed earlier, is a reasonable assumption even for DCO⁺.

Since DCO⁺ is optically thin, it is difficult to disentangle the density and temperature, and we choose to fix the temperature structure while allowing the density to vary. We use the temperature parameters derived in Flaherty et al. (2015) from the CO(3-2) emission (q = −0.216, R_c = 194 au, T_atm0 = 94 K, γ = 1, T_mid0 = 17.5 K). The uncertainty on these parameters will feed into the uncertainties on the density and turbulence, an effect that we consider in more detail below. The DCO⁺ molecules are assumed to be distributed in three rings whose radial abundance profile is given by

$$\begin{eqnarray}\left[\displaystyle \frac{{\mathrm{DCO}}^{+}}{{{\rm{H}}}_{2}}\right] & = & {\left[\displaystyle \frac{{\mathrm{DCO}}^{+}}{{{\rm{H}}}_{2}}\right]}_{\mathrm{in}}\exp \left(-\displaystyle \frac{{(r-{R}_{\mathrm{in}})}^{2}}{{\sigma }_{\mathrm{in}}^{2}}\right)\\ & & +{\left[\displaystyle \frac{{\mathrm{DCO}}^{+}}{{{\rm{H}}}_{2}}\right]}_{\mathrm{mid}}\exp \left(-\displaystyle \frac{{(r-{R}_{\mathrm{mid}})}^{2}}{{\sigma }_{\mathrm{mid}}^{2}}\right)\\ & & +{\left[\displaystyle \frac{{\mathrm{DCO}}^{+}}{{{\rm{H}}}_{2}}\right]}_{\mathrm{out}}\exp \left(-\displaystyle \frac{{(r-{R}_{\mathrm{out}})}^{2}}{{\sigma }_{\mathrm{out}}^{2}}\right),\end{eqnarray} \tag{ 8 }$$

where R_in, R_mid, and R_out are the central locations of the three rings and [DCO⁺/H₂]_in, [DCO⁺/H₂]_mid, and [DCO⁺/H₂]_out are the peak abundances of DCO⁺ in the three rings. We assume σ_in, σ_mid, and σ_out are equal to one-tenth of their respective central radii; we examine the consequences of this assumption later. The abundance is assumed to be constant between the vertical column densities Σ₂₁ = 0.79, 1000. Due to the low excitation temperature derived above, the complex radial structure with emission at radii where CO (a key molecule in the creation of DCO⁺) is expected be frozen out, and the evidence for the modest depletion of CO gas during freeze-out (Qi et al. 2015), we exclude the effects of freeze-out from our DCO⁺ model. Including DCO⁺ freeze-out for T < 19 K in an additional model fit, we find that only [DCO⁺/H₂]_mid and [DCO⁺/H₂]_out move from their fiducial values, increasing by 0.1 and 0.3 dex, respectively. For our fiducial DCO⁺ model, the free parameters are R_in, R_mid, R_out, [DCO⁺/H₂]_in, [DCO⁺/H₂]_mid, [DCO⁺/H₂]_out, as well as v_turb and the inclination.

Our final CO(2-1) models have free parameters q, R_c, T_atm0, T_mid0, R_in, v_turb, and inclination, while for C¹⁸O(2-1) we fix R_c, T_atm0, and T_mid0 and add X_CO and γ as free parameters. Allowing the temperature to vary when fitting CO(2-1), while fixing it for C¹⁸O(2-1) and DCO⁺ opens up the possibility of different temperature structures for the different molecules. After fitting, we confirm that the differences in temperature are small, especially compared to the underlying systematic uncertainty in amplitude calibration. After finding approximately accurate parameters, we refine the center of the disk relative to the phase center as well as the velocity center of the line using a simple grid search. We find these to be [−003, 002] and 5.74 km s⁻¹ consistently between all of the lines. The position angle is fixed (P.A. = 312°) but inclination is allowed to vary between the three molecules since its measurement will depend in part on the vertical location of the emission within the disk. As with the temperature, we find consistent results for all three molecules. We also assume a distance (d = 122 pc), stellar mass (M_* = 2.3 M_⊙; Montesinos et al. 2009), and disk mass (M_gas = 0.09 M_⊙). While not included in our modeling, the ∼10% uncertainty on the distance to HD 163296 (van Leeuwen 2007) will directly translate into a ∼10% systematic uncertainty on any radial distances measured in this system. The disk mass was taken from Isella et al. (2007), who used continuum observations of submillimeter dust emission to derive a dust mass and converted it to a gas mass assuming a gas-to-dust mass ratio of 100. Recent work has brought into question standard assumptions about the gas-to-dust mass ratio, CO/H₂, and C¹⁸O/CO (Bergin et al. 2013; Miotello et al. 2014; Ansdell et al. 2016; Kama et al. 2016), raising the possibility of a significant deviation from our assumed disk mass. To accommodate this, we allow X_CO to vary, but note that any deviations from the fiducial value of 10⁻⁴ may reflect the fact that the assumed dust-to-gas ratio or ¹⁸O/¹⁶O abundance is incorrect.

Once a suite of parameters has been specified, model images are generated and converted to visibilities using the MIRIAD task UVMODEL. The model visibilities are then directly compared to the data, with the goodness of fit calculated using the chi-squared statistic. To derive the posterior distributions for each parameter, we employ the MCMC routine EMCEE (Foreman-Mackey et al. 2013) based on the affine-invariant algorithm originally proposed in Goodman & Weare (2010). Our chains typically consist of 80 walkers and 1000 steps, with the first 500 removed as burn-in. Simple, physical priors are employed (e.g., ${v}_{\mathrm{turb}}\gt 0$), and we employ linear spacing in all parameters except R_c, where we fit $\mathrm{log}({R}_{c})$.

Systematic uncertainty in the amplitude calibration of the data is an important source of error since it dominates the statistical noise in these high S/N data. Since systematic uncertainty can directly affect the temperature derived from optically thick lines, and hence the thermal broadening, it potentially has a large effect on a derived non-thermal broadening (Teague et al. 2016), although high spatial resolution imaging can help to break this degeneracy (Simon et al. 2015). To account for systematic uncertainty in our MCMC modeling framework, we introduce an additional parameter, sys, that is a multiplicative scale factor on the amplitude of the model emission. The model is scaled down by sys to force the other parameters (e.g., temperature, abundance, turbulence) to adjust upwards accordingly. A value of sys = 0.8 corresponds to a true disk flux 20% lower than is derived from our data, while sys = 1.2 simulates a true disk flux 20% higher than what we derived from our data. By varying sys between 0.8 and 1.2, we can sample the full size of the posterior distribution functions allowed by the systematic uncertainty.

4.1. The Influence of Model Assumptions

5.1. DCO⁺

Using our ray-tracing disk model and MCMC-fitting routine, we are able to accurately confine the central radii of the three rings to 65.7 ± 0.9 au, ${149.9}_{-0.7}^{+0.5}\,\mathrm{au}$, and 259 ± 1 au (Table 2). Posterior distribution functions (PDFs), generated using the Python module corner (Foreman-Mackey et al. 2014), are shown in Figure 4; each parameter shows a Gaussian PDF with the exception of turbulence, which is an upper limit. The model spectrum is very well matched to the data, and a moment 0 (total-intensity) map of the residuals (subtracted in the visibility domain and imaged using the clean algorithm) only shows small 3σ features (Figure 5). This model is also able to match much of the DCO⁺(5-4) emission, despite the fact that this line was not included in the fit.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 2.  DCO⁺(3-2) Model Results

Model	vturb/cs	Inclination	Rin	Rmid	Rout	log([DCO+/H2]in)	log([DCO+/H2]mid)	log([DCO+/H2]out)	${\chi }_{\nu }^{2}$ a
		(◦)	(au)	(au)	(au)
Fiducial	<0.04	${48.8}_{-0.1}^{+0.2}$	65.7 ± 0.9	${149.9}_{-0.7}^{+0.5}$	259 ± 1	−10.83 ± 0.03	−10.79 ± 0.01	−10.99 ± 0.01	1.1076
Narrow Rings	<0.6	48.7 ± 0.2	77.9 ± 0.7	151.7 ± 0.3	${252.7}_{-0.7}^{+0.9}$	$-{9.43}_{-0.04}^{+0.05}$	$-{9.64}_{-0.02}^{+0.07}$	$-{10.02}_{-0.01}^{+0.03}$	1.1084
Variable ring widthb	<0.06	${48.6}_{-0.2}^{+0.1}$	${60.6}_{-1.1}^{+0.8}$	${148.0}_{-1.3}^{+1.4}$	${266}_{-2}^{+1}$	−10.86 ± 0.04	$-{11.06}_{-0.02}^{+0.04}$	−11.15 ± 0.02	1.1074
sys = 1.2	<0.05	${48.9}_{-0.02}^{+0.1}$	66.9 ± 1.1	149.8±0.7	${259.7}_{-0.9}^{+1.0}$	−10.62 ± 0.04	−10.62 ± 0.01	−10.89 ± 0.01	1.1075
sys = 0.8	<0.04	${48.7}_{-0.2}^{+0.1}$	${65.8}_{-1.3}^{+1.5}$	149.8 ± 0.7	259 ± 1	−11.02 ± 0.03	−10.95 ± 0.01	−11.11 ± 0.01	1.1076
Tmid0 variablec	<0.12	${48.8}_{-0.2}^{+0.1}$	${65.9}_{-0.9}^{+1.0}$	149.4 ± 0.3	260.6 ± 1.1	−10.30 ± 0.07	$-{10.44}_{-0.04}^{+0.06}$	$-{10.85}_{-0.02}^{+0.03}$	1.1075

Notes. Median, plus 3σ ranges, of the posterior distribution functions derived for each model. No significant degeneracies were found between the listed parameters.

^aDCO⁺(3-2) reduced chi-squared for the model defined by the median of the PDFs. Given the large number of degrees of freedom in the data, the differences in reduced chi-squared between the different models are highly significant. ^bThe radial widths of the middle and outer rings are no longer fixed at 10% of the radial location, but are instead allowed to vary as free parameters. We derive FWHMs of the middle and outer rings of 60 ± 3 au and 94 ± 7 au, respectively. ^cT_mid0 is allowed to vary, with DCO⁺(5-4) included in the fitting. T_mid0 is constrained to ${12.2}_{-0.5}^{+0.4}\,{\rm{K}}$.

Download table as:  ASCIITypeset image

Mathews et al. (2013) found that the DCO⁺(5-4) emission was confined to a ring from 110 to 160 au, consistent with the middle ring in our models. Qi et al. (2015) fit DCO⁺(4-3) with a ring of emission stretching from ${40}_{-3}^{+6}\,\mathrm{au}$ to ${290}_{-8}^{+6}\,\mathrm{au}$, similar to the total radial coverage of the three rings. Our derived radial structure confirms the finding by Qi et al. (2015) that DCO⁺ does not trace the CO condensation front, which is located at 90 au based on N₂H⁺. The CO snow line is bracketed by the two innermost rings, with neither contributing significant emission at 90 au.

In addition to the location of the rings, we place strong constraints on the other model parameters (v_turb, inclination and abundance). We measure peak abundances within the three rings of log([DCO⁺/H₂]) = −10.83 ± 0.03, −10.79 ± 0.01, and −10.99 ± 0.01 for the inner, middle, and outer rings, respectively. These values are similar to those predicted by chemical models (Willacy 2007; Öberg et al. 2015; Teague et al. 2015). We are able to reproduce the data almost entirely with Keplerian and thermal motion, with any non-thermal velocity dispersion limited to <0.04c_s. In the discussion below, we interpret this as a limit on the turbulence, but it also applies to other non-thermal non-Keplerian effects, such as spiral arms, disk winds, or warps. These processes do not contribute significantly in the DCO⁺ emitting region of the disk around HD 163296.

The statistical errors on the derived parameters are small (<1%) due to the large S/N of these data. The largest source of error is likely systematic effects that are not accounted for in our MCMC modeling. Here we consider three sources of uncertainty: (1) the assumption about the ring widths, (2) the uncertainty in the underlying temperature structure, and (3) the 20% systematic uncertainty in the amplitude calibration.

Narrow versus Wide Rings. In modeling the rings, we assume that each is a Gaussian with a radial width, σ, equal to 10% of the rings' radial location. If our assumption of σ is incorrect, then other model parameters, such as turbulence, may be driven away from their true values in order to compensate for the incorrect assumption. To test the possibility of much narrower rings, we rerun our MCMC trial with the widths set at 1% of the rings' radii. We find that the fit is significantly worse, a >10σ difference according to the Akaike Information Criterion (Akaike 1974), with the narrow ring model severely underpredicting the emission at line peak and overpredicting the emission at line center (Figure 6). While the reduced chi-squared for these two models are similar (Table 2), the large number of degrees of freedom in the data makes this difference highly significant. It appears that the walkers have attempted to accommodate the narrow rings with increased turbulence (${v}_{\mathrm{turb}}\lt 0.6{c}_{s}$; Table 2), but the change in turbulence cannot make up for the narrow rings. Based on these results, we can assume that we have not substantially overestimated the radial width of the DCO⁺ rings. This conclusion is consistent with the FWHM of the rings derived using the Richardson–Lucy algorithm: ∼40 au, ∼70 au, and ∼100 au compared to our assumption of ∼20 au, ∼30 au, and ∼60 au for the inner, middle, and outer rings respectively. When allowing the radial width of two outer rings to vary, we measure an FWHM of 60 ± 3 au and 94 ± 7 au for the middle and outer rings respectively, with a modest change in our turbulence upper limit (Table 2). This indicates that assuming a radial width for the rings likely does not substantially bias our estimate of the turbulence.

Zoom In Zoom Out Reset image size

Temperature Uncertainty. In our fiducial model, we assumed an underlying temperature structure, but this structure is not known perfectly and its uncertainty will propagate into other quantities. In particular, the midplane temperature (T_mid0) will strongly influence the derived abundances for this optically thin emission. Our knowledge of the midplane temperature is constrained in part by the CO(3-2) emission (T_mid0 = 17.5 ± 0.25 K; Flaherty et al. 2015) as well as by the ratio of the DCO⁺(3-2) emission to the DCO⁺(5-4) emission (16 ± 4 K). To account for the uncertainty on T_mid0, we run an additional MCMC trial with T_mid0 as a free parameter. This trial models both DCO⁺ emission lines and includes a Gaussian prior on T_mid0 based on the previous CO(3-2) result.

Despite the prior on ${T}_{\mathrm{mid}0}$, we find that the walkers have converged on ${T}_{\mathrm{mid}0}={12.2}_{-0.5}^{+0.4}\,{\rm{K}}$ (Table 2). With this lower temperature comes an increase in the turbulence limit to v_turb < 0.12c_s and an increase in the rings' abundances by 0.1–0.5 dex. This lower temperature model is as effective at fitting the DCO⁺(3-2) model as the fiducial model and is a good match to the DCO⁺(5-4) spectrum (Figure 7). As expected, there is a strong degeneracy between the midplane temperature and the abundance of the three rings (Figure 8), with the anticorrelation being strongest for the two outer rings. After the walkers have converged around the best fit, there is no strong degeneracy between the midplane temperature and turbulence, but the rise in the turbulence limit relative to the fiducial model indicates that there is some anticorrelation between the two quantities. Even when accounting for this degeneracy, we find that the turbulence PDF still shows signs that it is an upper limit (with a tail that continues to zero turbulence).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The midplane temperature derived from the DCO⁺ modeling is significantly smaller than that derived from CO(3-2), and there are a number of factors that could contribute to this discrepancy. The DCO⁺ measurement of T_mid0 relies on the relative fluxes of the J = (3-2) and J = (5-4) lines, and as a result will be strongly influenced by any uncertainty in the amplitude calibration for these data. The CO(3-2) measurement relies on the geometric separation of the near and far side of the disk, which is less severely influenced by amplitude calibration, although it is a more indirect measure of T_mid0 since little to no CO(3-2) emission arises from the cold midplane. Deviations of the midplane temperature from our simple power-law prescription may also influence our measurements. More detailed modeling is needed to understand the temperature structure and to reconcile the different measurements.

Systematic Uncertainty. As discussed above, the uncertainty in the amplitude calibration of radio interferometric data is larger than the statistical uncertainties in these high S/N data sets. Here we add in the systematic uncertainty parameter using the method described above and allow it to vary freely between 0.8 and 1.2. We find that sys very quickly moves toward the upper bound of 1.2, restricting itself to a narrow range of parameter space. During this move, the turbulence, inclination, and location of the three rings do not vary from the fiducial model (Table 2), indicating that they are not strongly degenerate with the amplitude calibration. The peak abundances all increase by 0.2 dex, since the emission is mostly optically thin and its total flux is strongly correlated with abundance, assuming a constant temperature.

To explore the low-brightness case, we fix sys at 0.8 and rerun the fit. We find that the turbulence, inclination, and location of the three rings again do not vary from the fiducial model, while the abundances are 0.2 dex lower (Table 2). Based on this analysis, we can ascribe a 0.2 dex uncertainty on the abundance for a 20% uncertainty on the flux calibration, assuming fixed temperature.

These two trials indicate that turbulence is not strongly degenerate with the systematic uncertainty. When sys = 1.2, where the model fits the data as though the disk is 20% brighter than recorded, the limit on v_turb increases to only 0.05c_s. Conversely, when sys = 0.8, we recover the same limit as in the fiducial model (v_turb < 0.04c_s). This lack of degeneracy arises because a change in amplitude calibration results in a uniform change in the surface brightness of the emission in each channel, while turbulence mainly changes the total surface area of the emission in each channel (Simon et al. 2015). The high spatial resolution is key to breaking the degeneracy between the turbulence and the amplitude calibration.

This effect is demonstrated in Figure 9. Here we show the intensity profile along a slice through the central velocity channel for the best-fit model (with zero turbulence) and three comparison models. The comparison models have individual parameters adjusted such that the peak flux increases by 20%, turbulence is increased to 0.4 c_s, DCO⁺ abundances are increased by 0.2 dex, or the midplane temperature is increased to 23 K. Only looking at the total flux, these three models would be indistinguishable from calibration uncertainty, but the spatial information breaks some of this degeneracy. In particular, the model with high turbulence substantially broadens the intensity profile. This difference, while small, is statistically significant after integrating over the entire 3D data set. Variations in the abundances do not significantly change the profile shape, explaining the strong degeneracy among these parameters and the systematic uncertainty. By modeling the entire spatial distribution of the emission, we are able to leverage the fact that neighboring pixels contain complementary information, e.g., pixels at similar distances from the central star will have similar temperatures and utilizing the information from all of these pixels provides a tighter constraint on temperature than is possible by modeling a single line of sight. This approach does require assumptions to be made about the underlying structure (e.g., the midplane temperature is well described by a radial power law) but the lack of residuals suggests that these assumptions are appropriate for explaining the HD 163296 DCO⁺ observations.

Zoom In Zoom Out Reset image size

5.1.1. Source of the DCO⁺ Rings

The disk around HD 163296 is unique in having three DCO⁺ rings and the origin of the DCO⁺ at these locations depends on the creation of DCO⁺ from CO. One of the main chemical pathways for generating DCO⁺ is

$$\begin{eqnarray}&&{{\rm{H}}}_{2}{{\rm{D}}}^{+}+\mathrm{CO}\to {{\rm{H}}}_{2}+{\mathrm{DCO}}^{+}.\end{eqnarray} \tag{ 9 }$$

The efficiency of this pathway increases as the temperature decreases, indicating that DCO⁺ can be found in locations with cold CO. An additional warm pathway through CH₂D⁺ has recently been explored by Favre et al. (2015), but the low temperature derived from the DCO⁺(5-4)/(4-3)/(3-2) line ratios and from the DCO⁺(5-4) and (3-2) line fitting favors the low-temperature pathway.

Assuming complete CO freeze-out at low temperatures, the coldest CO will be found just inside the CO condensation front. Öberg et al. (2015) find that the inner DCO⁺ ring in the disk around IM Lup lies just inside the CO snow line. In the disk around HD 163296, we find that the innermost ring, at 65.7 ± 0.9 au, also lies inside the CO condensation front at 90 au (Qi et al. 2015).

The presence of DCO⁺ at larger radii, and at temperatures below 19–25 K, requires a return of CO to the gas phase. One way for this to occur is due to the inward radial migration of dust (Cleeves 2016). As dust grains migrate inward, the temperature of the outer disk rises due to the removal of opacity that would otherwise shield the midplane, and this increase in temperature may be enough to push the CO above its freeze-out temperature. We find that the outermost DCO⁺ ring, at 259 ± 1 au, lies just outside the edge of the continuum emission (250 au Isella et al. 2016), consistent with this model. To rise above the CO freeze-out temperature at this radius would require a ∼20% increase in temperature above that predicted by our CO(2-1) model fit, consistent with the 10%−30% increase in temperature predicted by Cleeves (2016). Huang et al. (2017) also find DCO⁺ emission beyond the dust outer edge around V4046 Sgr, similar to that observed here.

The origin of the middle DCO⁺ ring is less clear. Öberg et al. (2015) account for the outer ring around IM Lup using non-thermal desorption of CO by UV photons and cosmic rays penetrating to the midplane, and a similar process may be at work in the disk around HD 163296. Additional DCO⁺ in the warm molecular layer, created through the CH₂D⁺ pathway (Favre et al. 2015), may also contribute at these distances. The diversity of DCO⁺ morphologies observed by Huang et al. (2017) suggest that between different protoplanetary disks, different molecular pathways can dominate, and it is certainly possible that different pathways dominate at different radii within a single system.

These chemical processes can generate a single broad ring covering hundreds of astronomical units (Willacy 2007; Öberg et al. 2015; Teague et al. 2015), and it is possible that the two outer rings are in fact one broad ring with a central depression in the middle. The depression between the two outer DCO⁺ rings corresponds to a bright, and possibly optically thick, ring of dust (Isella et al. 2016). Optically thick dust will obscure our view of DCO⁺ emission arising from the midplane, creating the appearance of a gap in an otherwise smooth distribution of DCO⁺. Our modeling of the C¹⁸O(2-1) and dust emission, described below, as well as the observations of Isella et al. (2016), does not indicate that this dust ring is optically thick, although we cannot rule out a very narrow ring of optically thick dust.

Conversely, localized heating from a planet (Cleeves et al. 2015) could lead to a localized increase in emission at the middle ring, by returning CO to the gas phase in a region where it would otherwise be frozen out onto dust grains, independent of the chemical origin of the outer ring. The middle ring is aligned with one of the dark dust rings in the high-resolution observations of Isella et al. (2016) and if the continuum gap is caused by a large planet, localized heating by this planet could lead to an increase in DCO⁺ emission.

All of the scenarios for explaining the two outer rings involve the return of CO from the solid phase to the gas phase before creating DCO⁺. Recent studies have found rings of CO around TW Hya (Schwarz et al. 2016) and AS 209 (Huang et al. 2016) beyond the CO snow line, suggesting that CO desorption in the cold outer disk may be common. Searching directly for such an effect with CO is challenging given that the midplane is hidden behind the bright warm molecular layer, but there is evidence of cold CO near the midplane in the outer disk around HD 163296. Qi et al. (2015) find that fitting the C¹⁸O(2-1) emission requires a CO depletion factor of only 5 in the freeze-out region. In our modeling of C¹⁸O(2-1), discussed below, even after accounting for this small depletion factor, we find residuals at ∼250 au, consistent with enhanced CO abundance and/or temperature near the location of the outermost DCO⁺ ring.

5.2. C¹⁸O and CO

DCO⁺ indicates that there is weak turbulence near the midplane. By combining this limit with those from other molecular lines whose emission arises from above the midplane, we can constrain the vertical gradient of the turbulence in the outer disk. In our previous analysis of the Science Verification data (Flaherty et al. 2015), we found an upper limit on the turbulence of <0.04c_s from optically thick CO(3-2), but only <0.31c_s from CO(2-1), that traces the kinematics of the upper edge of the molecular layer. C¹⁸O(2-1), being more optically thin than CO(3-2) but still subject to freeze-out at the midplane, is sensitive to an intermediate layer between CO(3-2) and DCO⁺(3-2). Using the SV ALMA data, we were only able to limit turbulence to <0.4c_s, and here we use the higher S/N cycle 2 data to improve upon these limits.

While the DCO⁺ data is well fit with narrow rings, the C¹⁸O and CO emission appears much smoother; image deconvolution does not reveal any prominent rings of material. We initially proceeded with fitting these data with the smooth model described above but found that the best-fit models consistently underestimated the flux in the central ∼100 au in both lines by ∼7% (Figure 10).

Zoom In Zoom Out Reset image size

Observing an excess in both an optically thick and optically thin line suggests that it is the temperature structure that needs to be modified, rather than the CO or C¹⁸O abundance or the surface density profile. A similar excess was observed in the gas around TW Hya (Rosenfeld et al. 2012), possibly due to a change in temperature or a warp in the gas disk. Here we consider a change in temperature, although we cannot completely rule out other possible explanations.

The inner excess can be accommodated if q changes with radius; a steeper temperature profile in the inner disk increases the temperature, and hence the emission, within this region. Such behavior is expected for an externally illuminated disk. D'Alessio et al. (2006) find that the midplane temperature beyond ∼5 au varies as R^−0.25, while between 0.1 au and ∼5 au it varies as R^−0.75 (within 0.1 au viscous heating begins to play a large role). The change in temperature profile shape with radius is associated with the flaring of the surface layers, and since the midplane is heated by reirradiation from the upper disk layers, its temperature profile depends on the shape of the surface layers. In the inner disk, the scale height is relatively flat, resulting in a steep drop off in temperature with radius. In the outer disk, the flaring increases dramatically, leading to a shallower temperature profile. While the radius for the transition from the D'Alessio models is much smaller than the size of the inner excess seen here, the model calculations were performed using a T Tauri central star, and the radial scales are expected to be larger for a more luminous Herbig star.

This feature may have been implicitly included in previous efforts to model this system. Boneberg et al. (2016) use models of an illuminated disk to reproduce C¹⁸O observations, and find T ∝ R^−0.5 out to 90 au. They also note that their models, designed to fit the C¹⁸O(2-1) emission within 90 au, are a poor match to the outer disk flux, consistent with a change in the temperature structure between the inner and outer disks. Qi et al. (2015) use D'Alessio models to define the temperature and density structure, and hence include this behavior. De Gregorio-Monsalvo et al. (2013) are also able to reproduce much of the CO emission with an irradiated disk model. Our parametric models require the addition of this feature since we do not calculate the temperature structure based on the heating by the central star.

To accommodate this behavior, we explore models with a double power-law temperature structure. Within R_break, the radial exponent on the temperature is q_in while outside of this radius the radial exponent on the temperature is q. This adds two parameters (R_break, q_in) to our model. In modeling the DCO⁺ emission, we did not include this effect since it was not apparent in these data. Given the ∼100 au size scale of the inner excess, it likely only affects the innermost DCO⁺ ring. We can estimate an approximate scale for this effect based on our DCO⁺ modeling when T_mid0 was allowed to vary. We found that the DCO⁺ abundance of the inner ring increased by 0.5 dex when the midplane temperature decreased from 17.5 K to 12.2 K (Table 2). In fitting C¹⁸O(2-1) and CO(2-1) below, we find that the resulting models predict that the midplane temperature is 10–15 K higher at the location of the inner ring, which we anticipate would result in a decrease in the inner ring DCO⁺ abundance of 1–1.5 dex.

5.2.1. C¹⁸O(2-1)

As noted earlier, we model the C¹⁸O(2-1) emission by allowing γ, v_turb, inclination, X_CO, and R_in to vary, along with R_break and q_in. We fix the remaining surface density and temperature structure parameters based on the CO(3-2) fit (q = −0.216, R_c = 194 au, T_mid0 = 17.5 K, T_atm0 = 93.8 K). Based on the results of Qi et al. (2015), we assume a CO depletion of a factor of five when CO is frozen out of the gas phase, which is assumed to occur where T < 19 K.

The C¹⁸O emission shows evidence for a resolved inner hole, which we initially constrain to ${R}_{\mathrm{in}}={35.6}_{-2.0}^{+1.7}\,\mathrm{au}$. Such a tight constraint comes partly from the spatial and spectral resolution of the data, and partially from the dependence of the total flux in the inner disk on R_in. Given the high S/N of this data, this latter constraint dominates the final uncertainty, but is highly dependent on the parameters defining the shape of the inner excess (R_break and q_in) and the model grid resolution. For these reasons we fix R_in = 35.6 au for the final MCMC trial to aid in the convergence of the other parameters. The final results are listed in Table 3 and marginalized PDFs are shown in Figure 11. We find that we are able to successfully reproduce much of the emission, as seen in both the spectra and the moment map of the residuals (Figure 12).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 3.  C¹⁸O(2-1) Model Results

Parameter	Value
vturb	<0.05cs
γ	0.62 ± 0.01
incl	${48.26}_{-0.05}^{+0.07}$
Rina	${35.6}_{-2.0}^{+1.7}\,\mathrm{au}$
log([CO/H2])	−4.610 ± 0.003
Rbreak	${102}_{-7}^{+4}\,\mathrm{au}$
qin	$-{1.4}_{-0.1}^{+0.2}$
${\chi }_{\nu }^{2}$ b	0.992

Notes. Median, plus 3σ ranges, of the posterior distribution functions derived from fitting to the C¹⁸O(2-1) data.

^aDuring the final MCMC trial, this parameter was fixed at 35.6 au to facilitate the convergence of the other parameters. ^bReduced chi-squared for the model defined by the median of the posterior distribution functions.

Download table as:  ASCIITypeset image

With this emission line, we limit the turbulence to <0.05c_s, similar to DCO⁺. This limit falls below that of Flaherty et al. (2015) due to the much higher S/N, and higher spatial resolution, of these data. In the channel maps, we find positive residuals at radii >250 au, suggesting that we are underestimating either the CO abundance or the temperature at these large radii. Cleeves (2016) find that the inward migration of dust can lead to an increase in gas temperatures of 10%–30% in the outer disk. The dust emission is confined to within 250 au (Isella et al. 2016) and the increased C¹⁸O emission beyond the dust is consistent with additional heating returning CO to the gas phase.

Qi et al. (2015), in addition to deriving a modest depletion of CO in the outer disk, also find evidence that the CO snow line occurs at a temperature of 25 K. We run an MCMC trial with the CO snow line occurring at 25 K, instead of at 19 K in the fiducial model, and find that the CO abundance has significantly increased (${X}_{\mathrm{CO}}=-{4.446}_{-0.008}^{+0.007}$) and the inner temperature profile has become shallower (${q}_{\mathrm{in}}=-{0.92}_{-0.24}^{+0.09}$) while the other parameters are consistent with their values from the fiducial model. In particular, turbulence (v_turb < 0.06c_s) has not substantially varied, indicating that our assumption about the location of the CO snow line does not substantially affect our turbulence measurement.

An assumption in the use of continuum-subtracted emission in the fiducial model is that the dust is optically thin. Optically thick dust can shield C¹⁸O emission from escaping the disk, giving the appearance of a hole or gap; Öberg et al. (2015) find an inner hole in the C¹⁸O emission around IM Lup, which they ascribe to optically thick dust. We can characterize the influence of dust on the C¹⁸O(2-1) radiative transfer by including dust in our models and fitting a C¹⁸O(2-1) spectrum that has not been continuum-subtracted. We parameterize the dust as a radially varying dust-to-gas ratio, with an underlying structure based on the observations of Isella et al. (2016). The base dust structure is a power law in radius (dust-to-gas ratio $\sim {R}^{-{\gamma }_{\mathrm{dust}}}$) that extends from 10 to 250 au, with the normalization on the dust-to-gas ratio, defined at 120 au, allowed to vary along with the radial power-law exponent. We fix the three gaps at 66, 100, and 160 au with radial widths of 33, 26, and 55 au. The depletion of the dust-to-gas ratio within the gaps is treated as multiplicative constants (dtg1/13, dtg1/7, and dtg1/3.6 for the inner, middle, and outer rings, respectively, where dtg1 is a free variable). This allows the model to explore an area of parameter space where the bright rings between the gaps are optically thick while the gaps themselves can still be optically thin. We also assume an opacity of κ = 2.3 cm² g⁻¹ and that the dust-to-gas ratio is a constant as a function of height. We verify after the fact that the majority of the dust emission arises from close to the midplane (z/r < 0.1) and our assumption of vertically uniform dust, without accounting for settling, will not substantially bias our results. Since the dust emission becomes most important in the inner disk, its emission will be highly degenerate with the inner excess. To alleviate this, we fix q_in at −0.5, consistent with Boneberg et al. (2016), and R_in = 35.6 au.

In our resulting MCMC trial, we find reasonable fits to the data (Figure 13) with dust that is optically thin throughout the outer disk. The dust-to-gas ratio is ∼0.01 with a peak optical depth of ∼0.5 in the inner disk and τ < 1 in the bright rings. Our results are consistent with those of Isella et al. (2016), who find that the dust only becomes optically thick within ∼50 au, which may explain the 35.6 au C¹⁸O inner hole but is unlikely to affect the bulk of the emission that arises from larger radii. Complete modeling of the gas and dust is beyond the scope of the paper, but the optically thin nature of the dust supports our use of a continuum-subtracted line when modeling the C¹⁸O(2-1) emission.

Zoom In Zoom Out Reset image size

5.2.2. CO(2-1)

Similar to the CO(3-2) fitting in Flaherty et al. (2015), we model CO(2-1) allowing q, T_mid0, T_atm0, R_c, inclination, and v_turb to vary, along with the addition of R_in, R_break, and q_in. Results are listed in Table 4 and PDFs are shown in Figure 14. The inner excess extends from ${11}_{-3}^{+5}\,\mathrm{au}$ out to ${70}_{-5}^{+17}\,\mathrm{au}$ with a temperature power-law slope of $-{0.57}_{-0.04}^{+0.05}$ in the inner disk and $-{0.27}_{-0.02}^{+0.01}$ in the outer disk. These slopes are similar to that predicted by radiative transfer models (D'Alessio et al. 2006; Boneberg et al. 2016). We are also able to limit turbulence to <0.06c_s.

Zoom In Zoom Out Reset image size

Table 4.  CO(2-1) Model Result

Parameter	Value
q	$-{0.27}_{-0.02}^{+0.01}$
log(Rc (au))	${2.30}_{-0.02}^{+0.03}$
vturb	<0.06cs
Tatm0	${87}_{-2}^{+1}\,{\rm{K}}$
Tmid0	${17.8}_{-0.6}^{+0.8}\,{\rm{K}}$
incl	${47.5}_{-0.2}^{+0.5}$
Rin	${11}_{-3}^{+5}\,\mathrm{au}$
Rbreak	${70}_{-5}^{+17}\,\mathrm{au}$
qin	$-{0.57}_{-0.04}^{+0.05}$
${\chi }_{\nu }^{2}$ a	1.16

Note. Median, plus 3σ ranges, of the posterior distribution functions derived from fitting to the CO(2-1) data.

^aReduced chi-squared for the model defined by the median of the posterior distribution functions.

Download table as:  ASCIITypeset image

There is a modest difference in temperature structure between the CO(2-1) fit (${T}_{\mathrm{atm}0}={87}_{-2}^{+1}\,{\rm{K}}$, $q=-{0.27}_{-0.02}^{+0.01}$) and the assumed structure in the modeling of DCO⁺ and C¹⁸O(2-1) (T_atm0 = 93.8 K and q = −0.216). This difference in temperature structure will most significantly affect the upper reaches of the atmosphere (e.g., z > 40 au at a radius of 100 au), which is above the emitting regions of these molecules. The difference in temperature in the CO, C¹⁸O, and DCO⁺ emitting regions between the structure derived with CO(2-1) and that derived previously from CO(3-2) is <5% and will not significantly affect our results.

We find a good match to the data, and the best fit has removed much of the inner excess, leaving only residuals at ∼3% of the peak flux. In the channel maps, we find residuals in the central velocity channels (Figure 15). It appears as though our model produces emission that is more flat-topped along the azimuthal direction than the data. A similar effect can be seen in Figure 13 of Simon et al. (2015), where they employ a density and temperature structure similar to that used here, but fed through the LIME radiative transfer code (Brinch & Hogerheijde 2010). This suggests that these residuals are not artifacts of the ray-tracing code, but are reflective of the assumed model structure.

Zoom In Zoom Out Reset image size

These residuals lead to a large discrepancy between the peak-to-trough ratio of the model spectrum and that of the data. Simon et al. (2015) found that the peak-to-trough ratio of an emission spectrum varies with the turbulence, with smaller peak-to-trough ratios associated with stronger turbulence. The large peak-to-trough ratio in the model, as compared to the data, could potentially be interpreted as a sign that the model has underestimated the turbulence, although a closer inspection reveals that this is not the case. Figure 15 shows a model with substantial turbulence (v_turb = 0.4c_s) that better fits the shape of the spectra. The channel maps reveal that this spectral match is misleading. The high-turbulence model has not eliminated the positive residuals near line center, but has simply introduced additional negative residuals that cancel the positive residuals when integrated over the entire image. The negative residuals fall along the edge of the observed emission, indicating that the model overestimates the spatial broadening in the image plane, which is a sign of an overestimate of the turbulence. This indicates that while our model is not a perfect fit to the data, we are likely not underestimating the turbulence.

5.2.3. Origin of the Inner Excess

Similar inner excess features have been observed in the gas surrounding TW Hya (Rosenfeld et al. 2012) and the dust surrounding IM Lup (Cleeves et al. 2016). Here we are able to model much of the emission by varying the slope of the radial temperature gradient, as predicted by radiative transfer models of flared protoplanetary disks (e.g., D'Alessio et al. 2006).

The exact shape of this inner excess depends on the vertical region of the disk being traced, as well as any other parameters with a radial profile (e.g., molecular abundances, for optically thin lines) that can potentially affect the slope. Differences between the midplane temperature profile, measured by C¹⁸O, and the atmosphere temperature profile, measured by CO(2-1), can lead to differences in q_in. D'Alessio et al. (2006) find that while the midplane temperature goes as R^−0.75 in the inner disk, the surface temperature profile is closer to R^−0.5. C¹⁸O may also be subject to selective photodissociation, which can deplete its abundance in the inner disk relative to the outer disk and relative to CO (Miotello et al. 2014). This can explain the steeper value of q_in and larger inner hole, derived for C¹⁸O relative to CO. Other sources of an inner excess, such as a warp (Rosenfeld et al. 2012), or high-altitude dust grains leading to additional heating of the disk (Cleeves et al. 2016), could potentially explain the observations. More detailed modeling that fully accounts for radiative transfer effects between different tracers is needed to constrain the combined inner and outer disk structures.

In Flaherty et al. (2015), we used CO(3-2) to constrain the non-thermal line width in the upper layers of the disk to < 0.04c_s. At the same time, we used Science Verification data of C¹⁸O(2-1) emission to constrain the motion closer to the midplane, but given the modest S/N of the data, we were only able to place weak constraints on the turbulence (<0.4c_s). Here, with higher S/N C¹⁸O(2-1) data and new DCO⁺(3-2) data, we can much more accurately characterize the turbulent motion toward the midplane of the disk. With DCO⁺(3-2) and C¹⁸O(2-1), we place upper limits of <0.04c_s and <0.05c_s, respectively, while the CO(2-1) upper limit of <0.06c_s confirms the weak turbulence in the upper layers found with CO(3-2). At the midplane in the outer disk, this corresponds to ≲10 m s⁻¹. Under the standard α-prescription for turbulent viscosity (ν = αc_sH Shakura & Sunyaev 1973), in which $\alpha \sim {({v}_{\mathrm{turb}}/{c}_{s})}^{2}$, our velocity limits imply α < 1.5 × 10⁻³, <2.4 × 10⁻³, and <3.2 × 10⁻³ for DCO+(3-2), C¹⁸O(2-1), and CO(2-1), respectively. Boneberg et al. (2016), in their modeling of the spectral energy distribution and the inner 90 au of the Science Verification C¹⁸O(2-1) data, are able to reproduce the data with α of (0.1–6.3) × 10⁻³, consistent with our observations. For DCO⁺ we have verified that our limits are robust against uncertainties in the amplitude calibration, the underlying temperature structure, and assumptions about the distribution of DCO⁺ throughout the disk, with these effects only leading to an increase of the upper limit to <0.12c_s.

A generic feature of MRI models of turbulence is the vertical gradient in velocity (Fromang & Nelson 2006; Simon et al. 2013, 2015). In these simulations, the turbulence reaches 0.3c_s–1c_s at >3H, where H is the pressure scale height, while even in full ideal MHD simulations (e.g., Fromang & Nelson 2006) it drops down to 0.05c_s–0.1c_s at the midplane. When accounting for ambipolar diffusion, an important non-ideal MHD effect in the outer disk regions probed by our observations, the midplane turbulence can drop as low as ∼0.01c_s (Simon et al. 2013, 2015; Bai 2015). The disk surface, within a column density of Σ = 0.01–0.1 g cm⁻², is still expected to be fully turbulent due to ionization by far-ultraviolet (FUV) photons (Perez-Becker & Chiang 2011), leading to motions similar to those in the ideal MHD simulations, if FUV radiation is able to penetrate into the outer disk.

To understand how our observations relate to these theoretical predictions, we first must determine the vertical origin of the CO, C¹⁸O, and DCO⁺ emission. To do this, along each line of sight through the disk we calculate the heights above which 5% and 95% of the emission arises in the context of our best-fit model, defining a band that encompasses 90% of the line emission. We then generated an intensity-weighted average emission band as a function of radius, shown in Figure 16. The optically thick CO(3-2) and CO(2-1) lines sample the upper edge of the molecular layer with emission originating from ∼5H within 100 au and ∼2H in the outer disk, where H is defined based on the midplane temperature. Our limit on the non-thermal motion is likely more reflective of the z ∼ 2H region beyond 100 au, rather than the z ∼ 5H region within 100 au because of the larger surface area of the outer disk. At these locations in the disk, our stringent limits stand in contrast to theoretical predictions. Changes in X-ray ionization or magnetic field strength may lead to these weak non-thermal motions (Simon et al. 2013; J. B. Simon et al. 2017, in preparation).

Zoom In Zoom Out Reset image size

C¹⁸O(2-1) has a much lower optical depth than CO(2-1), allowing emission from closer to the midplane to escape the disk, although it is still subject to freeze-out in the outer disk. Within 100 au, C¹⁸O(2-1) originates from below 3H while at larger radii it traces the edges of the CO condensation front. DCO⁺, due to its low optical depth, traces the cold gas near the midplane throughout much of the disk. Our limits at these small heights are more similar to the theoretical predictions, especially when accounting for damping by ambipolar diffusion and the uncertainty in temperature on our DCO⁺ measurement.

If the outer disk is massive enough and cool enough, it may be subject to gravitoturbulence, which predicts motions of 0.2c_s–0.4c_s with little variation with height (Shi & Chiang 2014), and large radial variations only among the most massive systems that are capable of driving spiral arms (Forgan et al. 2012). The lack of strong non-thermal motion in HD 163296's disk suggest that this is not happening at a detectable level. Complex motions can also be driven by Jupiter-mass planets (e.g., Fung & Chiang 2016). The recent discovery of multiple dark rings in the dust continuum (Isella et al. 2016) hints at the presence at Saturn-mass planets. The lack of strong non-thermal, non-Keplerian motion within the gas disk is consistent with these planets being sub-Jupiter mass.

Regardless of its nature, the lack of strong turbulence has important implications for chemical mixing and dust settling. Turbulent eddies can lead to mixing of chemical species in the vertical direction (Semenov & Wiebe 2011; Furuya & Aikawa 2014). Sharp chemical boundaries, such as condensation fronts, can be softened with mixing. Xu et al. (2017) also find that CO gas diffuses downward into the freeze-out zone, and without strong turbulent mixing the CO ice-coated dust particles are not brought up into the warmer atmosphere where they can sublimate, eventually leading to a depletion of CO from the gas phase in low-turbulence systems. Furuya & Aikawa (2014) find a similar result, depleted CO abundance in a low-turbulence disk, due to the build up of complex ices as a sink of carbon and oxygen in the midplane. When modeling C¹⁸O(2-1), we find that X_CO is diminished from the standard ISM value by a factor of 4.06 ± 0.04, consistent with these predictions. We do caution that given the degeneracies in our modeling approach, the decrease in X_CO may actually reflect a decrease in ¹⁸O abundance, dust-to-gas mass ratio, or a combination of all three.

Dust grains are also affected by turbulence through their interactions with the gas. In the vertical direction, dust grain motion is subject to the competing effects of settling toward the midplane and turbulence lifting them back up into the disk atmosphere (Dullemond & Dominik 2004). This lifting can occur even for modest levels of turbulence, v_turb ∼ 0.1c_s (Fromang & Papaloizou 2006), and is especially prominent for the smallest grains (Fromang & Nelson 2009). The dust settles to a scale height of ${h}_{p}/H\sim \sqrt{(}\alpha /\tau )$ when α < τ, where τ is the optical depth (Youdin & Lithwick 2007). With our disk model, millimeter-sized grains in the outer disk roughly correspond to τ ∼ 0.1, which, when combined with α < 1.5 × 10⁻³ from DCO⁺, implies h_p/H<0.12, or h_p/R < 0.007. Dust settling has been used as circumstantial evidence for α ∼ 10⁻⁴ (Mulders & Dominik 2012; Pinte et al. 2016) and settling has been inferred in the disk around HD 163296 from its diminished infrared excess (Meeus et al. 2001; Juhász et al. 2010), although these infrared measures are more sensitive to the inner disk, while our ALMA constraints apply to the outer (≳60 au) disk. In modeling the C¹⁸O and dust emission, Boneberg et al. (2016) require a level of dust settling associated with α ∼ (0.1–6.3) × 10⁻³, consistent with our directly measured constraints. Enhanced settling associated with low turbulence may also lead to less grain growth across condensation fronts (Ros & Johansen 2013). In general, smaller α leads to weaker collisional velocities (Ormel & Cuzzi 2007) and subsequently to enhanced grain growth, as has been inferred in the disk around HD 163296 based on the submillimeter spectral index (Natta et al. 2004; Isella et al. 2007; Guidi et al. 2016). Given the lack of strong turbulence in the disk around HD 163296, radial drift may dominate over the vertical motion in setting the trajectories of the grains (Turner et al. 2006; Zhu et al. 2015). The observational evidence for radial drift in the disk around HD 163296 (Guidi et al. 2016) supports the importance of this process in this system.

We have presented ALMA observations of DCO⁺(3-2), CO(2-1), and C¹⁸O(2-1) that are able to constrain the non-thermal line width in the disk around the Herbig star HD 163296 (<0.04c_s, <0.06c_s, and <0.05c_s respectively). For DCO⁺, we found that the constraint on turbulence is robust against uncertainties in the radial width of the rings and the uncertainty in the amplitude calibration, with a modest degeneracy due to uncertainty in the midplane temperature. In general, our results are relatively robust to our choice of temperature and density functional forms because each emission line traces a relatively narrow vertical region of the disk, and many different model structures can be made to pass through the constraints on temperature and density in this region. The complementary emitting regions of the three lines allow us to map the vertical structure of the turbulence, and we can rule out strong turbulence from very close to the midplane up to the surface layers between radii of ∼30 au (the inner edge of the C¹⁸O(2-1) emission) and ∼300 au (the outer edge of detectable DCO⁺ emission). The limits on the turbulence in the upper layers of the disk in particular fall below the motions predicted by typical models of both MRI and gravitoturbulence.

In modeling the disk, we also find that the DCO⁺ emission is confined to three distinct rings. These rings are likely a reflection of a complex chemical structure, possibly due to non-thermal desorption of CO or the inward migration of dust grains. Bright rings of dust (Isella et al. 2016; Zhang et al. 2016) are not collocated with the DCO⁺ rings, consistent with the DCO⁺ structure arising due to complex chemical processes rather than variations in the underlying gas surface density. In the CO and C¹⁸O emission, we find an inner excess that can be explained by a change in the temperature structure toward the inner disk. These observations highlight the complex structure of planet-forming disks around young stars.

We thank the referee for thoughtful comments that greatly improved the manuscript. K.F. and A.M.H. are supported by NSF grant AST-1412647. S.R. was supported by the Keck Northeast Astronomy Consortium REU program through NSF grant AST-1005024. J.B.S. and P.J.A. are supported by NSF grant AST-1313021. We thank I. Cleeves, J. Huang, R. Teague, and A. Isella for valuable discussions. This paper makes use of the following ALMA data: ADS/JAO.ALMA#2013.1.00366.S and ADS/JAO.ALMA#2011.0.00010.SV. ALMA is a partnership of ESO (representing its member states), NSF (USA), and NINS (Japan), together with NRC (Canada), NSC and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The 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. The research made use of Astropy, a community-developed core Python package for Astronomy (Astropy Collaboration et al. 2013). We thank Wesleyan University for time on its high performance computing cluster supported by the NSF under grant number CNS-0619508. This work was supported by the Momentum grant of the MTA CSFK Lendület Disk Research Group.

Figure 17 shows the progress of the chains for our DCO⁺(3-2) fiducial model fit, described in Section 5.1. The walkers quickly converge and settle around the best fit by step 500, with most parameters reaching the best fit by step 100. The average acceptance fraction for the walkers is 0.51.

Zoom In Zoom Out Reset image size