WEAK TURBULENCE IN THE HD 163296 PROTOPLANETARY DISK REVEALED BY ALMA CO OBSERVATIONS

In hydrostatic equilibrium, the gas density and temperature are related by:

$$\begin{eqnarray}&&-\displaystyle \frac{\partial \mathrm{ln}{\rho }_{{\rm{gas}}}}{\partial z}=\displaystyle \frac{\partial \mathrm{ln}{T}_{{\rm{gas}}}}{\partial z}+\displaystyle \frac{1}{{c}_{{\rm{s}}}^{2}}\left[\displaystyle \frac{{{GM}}_{*}z}{{({r}^{2}+{z}^{2})}^{3/2}}\right],\end{eqnarray} \tag{ 1 }$$

where ρ_gas is the gas density, T_gas is the gas temperature, M_* is the stellar mass, and r, z are the radial and vertical position within the disk. The sound speed is given by ${c}_{{\rm{s}}}^{2}={k}_{{\rm{B}}}{T}_{{\rm{gas}}}/\mu {m}_{{\rm{h}}}$, where k_B is Boltzmann's constant and m_h is the mass of hydrogen. The disk is assumed to be made of 80% molecular hydrogen with μ = 2.37.

At each radius the vertical density structure is calculated from the temperature profile. For the temperature structure we use functional form laid out in Rosenfeld et al. (2013a),

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

