Weak Turbulence in the HD 163296 Protoplanetary Disk Revealed by ALMA CO Observations
Kevin M. Flaherty, A. Meredith Hughes, Katherine A. Rosenfeld, Sean M. Andrews, Eugene Chiang, Jacob B. Simon, Skylar Kerzner, David J. Wilner
aa r X i v : . [ a s t r o - ph . S R ] O c t Weak Turbulence in the HD 163296 Protoplanetary Disk Revealed by ALMACO Observations
Kevin M. Flaherty , A. Meredith Hughes , Katherine A. Rosenfeld , Sean M. Andrews , EugeneChiang , , Jacob B. Simon , , , Skylar Kerzner , David J. Wilner ABSTRACT
Turbulence can transport angular momentum in protoplanetary disks and influencethe growth and evolution of planets. With spatially and spectrally resolved molecularemission line measurements provided by (sub)millimeter interferometric observations,it is possible to directly measure non-thermal motions in the disk gas that can be at-tributed to this turbulence. We report a new constraint on the turbulence in the diskaround HD 163296, a nearby young A star, determined from ALMA Science Verifi-cation observations of four CO emission lines (the CO(3-2), CO(2-1), CO(2-1), andC O(2-1) transitions). The different optical depths for these lines permit probes ofnon-thermal line-widths at a range of physical conditions (temperature and density)and depths into the disk interior. We derive stringent limits on the non-thermal mo-tions in the upper layers of the outer disk such that any contribution to the line-widthsfrom turbulence is <
3% of the local sound speed. These limits are approximately anorder of magnitude lower than theoretical predictions for full-blown MHD turbulencedriven by the magneto-rotational instability, potentially suggesting that this mechanismis less efficient in the outer (R &
1. Introduction
Planets form within the disks of gas and dust that surround young stars and are subject tothe dynamics within these systems. In particular, magneto-hydrodynamic (MHD) turbulence gen-erated by the magneto-rotational instability (MRI), which is the leading theoretical mechanismdriving angular momentum transport and disk accretion (Balbus & Hawley 1998; Armitage 2011; Van Vleck Observatory, Astronomy Department, Wesleyan University, 96 Foss Hill Drive, Middletown, CT 06459 Harvard-Smithsonian Center for Astrophysics, 60 Garden Street, Cambridge, MA 02138 Department of Earth and Planetary Science, 307 McCone Hall, University of California, Berkeley, CA 94720 Department of Astronomy, 501 Campbell Hall, University of California, Berkeley, CA 94720 Department of Space Studies, Southwest Research Institute, Boulder, CO 80302 Sagan Fellow JILA, University of Colorado and NIST, 440 HCB, Bouler, CO 80309 ∼ z > α ∼ .
01, where α relates the sound speed c s and pressure scale height to viscosity ( ν = αc s H ),are consistent with the observed decrease in accretion rate with time and measurements of diskmass. This method does not directly measure turbulent motion, but instead relies on scaling rela-tions among physically relevant parameters. More direct methods use the motion of the gas itselfand its influence on molecular emission lines, in particular the broadening of the line profile. Linebroadening consistent with transonic turbulent motion has been needed to fit near-infrared COovertone emission, sensitive to the upper layers of the inner AU of the disk, in a number of sources 3 –(Carr et al. 2004; Hartmann et al. 2004; Najita et al. 2009, 1996). Sub-mm interferometric obser-vations are also sensitive to line broadening, but can also trace other effects of turbulence such aschanges in the peak-to-trough ratio of the spectral line profile, as well as broadening of the images(Simon et al. 2015). Resolved images also help break some of the degeneracies with e.g. Keplerianrotation that influence unresolved observations. Hughes et al. (2011) model high spectral resolutionSubmillimeter Array (SMA) observations of CO(3-2) emission and find a tentative (3 σ ) detectionof 300m s − turbulence within the HD 163296 system but only an upper limit ( <
40m s − ) in TWHya. Guilloteau et al. (2012) use CS, a heavier molecule that is less affected than CO by thermalbroadening, to measure 130m s − turbulence within DM Tau, which suggests turbulence of ∼ ⊙ , ∼ CO(2-1), and C O(2-1)) both together and sepa-rately to probe different layers within the disk. In sections 2,3 we describe the data, our model, andour fitting procedure. In sections 4,5 we discuss the results of this fitting, where we find an upperlimit on turbulence that is an order of magnitude lower than predicted by numerical simulations ofMRI, as well as the implications regarding this weak turbulence.
2. The Data
ALMA observed the HD 163296 system as part of early Science Verification operations. Band6 data were taken on 2012 June 9, 23, and July 7th when 24 12m antennas were available coveringbaselines of 20 to 400m. Juno, Neptune, and Mars were used as amplitude calibrators, while thequasar J1924-292 was used as the bandpass calibrator. Science observations were interleaved withmeasurements of the phase calibrator quasar J1733-130 with a total on-source integration time of 84minutes. The spectral coverage was split into four spectral windows. Spectral windows 0 and 3 arecomposed of mostly line-free continuum while spectral window 1 included the CO and C O(2-1)lines and spectral window 2 included the CO(2-1) line. The four spectral windows were centeredon 217.105, 219.969, 230.969 and 233.999 GHz with bandwidths of 1.875, 0.937, 0.937, and 1.875GHz and spectral resolutions of 0.67, 0.33, 0.32, and 0.62 km s − . The beam size is 0.81x0.65”, 4 –while the primary beam has a FWHM of 22”.Band 7 data were taken on 2012 June 9, 11, 22, and July 6 with the same antenna configurationas the band 6 data. Juno and Neptune were used as amplitude calibrators with the quasarsJ1733-130 and J1924-292 used as the phase and bandpass calibrators respectively. The total on-source integration time was 140 minutes. The four spectral windows were centered at 360.169,356.734, 345.796, and 346.998 GHz with bandwidths of 117.3, 468.8, 468.8, 937.5 MHz and spectralresolutions of 0.025, 0.10, 0.10, 0.21 km s − . Spectral windows 0 and 3 cover line-free continuumwhile spectral window 1 contains HCO+(4-3) and spectral window 2 focuses on the CO(3-2) line.We do not include HCO+ in our analysis to simplify the chemical structure that must be accountedfor in our model. The beam size is 0.67x0.43” with a primary beam size of 14”.We use the calibrated Science Verification data; more details on the calibration procedure canbe found on the NRAO website . We do not perform self-calibration which likely leads to an over-estimate of the statistical uncertainties although, as discussed below, our results are dominated bysystematic uncertainties and degneracies betwenn various model parameters rather than statisticaluncertainties. For our analysis we use a low-resolution and high-resolution version of the CO(3-2) line. The low-resolution spectrum has been binned to a resolution of 0.3 km/s and includesonly the central 15 channels. The high-resolution spectrum is at the native resolution of 0.1 km/sand included 101 channels covering the entire line profile. The J=2-1 lines were extracted fromthe calibrated data set, including only the channels with significant line emission (typically ∼ http://almascience.nrao.edu/alma-data/science-verification
3. The Model
To measure the turbulence within the disk, we employ a simple parameteric model to fit thedata. We start with a temperature and surface density structure, compute the vertical structureusing hydrostatic equilibrium, calculate the velocity field taking into account Keplerian rotation aswell as pressure support, and compute the radiative transfer of the line through the disk, accountingfor both thermal and turbulent broadening. This model is based on the work of Rosenfeld et al.(2012) and Rosenfeld et al. (2013), itself based on the structure laid out in Dartois et al. (2003),and has successfully been used to model CO emission from other disks (Rosenfeld et al. 2012;Andrews et al. 2012; Rosenfeld et al. 2013; Williams & Best 2014; Rosenfeld et al. 2014). Wechoose not to derive the gas temperature structure based on a radiative transfer calculation account-ing for stellar irradiation and gas/dust interactions in order to speed up the model calculations andallow for a Markov-Chain Monte-Carlo (MCMC) approach for constraining the model parameters(described below).
In hydrostatic equilibrium, the gas density and temperature are related by: − ∂ ln ρ gas ∂z = ∂ ln T gas ∂z + 1 c s (cid:20) GM ∗ z ( r + z ) / (cid:21) , (1)where ρ gas is the gas density, T gas is the gas temperature, M ∗ is the stellar mass and r, z are theradial and vertical position within the disk. The sound speed is given by c s = k B T gas /µm h where k B is Boltzmann’s constant and m h is the mass of hydrogen. The disk is assumed to be made of80% molecular hydrogen with µ = 2 . T gas ( r, z ) = ( T atm + ( T mid − T atm )(cos πz Z q ) δ if z < Z q T atm if z ≥ Z q , (2), which is based on the Dartois et al. (2003) Type II structure and has been used to successfullymodel the HD 163296 disk CO emission previously. This temperature profile increases from T mid atthe midplane to T atm at a height of Z q . The midplane remains cold because it is shielded from thestellar flux that heats the surface layers. Mechanical heating from the accretion flow can becomeimportant in the midplane of the inner disk (D’Alessio et al. 2006) but we ignore this since ourmeasurements are most sensitive to the outer disk where stellar irradiation dominates. The valuesof T mid , T atm and Z q vary with radius according to: T mid = T mid0 (cid:16) r (cid:17) q T atm = T atm0 (cid:16) r (cid:17) q (3) Z q = Z q (cid:16) r (cid:17) . . This type of profile for the vertical and radial temperature structure is consistent with detailedmodels 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 possibledifferences 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 wetake to follow the form: Σ gas ( r ) = Σ c (cid:18) rR c (cid:19) − γ exp " − (cid:18) rR c (cid:19) − γ . (4)The value of Σ c depends on the total gas mass, the radial power law index of the surface density andthe critical radius (Σ c = M gas (2 − γ ) / (2 πR c )). Unless otherwise specified we assume M gas =0.09M ⊙ (Isella et al. 2007). The critical radius R c defines the turnover between a power law and a sharpexponential decay. This is motivated theoretically by a similarity solution to the disk structurecalculation 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 resolvedsub-mm dust emission from a collection of young solar-type stars and find γ typically close to 1with R c ranging from 15 to 200 AU. This functional form has also been successful in fitting boththe dust and the gas (Hughes et al. 2008) including the different apparent radii for the gas and dustdisk, although recent evidence suggests that in some systems the dust disks are smaller than thegaseous 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 asthe 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 comparisonwith literature models. When used we define H as c s / Ω, where the midplane temperature sets thesound speed.The velocity field is taken to be Keplerian, with modifications due to the height above themidplane and the pressure support of the gas: v r = GM ∗ r ( r + z ) / + 1 ρ gas ∂P gas ∂r , (5)where ρ gas and P gas and the gas mass density and pressure respectively. Rosenfeld et al. (2013)estimate that these modifications change the velocities by a few percent in the upper layers of thedisk. We do not include the self-gravity of the disk since this has a negligible effect on the velocityfield, especially for our small gas mass (Rosenfeld et al. 2013). 7 –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 emissionfrom the HD 163296 system, Qi et al. (2011) find a freeze out temperature of 19 K, which we takeas the freeze-out boundary in our models. We simulate freeze-out by decreasing the CO abundanceby 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 is0 . Z ∞ z phot n gas ( r, z ) dz < σ s . (6)This is the same boundary used in Qi et al. (2011) with σ s =0.79 × × cm − . This boundaryis defined by the absorption of UV photons from the central source (Aikawa & Herbst 1999) andis included in our model as a decrease in the abundance of CO by eight orders of magnitude. Theexact location of the photo-dissociation boundary may vary between different isotopologues due toselective self-shielding (Visser et al. 2009; Miotello et al. 2014). When modeling CO and C Oemission we initially assume ISM values of C/ C = 69 and O/ O = 557 (Wilson 1999), al-though 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 chemicalmodel for CO ignores much of the complexity of disk chemistry (Walsh et al. 2010; Jonkheid et al.2007) but has proven successful in fitting the CO emission from the HD 163296 system in the past(Rosenfeld et al. 2013; Qi et al. 2011). Once the structure of the disk is set, we calculate the flux through the disk. The line intensityis: I ν = Z ∞ S ν ( s ) exp[ − τ ν ( s )] K ν ( s ) ds, (7)where s is the linear coordinate along the line of sight increasing outward from the observer. Theoptical depth is τ ν ( s ) = R s K ν ( s ′ ) ds ′ where K ν ( s ) is the absorption coefficient and S ν ( s ) is thesource 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 temperatureand the source function is approximated by the Planck function. The CO emitting region is ata density that is much higher than the critical density ( ∼ cm − ), allowing us to assume LTEwithout much loss of accuracy (Pavlyuchenkov et al. 2007)The line is assumed to be a Gaussian with width∆ V = q (2 k B T ( r, z ) /m CO ) + v , (8)and a FWHM of 2 √ log 2∆ V . Throughout much of our modeling we scale v turb to the local soundspeed (e.g. ∆ V = p (2 k B T ( r, z ) /m CO ) + ( v turb /c s ) ). The sound speed, which is a factor of 8 – p m CO / ( µm H ) ≈ v turb , associated with the eddies but instead measurethe 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 meanvelocity are the same to within factors close to unity. Given the similarity of these factors, and theprevalence of the notation v turb in the literature, we use v turb to denote turbulence even thoughstrictly speaking we are measuring the RMS width, δv turb , of the velocity distribution.Finally, we keep the distance (d=122pc), 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 notstrongly depend on the other model parameters, and were adjusted to minimize residuals usinginitial model fits. Using the model described above, we can generate synthetic images for a given set of param-eters. We use the MIRIAD task UVMODEL to generate model visibilities sampled at the samespatial frequencies as the data. Even though the data contain both XX and YY polarizations, weonly use the total intensity (I=(XX+YY)/2) to compare with the model because the line emissionis 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 isderived from the dispersion in the data, as described earlier.To sample the posterior probability distribution for each of the parameters we use a MCMCtechnique. 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 toa Metropolis-Hastings chain, the affine-invariant method uses a large number of walkers whosemovements 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 distributionfunction (PDF) and can efficiently probe the PDF even when there are strong degeneracies betweenparameters, as is sometimes the case with our models. After an initial burn-in period the positionsof the walkers along the chains represent independent samples from the PDF and the distributionof the walkers in parameters space can be used as an estimate of the PDF. Our chains typicallyconsist of 40 walkers and 3000 steps. We only examine the last 1000-2000 steps in each chain,corresponding to >
10 autocorrelation times, resulting in 40000-80000 samples of the PDF. Theinitial positions of the walkers are chosen to be evenly distributed throughout parameter space. Werequire 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 9 –the results. ALMA line fluxes differ from those previously measured by the SMA for these samelines (Qi et al. 2011; Hughes et al. 2011) by ∼ τ =1 surface and a systematic uncertainty in thetotal line flux will directly translate into an uncertainty in the temperature. Since both temperatureand turbulence contribute to broadening of the line, an uncertainty in one will translate into anuncertainty in the other (although, as discussed below and in Simon et al. (2015), the spatialresolution of the images and the peak-to-trough ratio of the line profile help to break some of thisdegeneracy). For the optically thin lines whose flux is proportional to both the surface density andthe temperature, the exact influence of the systematic uncertainty is less obvious. This is especiallyimportant when combining information from multiple lines, which are subject to different systematicuncertainties.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 beingbrighter than the data (e.g. s = 1 . s =1, but for CO(3-2) we consider fitting when s is fixedat 0.8 or 1.2, and we treat it as a free parameter when fitting the four lines simultaneously. Despitethe large number of degrees of freedom, ∼ millions for the single line fits, including too many freeparameters prevents the walkers from converging on a best fit model, especially when there arestrong degeneracies involved. For this reason we fix s in the single line fits. In the multi-linefit there is enough additional information to allow for a constraint on s , and there we considersituations where s is left as a free parameter. The total number of free parameters depends on theline being fit and ranges from four ( v turb , γ , inclination, N(X)/N( C O)) when fitting CO andC O to 10 ( q , R c , v turb , T atm0 , γ , T mid0 , inclination, log(CO/H ), s and s ) when fitting all fouremissions lines simultaneously while accounting for the systematic uncertainty.
4. Results4.1. CO(3-2)
We start by modeling the high-resolution CO(3-2) line using 40 walkers over 1500 steps. Thefirst 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.09M ⊙ and γ =1, consistent with previous models (Isella et al. 2007; Qi et al. 2011). We allow R c to vary tocontrol the radial extent of the emission. By spatially resolving the front ( z >
0) and back ( z < T mid , T atm , and 10 – Z q ). Our initial attempts at MCMC modeling found that T atm0 and Z q were highly degenerate,which slowed significantly the convergence of the walkers. Given the preference of the models forhigher Z q we fix Z q =70AU. With many resolution elements across the radius of the disk, we canconstrain the radial profile of the temperature, and allow q to vary. Finally, we include turbulenceas 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. ∆ I V = ( V obs ( u, v ) − V mod ( u, v ))) in Figure 2. This modelreproduces much of the morphology of the disk, including the front/back emission with a drop inemission at the cold midplane, and has a reduced chi-squared close to one. There are still someresiduals between the model and the data but these may be due to deviations from our simplefunctional form for the temperature and density structure, or even non-circular motion, that willnot be perfectly captured in our model. A fit directly to the low-resolution data returns a nearlyidentical, although slightly weaker limit on turbulence (Table 1,Figure 3), while also producing agood fit to the data. Slight variations in the temperature structure are consistent with degeneraciesand are well below the systematic uncertainties, discussed below. The similarity in the fit to thesespectra with different velocity resolutions indicates that the spatial information is providing muchof the constraint on the models. Given the consistency between the high and low-resolution dataresults, while we report the structure from the high-resolution PDFs, we explore various aspects ofthe fitting using the low-resolution spectra.The temperature and density structure for the best fit model are shown in Figure 4. TheCO is confined to a thin layer that, at 200 AU, stretches from 10 to 60AU above the midplanewith temperatures ranging from 19 K 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 transfercalculations of the equilibrium disk temperature, assuming a dust distribution and a luminosity forthe central source, and found similar temperatures to those in our models. This is not surprisingsince CO(3-2) is optically thick and strongly sensitive to the temperature. Our radial temperatureprofile falls off with a power law index of -0.278 +0 . − . (-0.216 +0 . − . for the low-resolution fit),significantly shallower than the canonical value of -0.5, but consistent with values expected fromradiative 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 atwhich the temperature falls off with distance from the central star.While the statistical errors on the model parameters are small, they do not account for thesystematic effects. For example, the uncertainty in distance strongly influences q since a moredistant disk requires more extended radial emission which can be accomplished with a flatter radialtemperature profile. Adjusting the distance by ± q varies by up to 20% ; q risesfrom -0.236 ± .
004 to -0.175 +0 . − . as distance increases, compared to -0.216 +0 . − . at the fiducialdistance. The critical radius also increases with distance, from 158 +6 − to 203.7 +19 − AU for the 11 –near and far distances respectively, while the other model parameters are consistent within thestatistical uncertainty. The differences between q and R c when fitting the high-resolution andlow-resolution data (Table 1) also highlight the degeneracy between these two parameters. Themaximum spatial extent of the CO emission can be matched with a narrow disk (low R c ) andshallow 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 tothe statistical uncertainties.Similarly, the statistical uncertainty on the atmosphere temperature is <
5% although this doesnot include the uncertainty on the amplitude calibration. As discussed earlier, we account for thesystematic uncertainty by applying a scale factor to the model and we run additional MCMC trialsfitting 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 +0 . − . to 101 ± ± T atm0 . The variations ofthe other parameters in response to the uncertainty in the flux calibration are within the statisticalerrors 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 atmospheretemperature, the midplane temperature is constrained through the emission from the faint lowerhalf of the disk and the lack of emission between these two features. In our fiducial model we find T mid0 =21.8 +0 . − . . The constraints on q and T mid0 in turn lead to a CO ice line that falls withinthe range R ice =245 +27 − AU while previous observations of chemical tracers of CO freeze-out findthat it is closer to 150-160AU (Qi et al. 2013b; Mathews et al. 2013). Overestimating the midplanetemperature results in a more puffed up disk, pushing the τ = 1 surface higher in the disk and theatmosphere 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 spectrawhich has different values of T mid0 and T atm0 , and predicts a CO ice line at R ice =103 +23 − AU, butproduces similar images.While our models do not treat the dust, assumptions about the dust distribution may still beimplicit in our models. Jonkheid et al. (2007) find that in their chemical models the depletion ofdust from the upper layers of the dust leads to a drop in the gas temperature. In the context ofthe D’Alessio et al. (2006) models, Qi et al. (2011) show that the settling of large dust grains willaffect the height at which the disk reaches the atmosphere temperature; in other words more settlingresults in lower Z q . Our assumed value of Z q is most similar to models with little dust settling(Qi et al. 2011). The significant degeneracy between Z q and T atm0 means that a model with lower Z q , when combined with a lower T atm0 , would fit the data as well as the models presented here.Rosenfeld et al. (2013), using the same functional form for the vertical temperature profile as we dohere, find that the data can be fit with T atm0 =64K, given Z q =43 AU. This atmosphere temperatureis much lower than our T atm0 =86 ±
1K but the difference in Z q means that both models produce a 12 –nearly identical gas temperature at the τ = 1 surface. Further data probing the uppermost layersin 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 previousradiative transfer models of the HD 163296 system. By carefully accounting for parameter degen-eracies and systematic effects due to assumptions about the distance and flux calibration, we canmeasure the statistical and systematic uncertainties in temperature and density structure of thedisk. This indicates that our derived disk structure is not going to strongly bias our characterizationof turbulence. When fitting either the high-resolution or low-resolution CO(3-2) spectra, we consistently findlow levels of turbulence (Table 1). We conservatively quote three-sigma upper limits of v turb < s and v turb < s respectively. In the uppermost layers of the CO region, these limits correspondto velocity dispersions of less than ∼ − , with higher velocities at smaller radii. While wehave some constraint on the inner disk from the high velocity channels, most of the informationabout the turbulence comes from the spatially resolved outer disk ( R > δ v=0.1km s − ). To test whether we can distinguish betweenturbulence at the resolution limit from much weaker non-thermal motion, we run an additionalMCMC fit to the low-resolution CO(3-2) data with v turb fixed at 0.1 km s − while allowing theother model parameters to vary. Figure 6 shows the best fit v turb = 0 . − model along withour fiducial low turbulence model, using a Gaussian fit to the short-baseline visbilities to derivethe 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 betterfit ( > σ significance) from the low-turbulence models. This behavior indicates that fitting a modelto the full data set has a stronger diagnostic power than the broadening of the line in the spectraldomain. Simon et al. (2015) have suggested the peak-to-trough ratio as a potential diagnosticof turbulence. They found in simulated observations of numerical MRI simulations, with typical v turb ∼ . − . s at the CO(3-2) emitting surface, that the peak-to-trough ratio decreases as theturbulence increases. This behavior can be seen in Figure 6 where the v turb = 0 . − modelshave a much lower peak-to-trough ratio than the low-turbulence models and the data, consistentwith the conclusion that the turbulence is weaker than 0.1km s − .Peak-to-trough ratio is a one-dimensional metric that provides a convenient way of diagnosingdifferences in the three-dimensional data set. Variations in the temperature can also affect thepeak-to-trough ratio, however the uncertainty in temperature is dominated by the uncertainty inthe 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 13 –full three-dimensional data set, we account for the various influences of the different parameterson the images as well as the spectral shape. The similarity between the turbulence results whenfitting the low-resolution and high-resolution data suggest that the shape of the full line profileis not the dominant factor but instead spatial resolution plays a large role in constraining thenon-thermal motion. In the individual channel maps turbulence broadens the width of the diskemission (Guilloteau et al. 2012; Simon et al. 2015). We demonstrate this effect for the centralvelocity channel of the HD 163296 system in Figure 7 by varying the turbulence while fixing theother model parameters. For low levels of turbulence there is little variation in the images, butthe disk broadens dramatically for larger values of turbulence. This occurs even below the spectralresolution of the data, allowing us to push to low levels of turbulence.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 thisway line flux could be used as a diagnostic of turbulence, but this leaves it highly degenerate withtemperature, and subject to the uncertainty in the calibration of the line flux. The systematicuncertainty 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 ascale factor to the model. When fixing s = 1 .
2, i.e. when the true disk flux is 20% brighter than theobserved data, we find a limit of v turb < . s , while if the true disk emission were 20% fainter thanour data the limit on turbulence is v turb < . s . Most of the change in disk flux is absorbed by theatmosphere temperature, with only small variations in turbulence. An uncertainty in distance couldalso affect our result, since the one-sigma uncertainty on distance of 10pc (Montesinos et al. 2009)leads to a 16% uncertainty in the total flux. In additional MCMC trials with distance fixed at 112and 132pc we find that the limit on turbulence increases to v turb < s ,0.06c s respectively. Evenaccounting for the systematic uncertainty and the distance uncertainty, we find that non-thermalmotion is limited to very weak levels.The high-resolution spectrum exhibits an asymmetry in the line profile, with the blue-shiftedpeak brighter than the red-shifted peak (Figure 6). One possible explanation for the behavior iswith an eccentric disk; the disk material ’piles-up’ on the slower moving side of the orbit leadingto a higher flux from one side of the disk versus the other (e.g. Reg´aly et al. 2011). Our modelincludes the dynamics of Keplerian rotation and thermal motion and any additional motion willlikely be accounted for by the turbulence parameter. For this reason we choose to report our tur-bulence 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 featuresbeing modeled as turbulence, as well as the weak influence of turbulence on the data at such lowlevels (Figure 7), we conservatively quote an upper limit. Also, as noted earlier, or implementationof turbulence within the models means that we are really constraining the velocity dispersion, orRMS velocity, associated with turbulence rather than the mean velocity within turbulent eddies. 14 – As with CO(3-2), our CO(2-1) MCMC trial uses 40 walkers over 3000 steps. With a lowerspatial resolution than the CO(3-2) data, the walkers here take longer to converge and the numberof steps has been increased accordingly. CO(2-1) is also optically thick, making it sensitive tothe temperature structure, but lacks the spatial resolution to distinguish the near and far side ofthe disk. For this reason we fix T mid0 =17.5, based on the low-resolution CO(3-2) results, whileallowing q , R c , v turb , T atm0 , and inclination to vary during the fitting. Results, after removing thefirst 1000 steps, are listed in Table 1 with posterior distributions shown in Figure 8. Channel mapsand 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 fluxthroughout the spectra.Differences in the model parameter PDFs exist between the CO(2-1) and CO(3-2) models, butthese are most likely the result of degeneracies between parameters. In the CO(2-1) data we finda strong anti-correlation between R c and q and compared to the low-resolution CO(3-2) fit, theCO(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 thetwo lines still has very similar temperatures at the top of the CO emitting surface, to which thesedata are particularly sensitive.The upper limit on turbulence from CO(2-1) ( v turb < . s ) is substantially higher thanderived from the CO(3-2) line, most likely due to the lower spatial resolution of these data comparedto CO(3-2). This limit corresponds to velocity widths of 200-150 m s − in the upper layers of thedisk. To see if a satisfactory fit can be found with low turbulence, we run an additional trial withv turb fixed at 0.04c s , ie. the upper limit determined by the low-resolution CO(3-2) data. Theresulting 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 showfeatures only on the north-east side of the disk (Figure 12). This asymmetry may be related tothe temperature structure of the disk, as probed by an optically thick tracer. As discussed inDartois et al. (2003), the line of sight to an inclined disk probes different heights for the area ofthe disk closer to the observer, in this case the south-west side, than on the far (north-east) side ofthe 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 lineof sight through the far side of the disk has its τ =1 surface in an area with a lower density andtemperature than where it entered the disk. While our ray-tracing code will account for photonscoming from different surfaces, our assumed functional form for the temperature may not capturethe 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 cannotadjust 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 15 –of the disk. Since the CO(2-1) data has a lower spatial resolution than CO(3-2), it has more roomto increase the broadening of the images before the model becomes significantly wider than thedata. The end result is a high value of turbulence, and a worse fit to the peak-to-trough ratio, dueto limitations of the assumed temperature structure rather than large non-thermal motions withinthe disk. For this reason we quote an upper limit on turbulence despite the shape of the posteriordistribution (Figure 8) CO(2-1)
The decreased abundance of C relative to C reduces the optical depth of CO relativeto CO which makes the line emission more sensitive to the total mass of the disk, but intro-duces a degeneracy between column density and temperature. To avoid this degeneracy we fixthe 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 MCMCmodeling, and in the multiline fits discussed below, we found the best possible fit to the data camewhen 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 toallowing CO/H to vary. Model fits, after removing the burn-in, are listed in Table 1 with posteriordistributions shown in Figure 13. These models are able to accurately reproduce much of the COemission, although there are some discrepancies near the line peaks (Figure 14). This is likey dueto 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 themodel and the data are small, with the strongest residuals only 10% of the peak flux.As with CO(2-1), we put a limit on the turbulent broadening ( v turb < s ) that is lessstrigent than CO(3-2). At the coldest layers close to the midplane, just above CO freeze-out, thiscorresponds to a non-thermal velocity width less than ∼ − . We also find CO/CO=233 +3 − ,significantly different from the ISM value of 69 (Wilson 1999). This may be a sign of unaccountedfor CO chemistry, such as selective photodissociation (Miotello et al. 2014), or a sign that theassumed disk mass, which is derived from the dust emission, is an overestimate. This result isdiscussed in more detail below in the context of the multiline fit. O(2-1) C O, being more heavily depleted than CO, is fully optically thin and is able to trace theentire vertical extent of the disk. As with CO, we break the degeneracy between density andtemperature by using the derived values of q , T mid0 , and T atm0 from CO(3-2), while allowing thedepletion and γ to vary. The resulting PDFs are shown in Figure 15 with channel maps and residuals 16 –in Figure 16. We find a significant degeneracy between the depletion factor and γ which is notsurprising since both control the surface density of the disk. This anti-correlation may contributeto the difference in γ between the fit to C O and the fit to CO; C O implies a steeper surfacedensity profile, but this can be brought in line with the CO result with a depletion factor of ∼ O line is fit with a higher inclination than the other lines. This difference issmall, but significant. Further work is needed to fully explore the nature of this discrepancy.The C O line, with the lowest signal to noise, places the weakest constraint on the turbulentbroadening, v turb < . s . At the coldest layers just above the CO freeze-out temperature, thisimplies random velocity dispersions of less than 150m s − . To check if a low turbulence modelcan adequately fit the data we ran an additional MCMC trial with the turbulence fixed at thelow-resolution CO(3-2) limit ( v turb =0.04c s ) and the results are shown in Figure 17. The differencebetween the high turbulence fit and the low turbulence fit is small despite the order of magnitudedifference in turbulence. While C O is more sensitive to the midplane, its optically thin naturemakes it less sensitive to velocity broadening from turbulence. In an optically thick line, the velocityshear provided by turbulence allows more photons to escape, increasing the flux, while in opticallythin lines this effect is much smaller (Horne & Marsh 1986). This is consistent with the behaviorseen in Figure 17 for C O, as well as the weak limit on turbulence in CO.
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 ∼ CO and C O probe materialcloser to the midplane. By fitting one model to all four lines we can leverage this complementaryinformation to more accurately constrain the temperature, density, and turbulence throughout thedisk. For this MCMC trial we employ 60 walkers over 1000 steps, with the first 500 removed asburn-in. We fix M gas =0.09M ⊙ but allow q , R c , v turb , γ , T atm0 , T mid0 , and inclination to vary, alongwith 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 reportthose results here. Numerical simulations of MRI predict strong vertical gradients in the turbulentmotion as a function of height within the disk (Fromang & Nelson 2006; Simon et al. 2013) butgiven the lack of detection found with the other lines we do not attempt to include any change inturbulence with height, beyond allowing it to vary with the local sound speed.The results are listed in Table 2 with visibility spectra shown in Figure 19, and maps of thecentral velocity channel for each line in Figure 20. In terms of the overall density and temperaturestructure, while some of the parameters change compared to where they stood with the individualline fits, the overall structure is similar. Figure 21 shows the density and temperature structurefor 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 . − . relative to the canonical ISM value. 17 –The individual fits to CO and C O suggest a depletion of only a factor of ∼
3, which is smallerthan that derived here because of the lower midplane temperature assumed in those fits. Recentobservations 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 theCO/H ratio (Furuya & Aikawa 2014; Reboussin et al. 2015). It is also possible that the deviationin CO/H is a sign that the assumed disk mass is incorrect. The mass was derived from the sub-mmemission of the dust (Isella et al. 2007) and could deviate from the true value if there is a changein the gas to dust mass ratio, a concentration of dust into very large grains that are invisible inthe sub-mm, or if the dust does not follow the gas distribution. We have run additional trialsfixing the CO/H ratio while allowing the disk mass to vary and find M gas =0.08 ± ⊙ , which issimilar to our assumed value, but is subject to large degeneracies, especially with R c and producesa worse overall fit to the data. Further work, including a more physically realistic treatment of COchemistry near the ice line, is needed to confirm and understand the low CO/H ratio.The multi-line fit demonstrates that one consistent model, with weak turbulence, can fit allfour lines. We find a three-sigma upper limit of < s , which corresponds to non-thermal motionof ∼ − in the upper layers and ∼
50 m s − at the CO freeze-out temperature. This againfalls below the spectral resolution of the data, but the high spatial resolution and high signal tonoise help to constrain turbulence to such low values. The weaker constraint on turbulence maybe due to the isotopologues, with their lower spatial resolution and increased tolerance for highturbulence models, countering the pull of CO(3-2) towards low turbulence.When fitting the CO(3-2) line by itself, we kept the CO/H abundance fixed at the ISM valueof 10 − , while our multiline fit finds that the abundance is best fit with a value a factor of 4.9 +0 . − . lower. To see if this changes our initial result, we rerun the low-resolution CO(3-2) fit using thenew abundance derived from the multiline fit. We find v turb < s , consistent with the single linefit. In fitting just the CO(2-1) line we found that the turbulence measurement could be inflateddue to unaccounted for complexities in the temperature structure. A similar effect can be seenhere when including parameters for the systematic uncertainty, or when varying disk mass insteadof CO abundance. We account for the systematic uncertainty by introducing two additional freeparameters, 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 Gaussianprior with a dispersion of 20%. Accounting for the systematic uncertainty results in a better overallfit to the data ( χ ν =0.9186 vs 0.9192 when s=1, ammounting to a > σ difference) and a strongerconstraint on the turbulence ( v turb < s vs < 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 aworse overall fit ( χ ν = 0.9210) and has a weaker constraint on turbulence ( v turb < s ). Goingfurther and removing the systematic uncertainty parameters from the variable disk mass modelresults in an even worse fit ( χ ν = 0.9226) and an even softer limit on turbulence ( v turb < s ).This progression suggests that as the models get worse at representing the data turbulence gets 18 –inflated in an effort to account for the deficiencies in the underlying structure. Such an effect caneven be seen in the CO(3-2) low-resolution single line fit, where reducing the CO/H abundanceproduces a better fit compared to the fiducial model (at > σ significance) with a stronger limit onthe turbulence ( < s vs < s for the fiducial fit). While our final best-fit model is certainlynot a perfect representation of the HD 163296 system, a better model would likely lead to evenmore stringent constraints on turbulence.
5. Implications of weak turbulence
The high(low)-resolution CO(3-2) line puts an upper limit of only 3.1%c s (3.8%c s ) on theturbulence velocity dispersion in the outer disk, with a small increase to ∼ CO(2-1), and C O(2-1) linesoffer less stringent upper limits ( < s , < s , < s respectively) due to their lower spatialresolution. Simple tests fitting CO(2-1) and C O(2-1) with turbulence fixed at 3.8%c s show thatboth lines can be well fit with low turbulence models, which we also found in a combined multilinefit that limited turbulence to <
16% c s . In the case of optically thick lines like CO(2-1) subtletiesin the temperature structure can be masked by turbulence and in the case of optically thin lineslike C O(2-1) the profile is not particularly sensitive to turbulence. The upper limit from CO(3-2) corresponds to ∼ − across the upper layers of the CO emitting region while the moreoptically thin CO and C O limit non-thermal motion close to the CO freeze-out zone at themidplane to < − , although the multiline fit pushes that limit down to <
50m s − . While thisfalls below the spectral resolution of the data, fitting the full three-dimensional data set includesthe spatial information that drives the limit to such low values.The viscosity associated with turbulence is often parameterized by α which encapsulates muchof the detailed physics relating the largest velocity and spatial scales, c s , and H, to the turbulentviscosity ( ν = αc s H ) (Shakura & Sunyaev 1973). In a turbulent cascade, energy is deposited onthe largest spatial scales and cascades down to smaller scales (Kolmogorov 1941). Even though atthe spatial resolution of our CO(3-2) data ( ∼ ∼ α entails the strength and size of the cascade relative to H and c s which represent the largestspatial and velocities scales for the turbulence. The parameter α can also be related directly tothe stress tensor W Rφ , which is the density weighted sum of the Maxwell and Reynolds stresses(Balbus & Hawley 1998). This tensor encapsulates the correlated fluctuations in the R- φ plane andis connected to the viscous evolution through α =W Rφ /c s . The form of this equation suggests that α ∼ ( v turb /c s ) , consistent with simulations that measure both turbulent motion and the stresses(e.g. Simon et al. 2013, 2015). With this conversion from turbulent velocity to α , our CO(3-2)observations correspond to an upper limit of α < . × − while the multiline fit provides a limitof α < . σ detection of 300m s − turbulence for the HD 163296 system based on SMA data with higher spectral, but lower spatial,resolution, while de Gregorio-Monsalvo et al. (2013) find that a model with 0.1km s − turbulencefits the ALMA dataset better than a model with no turbulence. Our CO(3-2) limit correspondsto ∼ − across the upper layers of the CO emitting region. While the best fit model inHughes et al. (2011) is able to fit the SMA data, the limited constraint on the temperature structuremeans that this data is subject to a stronger degeneracy between temperature and turbulence. Inparticular, the low spatial resolution SMA data has difficulty distinguishing between a cold, flatdisk with large turbulent broadening and a warm, puffy disk with small turbulent broadening. Byseparating the front and back side of the disk we can limit the models to warmer disks with weakerturbulence.The difference between our upper limit and de Gregorio-Monsalvo et al. (2013) likely is a resultof different disk models. They employ a radiative transfer code to calculate the dust temperatureand then use this to calculate the gas emission subject to assumptions about equal gas and dusttemperatures, constant gas to dust mass ratio, the grain size distribution and the abundance of COrelative to H , subject to freeze-out and photodissociation. Our disk models focus solely on the gasstructure making them simpler although at the expense of the additional information about diskstructure that comes from the dust emission. The dust will mainly influence the gas temperature;the dust dominates the opacity and absorbs thermal energy from stellar radiation transferring itthrough collisions to the gas. By assuming a functional form for the temperature structure we avoidthe complicated radiative transfer problem relating the dust and stellar emission. Efforts to use thedust to constrain the gas structure have been complicated by recent observations that have foundevidence of different gas and dust distributions in protoplanetary disks. Andrews et al. (2012) findthat the dust disk within TW Hya is much more compact than the gas, even when accounting fordifferences in optical depth. de Gregorio-Monsalvo et al. (2013) note that their fit to the gas anddust in the HD 163296 system requires different surface density gradients. This suggests that amore detailed treatment is needed to fully reconcile the dust and gas emission. The complexity ofthe models used in de Gregorio-Monsalvo et al. (2013) do not allow for a full characterization ofthe uncertainties on the various parameters, which allows for the possibility that our results areconsistent with theirs to within the uncertainties.Our limits also fall below measurements in other stars. Hartmann et al. (1998) find α ∼ . − , which, at the typical height of the CS emission(z/r ∼ CO and C O ( . − ) probing a similar areawithin the disk, although CS may be subject to the same optically thin sensitivity issues as C O.Numerical simulations of the MRI, with background fields optimized for full-blown MRI tur-bulence in far-ultraviolet ionized layers, find that turbulent motions vary from a few percent of thesound speed near the midplane up to values comparable to the sound speed at 3-5 scale heightsabove the midplane, at surface densities . − (Simon et al. 2013; Fromang & Nelson 2006;Perez-Becker & Chiang 2011; Simon et al. 2015). The CO(3-2) τ = 1 surface in our model diskslies at z=30-100AU (Figure 18), which corresponds to 3H far from the star and 5H close to thestar, or equivalently surface densities ∼ × − g cm − . Simon et al. (2015), in their numericalsimulations of a disk modeled after the HD 163296 system, find that at these heights/surface den-sities at 100 AU, far-ultraviolet (FUV) photons are able to effectively ionize the disk leading tostrong turbulence. The distribution of velocities peaks at v turb /c s ∼ ∼ v turb /c s ∼ τ = 1 C O surface pierces through the disk toreach z=H on the bottom half in the disk, making it more sensitive to the midplane than CO(3-2).The less stringent limit provided by this line is consistent with the expected weak turbulence nearthe midplane, although as discussed earlier optically thin lines have a weaker response to velocityshear, limiting the ability of C O to constrain turbulence. Simulations have also found that a netvertical magnetic field is required to produce strong turbulence (Bai & Stone 2011; Simon et al.2013); in the HD 163296 system the MRI may still operate but simply with a weak magneticfield that results in weak turbulence. An MRI ‘dead zone’ due to Ohmic diffusion may be presentnear the midplane where the ionization is insufficient to effectively couple the gas to the magneticfield, although these zones are expected to extend out to R ∼ − R > × − M ⊙
21 –yr − (Mendigutia et al. 2013). Assuming this accretion rate applies throughout the disk, and thatthere are no external torques, α is directly related to ˙ M and the temperature and density structureof the disk (Armitage 2011). This accretion rate, combined with a disk mass of 0.09M ⊙ and acritical radius of 175AU, implies α increasing from 0.06 at 10 AU to 0.3 at 500 AU while ourobservations suggest α is less than 9.6 × − . The low disk mass and low turbulence appear to beincompatible with the high accretion rate onto the star. An uncertainty in the disk mass will factorinto this estimate of α , but the disk mass would need to be ∼ ⊙ to bring α in line with ourmeasurements. One possible explanation is that the viscosity increases strongly toward the innerdisk, leading to a higher accretion rate onto the star than in the outer disk. At the low densitiesin the outer disk ambipolar diffusion dominates while closer to the star ohmic dissipation and thehall effect play a larger role (Armitage 2011; Turner et al. 2014). It may also be the case that thedisk is not in steady state equilibrium. Mendigutia et al. (2013) find that the current accretionrate is at least an order of magnitude higher than in observations taken 15 years earlier. Given adisk mass of 0.09M ⊙ , an accretion rate of 5 × − M ⊙ yr − , two orders of magnitude lower thanthe recent observations, could reproduce the observed α in the outer disk.Another explanation for the high accretion rate with low turbulence is if there is a form ofangular momentum removal that does not lead to an observable velocity dispersion, such as a diskwind (e.g. Blandford & Payne 1982). Numerical models find that magneto-hydrodynamics windscan be driven from the disk, although the exact results depend strongly on the numerical setup(Fromang et al. 2013). Klaassen et al. (2013) observe evidence for a disk wind from HD 163296based on the morphology of large scale CO emission. Further observations are needed to determineif this wind originates in the outer disk, and if it removes enough angular momentum to explainthe accretion rate.Gravito-turbulence has been suggested as an alternative to MRI as the driver of turbulencewithin the disk (Armitage 2011), although as with MRI the predicted magnitude is still largerthan our upper limits. Shi & Chiang (2014) predict v turb /c s ∼ turb /c s ∼ α ∼ − (Nelson et al. 2013). Other hydrodynamicinstabilities (e.g. Lyra et al. 2014; Marcus et al. 2014) may also play a role when MRI is weak, 22 –depending on the exact structure of the disk. Further modeling is needed to determine if thephysical conditions necessary for these instabilities are present in this disk.Such low turbulence can lead to measurable effects on the dust distribution and the chemistry inthe disk. The dust vertical density distribution in a turbulent disk is particularly sensitive to the gasvelocities in the upper layers of the atmosphere (Fromang & Nelson 2009) with diminished velocitiesleading to increased settling (Dubrulle et al. 1995; Johansen & Klahr 2005; Carballido et al. 2006;Ciesla 2007). Our CO(3-2) measurements find weak turbulent motion in the upper atmosphere,suggesting that dust settling should be significant and indeed HD 163296 has been observed to havea diminished mid-infrared flux consistent with substantial settling of small dust grains (Meeus et al.2001; Juh´asz et al. 2010). The chemical structure may contain imprints of turbulence since largeturbulent velocities lead to enhanced mixing and a smoothing of chemical gradients. This canbe especially important near sharp boundaries, such as ice lines. Given that turbulence operatesover a typical length scale of H, it more strongly affects the steeper vertical chemical gradientsthan the shallower radial gradients (Semenov & Wiebe 2011) and is unlikely to strongly affectthe radial location of e.g. the CO snow line. Furuya & Aikawa (2014) find that models withweak turbulence can lead to substantial depletion of CO, since a more laminar structure allowsfor the slow sink of carbon out of CO and into carbon-chain molecules. This is consistent withour evidence of diminished CO/H although more detailed modeling is needed to distinguish thiscomplex chemistry from other effects such as selective photodissociation.
6. Conclusions
Using ALMA science verification data of CO(3-2), CO(2-1), CO(2-1), and C O(2-1) we con-strain the turbulent velocity dispersion within the protoplanetary disk surrounding the young starHD 163296. By employing physically realistic models and a Bayesian analysis we find a limit fromthe high-resolution CO(3-2) data of v turb < s , which implies α < × − , with consistentresults when fitting the other lines individually, or together with the CO(3-2) data. This tightconstraint comes from the high spatial resolution, which can both probe the temperature structureand directly constrain the turbulence in a way that was not possible with previous instruments.Our upper limit is well below that expected from the strong accretion onto the star itself, suggestingeither that angular momentum transport is strongly radially dependent or that angular momentumis removed in a way that does not generate significant non-thermal motion. Our limits also fallbelow predictions from numerical simulations of the MRI, suggesting that MRI is less efficient inthe outer disk than previously thought. This demonstrates that with the high spatial resolutionand sensitivity of ALMA we can place physically meaningful constraints on the turbulence.We thank the referee for comments that helped improve the quality of the paper. We gratefullyacknowledge support from NASA Origins of Solar Systems Grant NNX13AI32G. We also thankChunhua Qi and Karin Oberg for helpful discussions regarding the structure of CO. This paper 23 –makes use of the following ALMA data: ADS/JAO.ALMA Table 1. Model Fitting ResultsLine q R c (AU) v turb /c s T atm0 (K) γ T mid0 (K) inclination N(X)/N( C O) χ ν a CO(3-2) (high-res) -0.278 +0 . − . +13 − < ± b +0 . − . ± . +0 . − . +12 − < ± b +0 . − . ± . +0 . − . +21 − < +1 . − . b b +0 . − . CO(2-1) -0.216 b b < b +0 . − . b ± . +3 − O(2-1) -0.216 b b < b ± .
08 17.5 b +0 . − . +120 − σ ranges for the PDFs defined by the walker positions after removing burn-in. a Reduced chi-squared for the model defined by the median of the PDFs b Held fixed during the model fitting. 25 –Table 2. Multi-line model fitParameter Result q -.298 +0 . − . R c +2 − v turb < . s T atm0 +1 . − . γ +0 . − . T mid0 +0 . − . incl 48.4 +0 . − . log(CO/H ) -4.69 +0 . − . χ ν a Note. — Median plus 3 σ ranges for the PDFs definedby the walker positions af-ter removing burn-in. a Reduced chi-squaredfrom the combination ofall four lines for the modeldefined by the median ofthe PDFs. 26 – −4−2024 2.85 km/s 100 AU 2.54 km/s 2.22 km/s 1.90 km/s 1.59 km/s−4−2024 1.27 km/s 0.95 km/s 0.63 km/s 0.32 km/s 0.00 km/s−4−2024 -0.32 km/s -0.63 km/s -0.95 km/s -1.27 km/s -1.59 km/s−4−2024 ∆RA (") −4−2024 ∆ D e c ( " ) -1.90 km/s −4−2024-2.22 km/s −4−2024-2.54 km/s −4−2024-2.85 km/s −4−2024-3.17 km/s 0.150.380.610.841.061.291.52 J y / b e a m Fig. 1.— Images of select channels from the high-resolution CO(3-2) data (colored filled contours)along with the best fit model (black contours). Contour levels are set at 10%,25%,40%,... of thepeak flux, where 10% peak flux=0.15Jy/beam ∼ σ , as marked on the scale. Overall the modelsuccessfully reproduces much of the emission. 27 – −4−2024 2.85 km/s 100 AU 2.54 km/s 2.22 km/s 1.90 km/s 1.59 km/s−4−2024 1.27 km/s 0.95 km/s 0.63 km/s 0.32 km/s 0.00 km/s−4−2024 -0.32 km/s -0.63 km/s -0.95 km/s -1.27 km/s -1.59 km/s−4−2024 ∆RA (") −4−2024 ∆ D e c ( " ) -1.90 km/s −4−2024-2.22 km/s −4−2024-2.54 km/s −4−2024-2.85 km/s −4−2024-3.17 km/s 0.150.380.610.841.061.291.52 J y / b e a m Fig. 2.— Images of select channels from the high-resolution CO(3-2) data (colored filled con-tours) along with images generated from the difference in visibilities (red and black contours at10%,25%,40%,... of the peak flux). Red contours mark where the model is brighter than the datawhile black contours mark where the data is brighter than the model. The model does match theoverall shape and strength of the emission, while there are still some residuals. 28 – − . − . − . − . − . q R c (AU) . . . . . . . . . . v turb /c s . . . . . . . . . . T atm0 (K) . . . . . . . . . T mid0 (K) . . . . . . . incl ( o ) Fig. 3.— PDFs, normalized so that the total area under each distribution is one, for each parameterderived from fitting the low-resolution CO(3-2) data. 29 –
R (AU) Z ( A U ) CO 3-2 τ=10.01g/cm l o g n R (AU) Z ( A U ) CO 3-2 τ=1 . . . . . . . l o g n Fig. 4.— Gas density (colored contours) for the best model fit to the low-resolution CO(3-2) data.The density is only marked in the region of the disk where CO exists. In the left panel the blackcontours show the gas temperature while in the right panel the black contours mark the soundspeed in units of km s − . In both panels the dashed line shows the location of the CO(3-2) τ =1surface, which effectively traces the photodissociation boundary, while the upper/lower dotted linesmark the height where the cumulative surface density, as measured from z=170AU, is equal to 0.01and 0.1 g cm − . T atm0 (K) T m i d ( K ) Fig. 5.— In fitting the low-resolution CO(3-2) data we see evidence for a degeneracy between T mid0 and T atm0 . As T mid0 increases the disk becomes puffier, raising the τ = 1 surface higher in thedisk where it is warmer, and T atm0 must then adjust downward to maintain the same flux from the τ = 1 surface. 30 – −6 −4 −2 0 2 4 6Velocity (km/sec)0.00.20.40.60.81.01.2 A m p li t u d e Datav turb =0.1km/sv turb <0.04c s Fig. 6.— CO(3-2) high resolution spectra (black line) compared to the median model when turbu-lence is allowed to move toward very low values (red dot-dashed lines) or when it is fixed at 0.1 kms − (blue dashed lines). All spectra have been normalized to their peak flux to better highlight thechange in shape. The models with weak turbulence provide a significantly better fit to the datadespite the fact that the turbulence is smaller than the spectral resolution of the data. 31 – -3 -2 -1 v turb /c s W i d t h ( " ) −4−2024 ∆RA (") −4−2024 ∆ D e c ( " ) Fig. 7.— [Top] Width of the emission through the central velocity channel (left panel) along a linethat runs perpendicular to the emission (marked on the right panel, which shows the data in thecentral velocity channel along with a v turb =0.4c s model). As turbulence is varied in the models,the emission broadens. The vertical lines mark the spectral resolution for the low-res (right) andhigh-res (left) CO(3-2) spectra. The images change substantially with turbulence, allowing us toprobe non-thermal motion below the spectral resolution of our data. [Bottom] Model channelmaps for turbulence varying between zero (left panel), the fiducial model result (center panel) andv turb =0.4c s . The flux in these channel maps has been normalized to its peak value to highlight thechange in morphology, which is substantial for large turbulence motion. 32 – − . − . − . − . − . − . − . q R c (AU) . . . . . . . . . v turb /c s T atm0 (K) . . . . . . incl ( o ) Fig. 8.— Posterior distributions, normalized so that the total area under each distribution is one,for the five parameters used in fitting to the CO(2-1) line. −4−2024 -2.22 km/s 100 AU -1.90 km/s -1.59 km/s -1.27 km/s -0.95 km/s−4−2024 -0.63 km/s -0.32 km/s 0.00 km/s 0.32 km/s 0.63 km/s−4−2024 ∆RA (") −4−2024 ∆ D e c ( " ) J y / b e a m Fig. 9.— Channel maps for the central 15 channels of the CO(2-1) data (colored contours) withthe best fit model (black contours). Contour levels are set at 10%,25%,40%,... of the peak flux,where 10% peak flux=0.09Jy/beam ∼ σ . 33 – −4−2024 -2.22 km/s 100 AU -1.90 km/s -1.59 km/s -1.27 km/s -0.95 km/s−4−2024 -0.63 km/s -0.32 km/s 0.00 km/s 0.32 km/s 0.63 km/s−4−2024 ∆RA (") −4−2024 ∆ D e c ( " ) J y / b e a m Fig. 10.— Channel maps for the central 15 channels of the CO(2-1) data (colored filled contours)along with images generated from the difference in visibilities (red and black contours). Redcontours mark where the model is brighter than the data while black contours mark where the datais brighter than the model. The model matches much of the images. −6 −4 −2 0 2 4 6Velocity (km/sec)−2024681012 A m p li t u d e Datahigh v turb low v turb −6 −4 −2 0 2 4 6Velocity (km/sec)0.00.20.40.60.81.0 A m p li t u d e Fig. 11.— Visibility spectra for CO(2-1) comparing the data (solid black line) to a model with lowturbulence (v turb =0.037c s , blue dotted line) and a model with high turbulence (v turb =0.3c s , reddashed line). The other parameters in these models have been allowed to adjust to find the best fit.The bottom panel shows the spectra normalized to their peak flux to better highlight the differencesin line shape. While the high turbulence model nominally provides a better fit to the CO(2-1) data,the low turbulence model better matches the line profile, in particular the peak-to-trough ratio,which is sensitive to the turbulence. 34 – −4−2024 ∆RA (") −4−2024 ∆ D e c ( " ) -1.27 km/s 100 AU −4−2024-0.95 km/s −4−2024-0.63 km/s −4−2024-0.32 km/s 0.09440.18880.28320.37770.47210.56650.66090.75530.8497 J y / b e a m Fig. 12.— Imaged residuals (black solid contours) compared to the CO(2-1) data (colored filledcontours) for the low-turbulence fit. The residuals, indicating that the model underestimates thedata, are only found on the north-east side of the disk. This asymmetry is likely a radiative transfereffect related to the different heights probed on the near and far side of the disk, and the resultingtemperature structure. When the models cannot fully adjust the temperature structure to accountfor this asymmetry, they increase the turbulence leading to an overestimate of its true value. 35 – . . . . . . v turb /c s . . . . . γ . . . . . . . . . . incl ( o ) dep dep v t u r b / c s Fig. 13.— Posterior distributions, normalized so that the area under each distribtion is one, forthe four parameters used in fitting to the CO(2-1) line. We also find evidence for a degeneracybetween turbulence and the depletion factor (=the CO/CO abundance ratio); as the depletionof CO increases more of the emission arises from close to the midplane diminishing the verticalextent of the disk in the images and to match the full spatial extent of the emission the turbulencemust be increased. 36 – −4−2024 1.66 km/s 100 AU 1.33 km/s 1.00 km/s 0.66 km/s 0.33 km/s−4−2024 0.00 km/s -0.33 km/s -0.66 km/s -1.00 km/s -1.33 km/s−4−2024 ∆RA (") −4−2024 ∆ D e c ( " ) -1.66 km/s −4−2024-1.99 km/s −4−2024-2.32 km/s −4−2024-2.66 km/s −4−2024-2.99 km/s 0.050.110.180.250.320.380.45 J y / b e a m −4−2024 1.66 km/s 100 AU 1.33 km/s 1.00 km/s 0.66 km/s 0.33 km/s−4−2024 0.00 km/s -0.33 km/s -0.66 km/s -1.00 km/s -1.33 km/s−4−2024 ∆RA (") −4−2024 ∆ D e c ( " ) -1.66 km/s −4−2024-1.99 km/s −4−2024-2.32 km/s −4−2024-2.66 km/s −4−2024-2.99 km/s 0.050.110.180.250.320.380.45 J y / b e a m Fig. 14.— Channel maps for the model (top, solid black contour) and the image generated fromthe difference in visibilities (bottom,red are model > data, black are data > model) for the fit to the CO(2-1) data (filled contours). Contours are set at 10%,25%,40%,... of the pak flux, where 10%peak flux=0.045 Jy/beam ∼ σ . This model is able to successfully reproduce much of the emission. 37 – . . . . . . . v turb /c s . . . . . . . γ . . . . . . . . incl ( o ) dep dep γ Fig. 15.— Posterior distributions, normalized so that the area under each distribution is one, forthe four parameters used in fitting to the C O(2-1) line. We find a strong anticorrelation with thedepletion factor and γ ; this is not surprising since both control the surface density of gas. 38 – −4−2024 2.17 km/s 100 AU 1.83 km/s 1.50 km/s 1.17 km/s 0.83 km/s−4−2024 0.50 km/s 0.17 km/s -0.17 km/s -0.50 km/s -0.83 km/s−4−2024 ∆RA (") −4−2024 ∆ D e c ( " ) -1.17 km/s −4−2024-1.50 km/s −4−2024-1.83 km/s −4−2024-2.17 km/s −4−2024-2.50 km/s 0.03620.07240.10860.14480.18100.2172 J y / b e a m −4−2024 2.17 km/s 100 AU 1.83 km/s 1.50 km/s 1.17 km/s 0.83 km/s−4−2024 0.50 km/s 0.17 km/s -0.17 km/s -0.50 km/s -0.83 km/s−4−2024 ∆RA (") −4−2024 ∆ D e c ( " ) -1.17 km/s −4−2024-1.50 km/s −4−2024-1.83 km/s −4−2024-2.17 km/s −4−2024-2.50 km/s 0.03620.07240.10860.14480.18100.2172 J y / b e a m Fig. 16.— Channel maps for the model (top, solid black contour) and the image generated fromthe difference in visibilities (bottom,red are model > data, black are data > model) for the fit to theC O(2-1) data (filled contours). Contours are in units of 0.036Jy/beam (=5 σ ). This model is ableto successfully reproduce much of the emission. 39 – −6 −4 −2 0 2 4 6Velocity (km/sec)0.00.20.40.60.81.01.21.41.6 A m p li t u d e Datalow v turb high v turb
Fig. 17.— Visibility spectra comparing the data (black solid line) to the best fit model withturbulence fixed at a low value (v turb =0.038c s , blue dotted lines) and the best fit model with highturbulence (v turb =0.4c s , red dashed line). The difference in spectra is small compared to the largedifference in turbulence. 40 –Fig. 18.— Height of the τ = 1 surface for the best fit model derived for each line (filled contours).Overplotted are the observations for each line (similar contour levels as in Figures 1,9,14,16).Regions where there is disk emission but no marked τ =1 surface are optically thin. The differentlines probe different heights within the disk, giving us a handle on turbulence and temperature asa function of height within the disk. 41 – −6 −4 −2 0 2 4 6Velocity (km/sec)0510152025 A m p li t u d e CO(3-2) −6 −4 −2 0 2 4 6Velocity (km/sec)0246810 A m p li t u d e CO(2-1) −6 −4 −2 0 2 4 6Velocity (km/sec)012345 A m p li t u d e CO(2-1) −6 −4 −2 0 2 4 6Velocity (km/sec)0.00.20.40.60.81.01.21.41.6 A m p li t u d e C O(2-1)
Fig. 19.— Visibility spectra for all four emission lines (dark) along with the median model (red-dashed). The model fits have weak turbulence (v turb < s ) and are able to consistently fit allfour lines, although it does underestimate CO. 42 – −4−2024 ∆RA (") −4−2024 ∆ D e c ( " ) CO(3-2) −4−2024 ∆RA (") −4−2024 ∆ D e c ( " ) CO(2-1) −4−2024 ∆RA (") −4−2024 ∆ D e c ( " ) CO(2-1) −4−2024 ∆RA (") −4−2024 ∆ D e c ( " ) C O(2-1)
Fig. 20.— Emission from the central velocity channel of all four lines (colored filled contours)compared to the model defined by the median of the PDFs (black lines). The model images arewell-matched to the data, demonstrating that a model with low-turbulence is able to accurately fitall four lines. 43 –
R (AU) Z ( A U ) CO 3-2 τ=1 l o g n Fig. 21.— Gas density (colored contours) and temperature (black contours) for the best model fitto all four CO lines. The density is only marked in the region of the disk where CO exists. 44 –
REFERENCES
Aikawa, Y. & Herbst, E. 1999, A&A, 351, 233Andrews, S.M., Wilner, D.J., Hughes, A.M., Qi, C., & Dullemond, C.P. 2010, ApJ, 723, 1241Andrews, S., et al. 2012, ApJ, 744, 162Armitage, P.J. 2011, ARA&A, 49, 195Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33Bai, X.-N., & Stone, J. M. 2011, ApJ, 736, 144Balbus, S.A. & Hawley, J.F. 1991, ApJ, 376, 214Balbus, S. A., & Hawley, J. F. 1998, Reviews of Modern Physics, 70, 1Bergin, E.A., Langer, W.D., & Goldsmith, P.F. 1995, ApJ, 441, 222Blandford, R. D., & Payne, D. G. 1982, MNRAS, 199, 883Blum, J., & Wurm, G. 2000, Icarus, 143, 138Blum, J., & Wurm, G. 2008, ARA&A, 46, 21Carballido, A., Fromang, S., & Papaloizou, J. 2006, MNRAS, 373, 1633Carr, J.S., Tokunaga, A.T., & Najita, J. 2004, ApJ, 603, 213Ciesla, F.J. 2007 ApJ, 654, L159Cleeves, L. I., Adams, F. C., & Bergin, E. A. 2013, ApJ, 772, 5Crida, A., Morbidelli, A., & Masset, F. 2006, Icarus, 181, 587D’Alessio, P., Calvet, N., Hartmann, L., Franco-Hern´andez, R., & Serv´ın, H. ApJ, 638, 314Dartois, E., Dutrey, A., & Guilloteau, S. 2003, A&A, 399, 773de Gregorio-Monsalvo, I., et al. 2013, A&A,, 557, 133Dubrulle, B., Morfill, G., & Sterzik, M. 1995, Icarus, 114, 237Dzyurkevich, N., Turner, N. J., Henning, T., & Kley, W. 2013, ApJ, 765, 114Favre, C., Cleeves, L.I., Bergin, E.A., Qi, C., & Blake, G.A. 2013, ApJ, 776L, 38Flock, M., Dzyurkevich, N., Klahr, H., Turner, N. J., & Henning, T. 2011, ApJ, 735, 122Foreman-Mackey, D., Hogg, D.W., Lang, D., & Goodman, J. 2013, PASP, 125, 306 45 –Forgan, D., Armitage, P. J., & Simon, J. B. 2012, MNRAS, 426, 2419France, K., Herczeg, G.J., McJunkin, M., & Pention, S.V. 2014, ApJ, 794, 160Fromang, S. & Nelson, R.P. 2006, A&A, 457, 343Fromang, S., & Nelson, R. P. 2009, A&A, 496, 597Fromang, S., Latter, H., Lesur, G., & Ogilvie, G.I. A&A, 552, 71Fung, J., Shi, J.-M., & Chiang, E. 2014, ApJ, 782, 88Furuya, K., & Aikawa, Y. 2014, ApJ, 790, 97Gammie, C.F. 1996, ApJ, 457, 355Goodman, J., & Weare, J. 2010, Commun. Appl. Math. Comput. Sci., 5, 65Gressel, O., Nelson, R. P., & Turner, N. J. 2011, MNRAS, 415, 3291Guilloteau, S., Dutrey, A., Wakelam, V., Hersant, F., Semenov, D., Chapillon, E., Henning, T., &Pi´etu, V. 2012, A&A, 548, 70Gundlach, B., Kilias, S., Beitz, E., & Blum, J. 2011, Icarus, 214, 717Hartmann, L., Calvet, N., Gullbring, E., & D’Alessio, P. 1998, ApJ, 495, 385Hartmann, L., Hinkle, K., & Calvet, N. 2004, ApJ, 609, 906Hughes, A.M., Wilner, D.J., Qi, C., & Hogerheijde, M.R. 2008, ApJ, 678, 1119Hughes, A.M., Wilner, D.J., Andrews, S.M., Qi, C., & Hogerheijde, M.R. 2011, ApJ, 727, 85Ida, S., Guillot, T., & Morbidelli, A. 2008, 686, 1292Isella, A., Testi, L., Natta, A., Neri, R., Wilner, D., & Qi, C. 2007, A&A, 469, 213Johansen, A., & Klahr, H. 2005, ApJ, 634, 1353Johansen, A., Oishi, J. S., Mac Low, M.-M., et al. 2007, Nature, 448, 1022Jonkheid, B., Dullemond, C.P., Hogerheijde, M.R., & van Dishoeck, E.F. 2007, A&A, 463, 203Juh´asz, A., et al. 2010, ApJ, 721, 431Klaassen, P.D., et al. 2013, A&A, 555, 73Klahr, H. H., & Henning, T. 1997, Icarus, 128, 213Klahr, H. H., & Bodenheimer, P. 2003, ApJ, 582, 869 46 –Kley, W., & Nelson, R. P. 2012, ARA&A, 50, 211Kolmogorov, A. 1941, Akademiia Nauk SSSR Doklady, 30, 301Horne, K., & Marsh, T.R. 1986, MNRAS, 218, 761Laughlin, G., Steinacker, A., & Adams, F.C. 2004, 608, 489Lynden-Bell, D. & Pringle, J.E. 1974, MNRAS, 168, 603Lyra, W., Turner, N., & McNally, C. 2014, arXiv:1410.8092Marcus, P., Pei, S., Jiang, C.-H., et al. 2014, arXiv:1410.8143Mathews, G.S., et al. 2013, A&A, 557, 132Meeus, G., Waters, L. B. F. M., Bouwman, J., et al. 2001, A&A, 365, 476Mendigutia, I., et al. 2013, ApJ, 776, 44Miller, K. A., & Stone, J. M. 2000, ApJ, 534, 398Miotello, A., Bruderer, S., & van Dishoeck, E. F. 2014, A&A, 572, A96Montesinos, B., Eiroa, C., Mora, A., & Merin, B. 2009, A&A, 495, 901Mori, S. & Okuzumi, S. 2015, arXiv:1505.04896Najita, J.R., Carr, J.S., Glassgold, A.E., Shu, F.H., & Tokunaga, A.T. 1996, ApJ, 462, 919Najita, J.R., Doppmann, G.W., Carr, J.S., Graham, J.R. & Eisner, J.A. 2009, ApJ, 691, 738Nelson, R. P., Gressel, O., & Umurhan, O. M. 2013, MNRAS, 435, 2610Pavlyuchenkov, Y., et al. 2007, ApJ, 669, 1262Perez-Becker, D., & Chiang, E. 2011, ApJ, 735, 8Poppe, T., Blum, J., & Henning, T. 2000, 533, 454Qi, C., et al. 2011, ApJ, 740, 84Qi, C., et al. 2013, Science, 341, 630Qi, C., Oberg, K., Wilner, D.J., & Rosenfeld, K.A. 2013, ApJ, 765, 14Reg´aly, Z., S´andor, Z., Dullemond, C. P., & Kiss, L. L. 2011, A&A, 528, A93Reboussin, L., Wakelam, V., Guilloteau, S., Hersant, F., & Dutrey, A. 2015, arXiv:1505.01309Rosenfeld, K. A., Qi, C., Andrews, S. M., et al. 2012, ApJ, 757, 129 47 –Rosenfeld, K.A., Andrews, S.M., Hughes, A.M., Wilner, D.J., & Qi, C. 2013, ApJ, 774, 16Rosenfeld, K. A., Andrews, S. M., Wilner, D. J., Kastner, J. H., & McClure, M. K. 2013, ApJ, 775,136Rosenfeld, K. A., Chiang, E., & Andrews, S. M. 2014, ApJ, 782, 62Sano, T., Miyama, S. M., Umebayashi, T., & Nakano, T. 2000, ApJ, 543, 486Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337Shi, J. & Chiang, E. 2014, 78, 34Simon, J.B., Bai, X., Armitage, P.J., Stone, J.M., & Beckwith, K. 2013, ApJ, 775, 73Simon, J.B., Hughes, A.M., Flaherty, K.M., Bai, X.-N., & Armitage, P.J. in prepSemenov, D., & Wiebe, D. 2011, ApJS, 196, 25Tilling, I., et al. 2012, A&A, 538, 20Turner, N. J., & Drake, J. F. 2009, ApJ, 703, 2152Turner, N.J., et al. 2014, PPVIVisser, R., van Dishoeck, E.F., & Black, J.H. 2009, A&A, 503, 323Walsh,C., Millar, T.J., & Nomura, H. 2010, ApJ, 722, 1607Walsh, C., Juh´asz, A., Pinilla, P., et al. 2014, ApJ, 791, L6Williams, J.P., & Best, W.M.J. 2014, ApJ, 788, 59Wilson, T. L. 1999, Reports on Progress in Physics, 62, 143Youdin, A. N., & Shu, F. H. 2002, ApJ, 580, 494