which is based on the Dartois et al. (2003) Type II structure and has been used to successfully model the HD 163296 disk CO emission previously. This temperature profile increases from T_mid at the midplane to T_atm at a height of Z_q. The midplane remains cold because it is shielded from the stellar flux that heats the surface layers. Mechanical heating from the accretion flow can become important in the midplane of the inner disk (D'Alessio et al. 2006) but we ignore this since our measurements are most sensitive to the outer disk where stellar irradiation dominates. The values of T_mid, T_atm, and Z_q vary with radius according to

$$\begin{eqnarray}{T}_{{\rm{mid}}} & = & {T}_{{\rm{mid0}}}{\left(\displaystyle \frac{r}{150\;\mathrm{AU}}\right)}^{q}\\ {T}_{{\rm{atm}}} & = & {T}_{{\rm{atm0}}}{\left(\displaystyle \frac{r}{150\;\mathrm{AU}}\right)}^{q}\\ {Z}_{q} & = & {Z}_{q0}{\left(\displaystyle \frac{r}{150\;\mathrm{AU}}\right)}^{1.3}.\end{eqnarray} \tag{ 3 }$$

This type of profile for the vertical and radial temperature structure is consistent with detailed models of the gas (Jonkheid et al. 2007; Walsh et al. 2010) and the dust (D'Alessio et al. 2006). We only consider the gas temperature in our models and are not concerned with any possible differences between the gas and dust temperature that can arise in the uppermost layers of the disk (Jonkheid et al. 2007).

At each radius the vertically integrated mass is set by the surface density profile, which we take to follow the form:

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

The value of Σ_c depends on the total gas mass, the radial power law index of the surface density, and the critical radius (${{\rm{\Sigma }}}_{c}={M}_{{\rm{gas}}}(2-\gamma )/(2\pi {R}_{{\rm{c}}}^{2})$). Unless otherwise specified, we assume M_gas = 0.09 M_⊙ (Isella et al. 2007). The critical radius R_c defines the turnover between a power law and a sharp exponential decay. This is motivated theoretically by a similarity solution to the disk structure calculation for a disk with viscosity that scales with radius (ν ∝ r^γ) (Lynden-Bell & Pringle 1974; Hartmann et al. 1998). Andrews et al. (2010) use this surface density profile to model resolved sub-millimeter dust emission from a collection of young solar-type stars and find γ typically close to 1 with R_c ranging from 15 to 200 AU. This functional form has also been successful in fitting both the dust and the gas (Hughes et al. 2008) including the different apparent radii for the gas and dust disk, although recent evidence suggests that in some systems the dust disks are smaller than the gaseous component (Andrews et al. 2012; de Gregorio-Monsalvo et al. 2013; Walsh et al. 2014).

In a disk with an isothermal vertical profile, the pressure scale height H is simply defined as the width of the Gaussian vertical density profile. With our smoothly varying temperature profile, the definition of H is not as straightforward, although it is still useful to estimate H for comparison with literature models. When used we define H as c_s/Ω, where the midplane temperature sets the sound speed.

The velocity field is taken to be Keplerian, with modifications due to the height above the midplane and the pressure support of the gas:

$$\begin{eqnarray}&&\displaystyle \frac{{v}^{2}}{r}=\displaystyle \frac{{{GM}}_{*}r}{{({r}^{2}+{z}^{2})}^{3/2}}+\displaystyle \frac{1}{{\rho }_{{\rm{gas}}}}\displaystyle \frac{\partial {P}_{{\rm{gas}}}}{\partial r},\end{eqnarray} \tag{ 5 }$$

where ρ_gas and P_gas are the gas mass density and pressure, respectively. Rosenfeld et al. (2013a) estimate that these modifications change the velocities by a few percent in the upper layers of the disk. We do not include the self-gravity of the disk since this has a negligible effect on the velocity field, especially for our small gas mass (Rosenfeld et al. 2013a).

Gas-phase CO is present in the regions of the disk where it is not frozen out onto dust grains, or photo-dissociated by high energy stellar photons. In modeling SMA observations of CO emission from the HD 163296 system, Qi et al. (2011) find a freeze out temperature of 19 K, which we take as the freeze-out boundary in our models. We simulate freeze-out by decreasing the CO abundance by eight orders of magnitude in regions where the gas temperature drops below 19 K.

Photo-dissociation occurs above a height where the vertically integrated column density is

$$\begin{eqnarray}&&0.706{\displaystyle \int }_{{z}_{\mathrm{phot}}}^{\infty }{n}_{{\rm{gas}}}(r,z){dz}\lt {\sigma }_{s}.\end{eqnarray} \tag{ 6 }$$

This is the same boundary used in Qi et al. (2011) with σ_s = 0.79 × 1.59 × 10²¹ cm⁻². This boundary is defined by the absorption of UV photons from the central source (Aikawa & Herbst 1999) and is included in our model as a decrease in the abundance of CO by eight orders of magnitude. The exact location of the photo-dissociation boundary may vary between different isotopologues due to selective self-shielding (Visser et al. 2009; Miotello et al. 2014). When modeling ¹³CO and C¹⁸O emission we initially assume ISM values of ¹²C/¹³C = 69 and ¹⁶O/¹⁸O = 557 (Wilson 1999), although we also consider models in which these isotope abundances, as well as the CO/H₂ abundance, are allowed to vary. Unless otherwise specified we assume CO/H₂ = 10⁻⁴. Our simple chemical model for CO ignores much of the complexity of disk chemistry (Jonkheid et al. 2007; Walsh et al. 2010) but has proven successful in fitting the CO emission from the HD 163296 system in the past (Qi et al. 2011; Rosenfeld et al. 2013a).

Once the structure of the disk is set, we calculate the flux through the disk. The line intensity is

$$\begin{eqnarray}&&{I}_{\nu }={\displaystyle \int }_{0}^{\infty }{S}_{\nu }(s)\mathrm{exp}[-{\tau }_{\nu }(s)]{K}_{\nu }(s){ds},\end{eqnarray} \tag{ 7 }$$

where s is the linear coordinate along the line of sight increasing outward from the observer. The optical depth is ${\tau }_{\nu }(s)={\displaystyle \int }_{0}^{s}{K}_{\nu }({s}^{\prime }){{ds}}^{\prime }$ where K_ν(s) is the absorption coefficient and S_ν(s) is the source function. For simplicity we assume the system is in local thermodynamic equilibrium (LTE) with the level populations given by the Boltzmann equation and the local gas temperature and the source function is approximated by the Planck function. The CO emitting region is at a density that is much higher than the critical density (∼10⁴ cm⁻³), allowing us to assume LTE without much loss of accuracy (Pavlyuchenkov et al. 2007)

The line is assumed to be a Gaussian with width

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

and a FWHM of $2\sqrt{\mathrm{log}2}$ΔV. Throughout much of our modeling, we scale v_turb to the local sound speed (e.g., ${\rm{\Delta }}V=\sqrt{(2{k}_{{\rm{B}}}T(r,z)/{m}_{\mathrm{CO}})+{({v}_{{\rm{turb}}}/{c}_{{\rm{s}}})}^{2}}$). The sound speed, which is a factor of $\sqrt{{m}_{\mathrm{CO}}/(\mu {m}_{{\rm{H}}})}\approx 3.4$ larger than the thermal broadening of CO, is expected to set the velocity scale for turbulence (Balbus & Hawley 1998) and this parameterization allows turbulence to slowly vary with radius and height within the disk. Since we cannot resolve the spatial scale of the turbulence, we do not detect the bulk motion, v_turb, associated with the eddies but instead measure the rms of the velocity distribution, δv_turb. For a Maxwell–Boltzmann distribution of velocities, which is expected from numerical simulations of turbulence (Simon et al. 2015), the rms and mean velocity are the same to within factors close to unity. Given the similarity of these factors, and the prevalence of the notation v_turb in the literature, we use v_turb to denote turbulence even though strictly speaking we are measuring the rms width, δv_turb, of the velocity distribution.

Finally, we keep the distance (d = 122 pc), stellar mass (M_* = 2.3 M_⊙) and position angle (PA = 312^o) fixed. The offset of the disk from the phase center and velocity center of the spectrum do not strongly depend on the other model parameters, and were adjusted to minimize residuals using initial model fits.

Using the model described above, we can generate synthetic images for a given set of parameters. We use the MIRIAD task UVMODEL to generate model visibilities sampled at the same spatial frequencies as the data. Even though the data contain both XX and YY polarizations, we only use the total intensity (I = (XX + YY)/2) to compare with the model because the line emission is not expected to be polarized. The goodness of fit between the model (V_mod) and observed (V_obs) visibilities is calculated using the chi-squared statistic where the uncertainty at each uv point is derived from the dispersion in the data, as described earlier.

To sample the posterior probability distribution for each of the parameters we use a MCMC technique. In particular, we employ the affine-invariant routine EMCEE (Foreman-Mackey et al. 2013) based on the algorithm originally presented in Goodman & Weare (2010). As opposed to a Metropolis–Hastings chain, the affine-invariant method uses a large number of walkers whose movements through parameter space are proposed along lines to the position of the other walkers. This has the advantage of requiring less fine-tuning to accurately sample the posterior distribution function (PDF) and can efficiently probe the PDF even when there are strong degeneracies between parameters, as is sometimes the case with our models. After an initial burn-in period the positions of the walkers along the chains represent independent samples from the PDF and the distribution of the walkers in parameters space can be used as an estimate of the PDF. Our chains typically consist of 40 walkers and 3000 steps. We only examine the last 1000–2000 steps in each chain, corresponding to >10 autocorrelation times, resulting in 40,000–80,000 samples of the PDF. The initial positions of the walkers are chosen to be evenly distributed throughout parameter space. We require that M_gas, R_c, and v_turb are positive with uniform priors across the parameter space.

Systematic uncertainties in the amplitude calibration of the data can have a large effect on the results. ALMA line fluxes differ from those previously measured by the SMA for these same lines (Hughes et al. 2011; Qi et al. 2011) by ∼5%, which is within the expected 20% limit on the systematic uncertainty of the ALMA and SMA data. Since CO(3-2) is optically thick its flux is directly proportional to the temperature of the τ = 1 surface and a systematic uncertainty in the total line flux will directly translate into an uncertainty in the temperature. Since both temperature and turbulence contribute to broadening of the line, an uncertainty in one will translate into an uncertainty in the other (although, as discussed below and in Simon et al. (2015), the spatial resolution of the images and the peak-to-trough ratio of the line profile help to break some of this degeneracy). For the optically thin lines whose flux is proportional to both the surface density and the temperature, the exact influence of the systematic uncertainty is less obvious. This is especially important when combining information from multiple lines, which are subject to different systematic uncertainties.

To account for the systematic uncertainty, we scale all of the model visibilities by a factor of s. This factor is defined such that a value larger than one corresponds to the true disk flux being brighter than the data (e.g., s = 1.2 implies the true flux is 20% brighter than the measured flux), while a value less than one corresponds to a true disk flux weaker than the measured flux. For much of our modeling this is kept fixed at s = 1, but for CO(3-2) we consider fitting when s is fixed at 0.8 or 1.2, and we treat it as a free parameter when fitting the four lines simultaneously. Despite the large number of degrees of freedom, ∼millions for the single line fits, including too many free parameters prevents the walkers from converging on a best fit model, especially when there are strong degeneracies involved. For this reason we fix s in the single line fits. In the multi-line fit there is enough additional information to allow for a constraint on s, and there we consider situations where s is left as a free parameter. The total number of free parameters depends on the line being fit and ranges from four (v_turb, γ, inclination, N(X)/N(¹²C¹⁶O)) when fitting ¹³CO and C¹⁸O to 10 (q, R_c, v_turb, T_atm0, γ, T_mid0, inclination, log(CO/H₂), s₃₂, and s₂₁) when fitting all four emissions lines simultaneously while accounting for the systematic uncertainty.

We start by modeling the high-resolution CO(3-2) line using 40 walkers over 1500 steps. The first 1000 steps are excluded to remove the burn-in phase of the MCMC chains; with typical auto-correlation times of 20–80 steps this corresponds to when the walkers have settled around the best fit. Since CO(3-2) is optically thick it is not sensitive to the surface density, and we fix M_gas = 0.09 M_⊙ and γ = 1, consistent with previous models (Isella et al. 2007; Qi et al. 2011). We allow R_c to vary to control the radial extent of the emission. By spatially resolving the front (z > 0) and back (z < 0) side of the disk, along with the cold midplane, we effectively have two data points to confine T(z), which is not enough to fully constrain the three parameters that control T(z) (T_mid0, T_atm0, and Z_q0). Our initial attempts at MCMC modeling found that T_atm0 and Z_q0 were highly degenerate, which slowed significantly the convergence of the walkers. Given the preference of the models for higher Z_q0 we fix Z_q0 = 70 AU. With many resolution elements across the radius of the disk, we can constrain the radial profile of the temperature, and allow q to vary. Finally, we include turbulence as well as inclination.

Results are listed in Table 1, along with the three-sigma uncertainties about the median. Channel maps for the best-fit model are shown in Figure 1. We also show residual channel maps (subtracted in the visibility domain, ie. ${\rm{\Delta }}{I}_{V}=({V}_{{\rm{obs}}}(u,v)-{V}_{{\rm{mod}}}(u,v))$) in Figure 2. This model reproduces much of the morphology of the disk, including the front/back emission with a drop in emission at the cold midplane, and has a reduced chi-squared close to one. There are still some residuals between the model and the data but these may be due to deviations from our simple functional form for the temperature and density structure, or even non-circular motion, that will not be perfectly captured in our model. A fit directly to the low-resolution data returns a nearly identical, although slightly weaker limit on turbulence (Table 1, Figure 3), while also producing a good fit to the data. Slight variations in the temperature structure are consistent with degeneracies and are well below the systematic uncertainties, discussed below. The similarity in the fit to these spectra with different velocity resolutions indicates that the spatial information is providing much of the constraint on the models. Given the consistency between the high and low-resolution data results, while we report the structure from the high-resolution PDFs, we explore various aspects of the fitting using the low-resolution spectra.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 1.  Model Fitting Results

Line	q	Rc(AU)	vturb/cs	Tatm0(K)	γ	Tmid0(K)	Inclination	N(X)/N(12C16O)	${{\boldsymbol{\chi }}}_{\nu }^{2}$ a
CO(3-2) (high-res)	$-{0.278}_{-0.008}^{+0.003}$	${175}_{-4}^{+13}$	<0.031	86 ± 1	1b	${21.8}_{-0.4}^{+0.7}$	47.6 ± 0.1	1	0.96
CO(3-2) (low-res)	$-{0.216}_{-0.008}^{+0.01}$	${194}_{-15}^{+12}$	<0.038	94 ± 2	1b	${17.5}_{-0.7}^{+0.8}$	48.3 ± 0.3	1	0.94
CO(2-1)	$-{0.27}_{-0.01}^{+0.02}$	${224}_{-16}^{+21}$	<0.31	${79.0}_{-1.5}^{+1.1}$	1b	17.5b	${46.1}_{-0.3}^{+0.4}$	1	0.92
13CO(2-1)	−0.216b	194b	<0.34	94b	${0.78}_{-0.01}^{+0.02}$	17.5b	47.5 ± 0.2	${233}_{-15}^{+3}$	0.94
C18O(2-1)	−0.216b	194b	<0.40	94b	1.13 ± 0.08	17.5b	${51.5}_{-0.3}^{+0.2}$	${1240}_{-110}^{+120}$	0.92

Notes. Median plus 3σ ranges for the PDFs defined by the walker positions after removing burn-in.

^aReduced chi-squared for the model defined by the median of the PDFs. ^bHeld fixed during the model fitting.

Download table as:  ASCIITypeset image

The temperature and density structure for the best-fit model are shown in Figure 4. The CO is confined to a thin layer that, at 200 AU, stretches from 10 to 60 AU above the midplane with temperatures ranging from 19 to ∼70 K. Previous models (Qi et al. 2011; Tilling et al. 2012; de Gregorio-Monsalvo et al. 2013) of the disk around HD 163296 employed radiative transfer calculations of the equilibrium disk temperature, assuming a dust distribution and a luminosity for the central source, and found similar temperatures to those in our models. This is not surprising since CO(3-2) is optically thick and strongly sensitive to the temperature. Our radial temperature profile falls off with a power-law index of $-{0.278}_{-0.008}^{+0.003}$ ($-{0.216}_{-0.008}^{+0.01}$ for the low-resolution fit), significantly shallower than the canonical value of −0.5, but consistent with values expected from radiative transfer models that take into account the flaring of the outer disk (D'Alessio et al. 2006). Flaring leads to more direct interception of starlight in the outer disk, thereby slowing the rate at which the temperature falls off with distance from the central star.

Zoom In Zoom Out Reset image size

While the statistical errors on the model parameters are small, they do not account for the systematic effects. For example, the uncertainty in distance strongly influences q since a more distant disk requires more extended radial emission which can be accomplished with a flatter radial temperature profile. Adjusting the distance by ± 10 pc, its 1σ uncertainty (Montesinos et al. 2009), and rerunning the fit to the low-resolution data we find that q varies by up to 20%; q rises from −0.236 ± 0.004 to $-{0.175}_{-0.02}^{+0.01}$ as distance increases, compared to $-{0.216}_{-0.008}^{+0.01}$ at the fiducial distance. The critical radius also increases with distance, from ${158}_{-5}^{+6}$ to ${203.7}_{-5}^{+19}$ AU for the near and far distances respectively, while the other model parameters are consistent within the statistical uncertainty. The differences between q and R_c when fitting the high-resolution and low-resolution data (Table 1) also highlight the degeneracy between these two parameters. The maximum spatial extent of the CO emission can be matched with a narrow disk (low R_c) and shallow temperature profile (high q) or a wide disk (high R_c) and a steeper temperature profile (low q). This effect is small, leading to 10%–20% differences in q and R_c, but is large compared to the statistical uncertainties.

Similarly, the statistical uncertainty on the atmosphere temperature is <5% although this does not include the uncertainty on the amplitude calibration. As discussed earlier, we account for the systematic uncertainty by applying a scale factor to the model and we run additional MCMC trials fitting the low-resolution data with s fixed at 1.2 and 0.8, mimicking a true disk flux that is 20% higher/lower than the data. The biggest effect is in T_atm0 which varies from ${84.3}_{-2.0}^{+0.7}$ to 101 ± 3 for low/high disk flux, respectively, compared to 94 ± 2 in the fiducial low-resolution model fit. The influence of systematic uncertainty dwarfs the statistical uncertainties for T_atm0. The variations of the other parameters in response to the uncertainty in the flux calibration are within the statistical errors derived from the PDF.

The atmosphere temperature is also subject to its degeneracy with the midplane temperature (Figure 5). While the emission from the bright upper half of the disk constrains the atmosphere temperature, the midplane temperature is constrained through the emission from the faint lower half of the disk and the lack of emission between these two features. In our fiducial model we find T_mid0 = ${21.8}_{-0.4}^{+0.7}.$ The constraints on q and T_mid0 in turn lead to a CO ice line that falls within the range R_ice = ${245}_{-13}^{+27}$ AU while previous observations of chemical tracers of CO freeze-out find that it is closer to 150–160 AU (Mathews et al. 2013; Qi et al. 2013b). Overestimating the midplane temperature results in a more puffed up disk, pushing the τ = 1 surface higher in the disk and the atmosphere temperature must be reduced accordingly to maintain the same temperature at the τ = 1 surface. This can be seen when looking at the results from fitting the low-resolution spectra which has different values of T_mid0 and T_atm0, and predicts a CO ice line at R_ice = ${103}_{-18}^{+23}$ AU, but produces similar images.

Zoom In Zoom Out Reset image size

While our models do not treat the dust, assumptions about the dust distribution may still be implicit in our models. Jonkheid et al. (2007) find that in their chemical models the depletion of dust from the upper layers of the dust leads to a drop in the gas temperature. In the context of the D'Alessio et al. (2006) models, Qi et al. (2011) show that the settling of large dust grains will affect the height at which the disk reaches the atmosphere temperature; in other words more settling results in lower Z_q0. Our assumed value of Z_q0 is most similar to models with little dust settling (Qi et al. 2011). The significant degeneracy between Z_q0 and T_atm0 means that a model with lower Z_q0, when combined with a lower T_atm0, would fit the data as well as the models presented here. Rosenfeld et al. (2013a), using the same functional form for the vertical temperature profile as we do here, find that the data can be fit with T_atm0 = 64 K, given Z_q0 = 43 AU. This atmosphere temperature is much lower than our T_atm0 = 86 ± 1 K but the difference in Z_q0 means that both models produce a nearly identical gas temperature at the τ = 1 surface. Further data probing the uppermost layers in the disk atmosphere, where CO is photodissociated, are needed to break this degeneracy.

Overall we are able to find a model that accurately fits the data and is consistent with previous radiative transfer models of the HD 163296 system. By carefully accounting for parameter degeneracies and systematic effects due to assumptions about the distance and flux calibration, we can measure the statistical and systematic uncertainties in temperature and density structure of the disk. This indicates that our derived disk structure is not going to strongly bias our characterization of turbulence.

4.1.1. Turbulence

When fitting either the high-resolution or low-resolution CO(3-2) spectra, we consistently find low levels of turbulence (Table 1). We conservatively quote three-sigma upper limits of v_turb < 0.03 c_s and v_turb < 0.04 c_s, respectively. In the uppermost layers of the CO region, these limits correspond to velocity dispersions of less than ∼9–19 m s⁻¹, with higher velocities at smaller radii. While we have some constraint on the inner disk from the high velocity channels, most of the information about the turbulence comes from the spatially resolved outer disk (R > 30 AU) and our upper limit is indicative of the behavior in this region. Our limit is well below the spectral resolution of even the high-resolution data (δv = 0.1 km s⁻¹). To test whether we can distinguish between turbulence at the resolution limit from much weaker non-thermal motion, we run an additional MCMC fit to the low-resolution CO(3-2) data with v_turb fixed at 0.1 km s⁻¹ while allowing the other model parameters to vary. Figure 6 shows the best-fit v_turb = 0.1 km s⁻¹ model along with our fiducial low turbulence model, using a Gaussian fit to the short-baseline visbilities to derive the spectra. The low turbulence models are a much better fit to the spectra, and the chi-squared, which captures the behavior of the full three-dimensional data set, also finds a significantly better fit (>10σ significance) from the low-turbulence models. This behavior indicates that fitting a model to the full data set has a stronger diagnostic power than the broadening of the line in the spectral domain. Simon et al. (2015) have suggested the peak-to-trough ratio as a potential diagnostic of turbulence. They found in simulated observations of numerical MRI simulations, with typical v_turb ∼ 0.1–0.5 c_s at the CO(3-2) emitting surface, that the peak-to-trough ratio decreases as the turbulence increases. This behavior can be seen in Figure 6 where the v_turb = 0.1 km s⁻¹ models have a much lower peak-to-trough ratio than the low-turbulence models and the data, consistent with the conclusion that the turbulence is weaker than 0.1 km s⁻¹.

Zoom In Zoom Out Reset image size

Peak-to-trough ratio is a one-dimensional metric that provides a convenient way of diagnosing differences in the three-dimensional data set. Variations in the temperature can also affect the peak-to-trough ratio, however the uncertainty in temperature is dominated by the uncertainty in the flux calibration and within the range of temperatures given by this uncertainty the peak-to-trough ratio is substantially more sensitive to turbulence than it is to temperature (Simon et al. 2015). By simultaneously fitting the temperature structure along with the turbulence across the full three-dimensional data set, we account for the various influences of the different parameters on the images as well as the spectral shape. The similarity between the turbulence results when fitting the low-resolution and high-resolution data suggest that the shape of the full line profile is not the dominant factor but instead spatial resolution plays a large role in constraining the non-thermal motion. In the individual channel maps turbulence broadens the width of the disk emission (Guilloteau et al. 2012; Simon et al. 2015). We demonstrate this effect for the central velocity channel of the HD 163296 system in Figure 7 by varying the turbulence while fixing the other model parameters. For low levels of turbulence there is little variation in the images, but the disk broadens dramatically for larger values of turbulence. This occurs even below the spectral resolution of the data, allowing us to push to low levels of turbulence.

Zoom In Zoom Out Reset image size

The broadening of the images not only changes their shape, but also changes their total flux. This can occur even when the broadening is below the spatial and spectral resolution. In this way line flux could be used as a diagnostic of turbulence, but this leaves it highly degenerate with temperature, and subject to the uncertainty in the calibration of the line flux. The systematic uncertainty is, conservatively, 20% at the ALMA wavelengths and dominates over the channel-to-channel uncertainties. As discussed earlier, we can account for systematic uncertainty by applying a scale factor to the model. When fixing s = 1.2, i.e., when the true disk flux is 20% brighter than the observed data, we find a limit of v_turb < 0.05 c_s, while if the true disk emission were 20% fainter than our data the limit on turbulence is v_turb < 0.08 c_s. Most of the change in disk flux is absorbed by the atmosphere temperature, with only small variations in turbulence. An uncertainty in distance could also affect our result, since the 1σ uncertainty on distance of 10 pc (Montesinos et al. 2009) leads to a 16% uncertainty in the total flux. In additional MCMC trials with distance fixed at 112 and 132 pc, we find that the limit on turbulence increases to v_turb < 0.05 c_s, 0.06 c_s, respectively. Even accounting for the systematic uncertainty and the distance uncertainty, we find that non-thermal motion is limited to very weak levels.

The high-resolution spectrum exhibits an asymmetry in the line profile, with the blueshifted peak brighter than the redshifted peak (Figure 6). One possible explanation for the behavior is that with an eccentric disk, the disk material "piles-up" on the slower moving side of the orbit leading to a higher flux from one side of the disk versus the other (e.g., Regály et al. 2011). Our model includes the dynamics of Keplerian rotation and thermal motion and any additional motion will likely be accounted for by the turbulence parameter. For this reason we choose to report our turbulence results as an upper limit rather than a detection. The posterior distribution for turbulence (Figure 3) has a shape characteristic of a detection but given the possibility of unaccounted features being modeled as turbulence, as well as the weak influence of turbulence on the data at such low levels (Figure 7), we conservatively quote an upper limit. Also, as noted earlier, our implementation of turbulence within the models means that we are really constraining the velocity dispersion, or rms velocity, associated with turbulence rather than the mean velocity within turbulent eddies.

As with CO(3-2), our CO(2-1) MCMC trial uses 40 walkers over 3000 steps. With a lower spatial resolution than the CO(3-2) data, the walkers here take longer to converge and the number of steps has been increased accordingly. CO(2-1) is also optically thick, making it sensitive to the temperature structure, but lacks the spatial resolution to distinguish the near and far side of the disk. For this reason we fix T_mid0 = 17.5, based on the low-resolution CO(3-2) results, while allowing q, R_c, v_turb, T_atm0, and inclination to vary during the fitting. Results, after removing the first 1000 steps, are listed in Table 1 with posterior distributions shown in Figure 8. Channel maps and imaged residual visibilities are shown in Figures 9 and 10, respectively. In the channel maps, outside of minor residuals, we are able to match the radial extent of the disk, as well as the flux throughout the spectra.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Differences in the model parameter PDFs exist between the CO(2-1) and CO(3-2) models, but these are most likely the result of degeneracies between parameters. In the CO(2-1) data we find a strong anti-correlation between R_c and q and compared to the low-resolution CO(3-2) fit, the CO(2-1) line settles on a model with steeper q and larger R_c, consistent with this anti-correlation. We also find a smaller inclination and lower T_atm0. The resulting disk structure derived from the two lines still has very similar temperatures at the top of the CO emitting surface, to which these data are particularly sensitive.

The upper limit on turbulence from CO(2-1) (v_turb < 0.31 c_s) is substantially higher than derived from the CO(3-2) line, most likely due to the lower spatial resolution of these data compared to CO(3-2). This limit corresponds to velocity widths of 200–150 m s⁻¹ in the upper layers of the disk. To see if a satisfactory fit can be found with low turbulence, we run an additional trial with v_turb fixed at 0.04 c_s, ie. the upper limit determined by the low-resolution CO(3-2) data. The resulting spectrum for this model is shown in Figure 11, along with the high turbulence CO(2-1) fit. The low turbulence fit underestimates the overall line flux, but better matches the line shape, in particular the peak to trough ratio. The imaged residuals for the low-turbulence model show features only on the north–east side of the disk (Figure 12). This asymmetry may be related to the temperature structure of the disk, as probed by an optically thick tracer. As discussed in Dartois et al. (2003), the line of sight to an inclined disk probes different heights for the area of the disk closer to the observer, in this case the south–west side, than on the far (north–east) side of the disk. This is because a line of sight through the near side of the disk ends at at smaller radius, where the temperature and densities are higher, than where it entered the disk. Conversely, a line of sight through the far side of the disk has its τ = 1 surface in an area with a lower density and temperature than where it entered the disk. While our ray-tracing code will account for photons coming from different surfaces, our assumed functional form for the temperature may not capture the difference in temperature between these surfaces. This would show up as asymmetric residuals, such as those seen in CO(2-1) (Figure 12) and in CO(3-2) (Figure 2). The MCMC routine cannot adjust the temperature parameters to account for the deviations from the assumed functional form, and tries to remove the residuals instead by increasing the turbulence, which raises the overall flux of the disk. Since the CO(2-1) data has a lower spatial resolution than CO(3-2), it has more room to increase the broadening of the images before the model becomes significantly wider than the data. The end result is a high value of turbulence, and a worse fit to the peak-to-trough ratio, due to limitations of the assumed temperature structure rather than large non-thermal motions within the disk. For this reason we quote an upper limit on turbulence despite the shape of the posterior distribution (Figure 8).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The decreased abundance of ¹³C relative to ¹²C reduces the optical depth of ¹³CO relative to ¹²CO which makes the line emission more sensitive to the total mass of the disk, but introduces a degeneracy between column density and temperature. To avoid this degeneracy we fix the temperature structure, defined by q, T_atm0, and T_mid0, based on the low-resolution CO(3-2) observations. We then allow γ to vary, along with v_turb and inclination. We fix the disk mass and R_c, while leaving the depletion of ¹³CO as a free parameter. During our initial attempts at MCMC modeling, and in the multiline fits discussed below, we found the best possible fit to the data came when the disk mass was fixed and the abundance was allowed to vary. Specifically we treat the ¹³CO/CO depletion as a free parameter, although in the case of a single line fit this is equivalent to allowing CO/H₂ to vary. Model fits, after removing the burn-in, are listed in Table 1 with posterior distributions shown in Figure 13. These models are able to accurately reproduce much of the ¹³CO emission, although there are some discrepancies near the line peaks (Figure 14). This is likey due to the temperature structure; the residuals show the same north-south asymmetry as in CO(2-1) and CO(3-2). In the channel maps and imaged residuals (Figure 14) the differences between the model and the data are small, with the strongest residuals only 10% of the peak flux.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As with CO(2-1), we put a limit on the turbulent broadening (v_turb < 0.34 c_s) that is less strigent than CO(3-2). At the coldest layers close to the midplane, just above CO freeze-out, this corresponds to a non-thermal velocity width less than ∼130 m s⁻¹. We also find ¹³CO/CO = ${233}_{-15}^{+3},$ significantly different from the ISM value of 69 (Wilson 1999). This may be a sign of unaccounted for CO chemistry, such as selective photodissociation (Miotello et al. 2014), or a sign that the assumed disk mass, which is derived from the dust emission, is an overestimate. This result is discussed in more detail below in the context of the multiline fit.

C¹⁸O, being more heavily depleted than ¹³CO, is fully optically thin and is able to trace the entire vertical extent of the disk. As with ¹³CO, we break the degeneracy between density and temperature by using the derived values of q, T_mid0, and T_atm0 from CO(3-2), while allowing the depletion and γ to vary. The resulting PDFs are shown in Figure 15 with channel maps and residuals in Figure 16. We find a significant degeneracy between the depletion factor and γ which is not surprising since both control the surface density of the disk. This anti-correlation may contribute to the difference in γ between the fit to C¹⁸O and the fit to ¹³CO; C¹⁸O implies a steeper surface density profile, but this can be brought in line with the ¹³CO result with a depletion factor of ∼1700. We note that the C¹⁸O line is fit with a higher inclination than the other lines. This difference is small, but significant. Further work is needed to fully explore the nature of this discrepancy.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The C¹⁸O line, with the lowest signal to noise, places the weakest constraint on the turbulent broadening, v_turb < 0.4 c_s. At the coldest layers just above the CO freeze-out temperature, this implies random velocity dispersions of less than 150 m s⁻¹. To check if a low turbulence model can adequately fit the data we ran an additional MCMC trial with the turbulence fixed at the low-resolution CO(3-2) limit (v_turb = 0.04 c_s) and the results are shown in Figure 17. The difference between the high turbulence fit and the low turbulence fit is small despite the order of magnitude difference in turbulence. While C¹⁸O is more sensitive to the midplane, its optically thin nature makes it less sensitive to velocity broadening from turbulence. In an optically thick line, the velocity shear provided by turbulence allows more photons to escape, increasing the flux, while in optically thin lines this effect is much smaller (Horne & Marsh 1986). This is consistent with the behavior seen in Figure 17 for C¹⁸O, as well as the weak limit on turbulence in ¹³CO.

Zoom In Zoom Out Reset image size

The four lines, with their different optical depths, probe different regions within the disk. Figure 18 shows the height of the τ = 1 surface for the best fit models to each of the four lines. At large radii, CO(3-2) is sensitive to materal at z ∼ 100 AU, while ¹³CO and C¹⁸O probe material closer to the midplane. By fitting one model to all four lines we can leverage this complementary information to more accurately constrain the temperature, density, and turbulence throughout the disk. For this MCMC trial we employ 60 walkers over 1000 steps, with the first 500 removed as burn-in. We fix M_gas = 0.09 M_⊙ but allow q, R_c, v_turb, γ, T_atm0, T_mid0, and inclination to vary, along with the CO/H₂ abundance. We find that allowing the CO/H₂ to vary, instead of the individual ¹³CO and C¹⁸O depletion factors or the gas mass, provides the best multiline fit and we report those results here. Numerical simulations of MRI predict strong vertical gradients in the turbulent motion as a function of height within the disk (Fromang & Nelson 2006; Simon et al. 2013) but given the lack of detection found with the other lines we do not attempt to include any change in turbulence with height, beyond allowing it to vary with the local sound speed.

Zoom In Zoom Out Reset image size

The results are listed in Table 2 with visibility spectra shown in Figure 19, and maps of the central velocity channel for each line in Figure 20. In terms of the overall density and temperature structure, while some of the parameters change compared to where they stood with the individual line fits, the overall structure is similar. Figure 21 shows the density and temperature structure for the model defined by the median of the PDFs, which is similar to the high-resolution CO(3-2) line result. We find that CO is depleted by a factor of ${4.9}_{-0.4}^{+0.6}$ relative to the canonical ISM value. The individual fits to ¹³CO and C¹⁸O suggest a depletion of only a factor of ∼3, which is smaller than that derived here because of the lower midplane temperature assumed in those fits. Recent observations have found mixed results when measuring CO/H₂ (e.g., Favre et al. 2013; France et al. 2014), while models suggest that complex gas-grain chemistry in the cold midplane can affect the CO/H₂ ratio (Furuya & Aikawa 2014; Reboussin et al. 2015). It is also possible that the deviation in CO/H₂ is a sign that the assumed disk mass is incorrect. The mass was derived from the sub-millimeter emission of the dust (Isella et al. 2007) and could deviate from the true value if there is a change in the gas to dust mass ratio, a concentration of dust into very large grains that are invisible in the sub-millimeter, or if the dust does not follow the gas distribution. We have run additional trials fixing the CO/H₂ ratio while allowing the disk mass to vary and find M_gas = 0.08 ± 0.01 M_⊙, which is similar to our assumed value, but is subject to large degeneracies, especially with R_c and produces a worse overall fit to the data. Further work, including a more physically realistic treatment of CO chemistry near the ice line, is needed to confirm and understand the low CO/H₂ ratio.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 2.  Multi-line Model Fit

Parameter	Result
q	${-.298}_{-0.008}^{+0.017}$
Rc	${195}_{-17}^{+2}$
${v}_{{\rm{turb}}}$	<0.16 cs
Tatm0	${85.5}_{-0.8}^{+1.2}$
γ	${0.89}_{-0.02}^{+0.07}$
Tmid0	${22.5}_{-0.6}^{+0.5}$
incl	${48.4}_{-0.4}^{+0.1}$
log(CO/H2)	$-{4.69}_{-0.04}^{+0.05}$
${{\chi }}_{\nu }^{2}$	0.92a

Note. Median plus 3σ ranges for the PDFs defined by the walker positions after removing burn-in.

^aReduced chi-squared from the combination of all four lines for the model defined by the median of the PDFs.

Download table as:  ASCIITypeset image

The multi-line fit demonstrates that one consistent model, with weak turbulence, can fit all four lines. We find a three-sigma upper limit of <0.16 c_s, which corresponds to non-thermal motion of ∼70–130 m s⁻¹ in the upper layers and ∼50 m s⁻¹ at the CO freeze-out temperature. This again falls below the spectral resolution of the data, but the high spatial resolution and high signal to noise help to constrain turbulence to such low values. The weaker constraint on turbulence may be due to the isotopologues, with their lower spatial resolution and increased tolerance for high turbulence models, countering the pull of CO(3-2) toward low turbulence.

When fitting the CO(3-2) line by itself, we kept the CO/H₂ abundance fixed at the ISM value of 10⁻⁴, while our multiline fit finds that the abundance is best fit with a value a factor of ${4.9}_{-0.4}^{+0.6}$ lower. To see if this changes our initial result, we rerun the low-resolution CO(3-2) fit using the new abundance derived from the multiline fit. We find v_turb < 0.02 c_s, consistent with the single line fit.

In fitting just the CO(2-1) line we found that the turbulence measurement could be inflated due to unaccounted for complexities in the temperature structure. A similar effect can be seen here when including parameters for the systematic uncertainty, or when varying disk mass instead of CO abundance. We account for the systematic uncertainty by introducing two additional free parameters, one for the band 6 lines (CO(2-1), ¹³CO(2-1), C¹⁸O(2-1)) and one for the band 7 line (CO(3-2)) since the two bands were calibrated separately. Each parameter is subject to a Gaussian prior with a dispersion of 20%. Accounting for the systematic uncertainty results in a better overall fit to the data (${\chi }_{\nu }^{2}$ = 0.9186 versus 0.9192 when s = 1, ammounting to a >10σ difference) and a stronger constraint on the turbulence (v_turb < 0.13 c_s versus <0.16 c_s when s = 1). Allowing the disk mass to vary, instead of the CO/H₂ ratio, and including systematic uncertainty ends up with a model that is a worse overall fit (${\chi }_{\nu }^{2}$ = 0.9210) and has a weaker constraint on turbulence (v_turb < 0.21 c_s). Going further and removing the systematic uncertainty parameters from the variable disk mass model results in an even worse fit (${\chi }_{\nu }^{2}$ = 0.9226) and an even softer limit on turbulence (v_turb < 0.24 c_s). This progression suggests that as the models get worse at representing the data turbulence gets inflated in an effort to account for the deficiencies in the underlying structure. Such an effect can even be seen in the CO(3-2) low-resolution single line fit, where reducing the CO/H₂ abundance produces a better fit compared to the fiducial model (at >10σ significance) with a stronger limit on the turbulence (<0.02 c_s versus <0.04 c_s for the fiducial fit). While our final best-fit model is certainly not a perfect representation of the HD 163296 system, a better model would likely lead to even more stringent constraints on turbulence.