EDP Sciences
Free Access
Issue
A&A
Volume 605, September 2017
Article Number A16
Number of page(s) 25
Section Interstellar and circumstellar matter
DOI https://doi.org/10.1051/0004-6361/201630329
Published online 01 September 2017

© ESO, 2017

1. Introduction

The advent of the Acatama Large Millimeter/submillimeter Array (ALMA) is providing an unprecedented level of angular resolution and sensitivity at (sub)mm wavelengths. This technological step forward allows us to determine more precisely fundamental quantities of protoplanetary disks, the cradles of planet formation. In particular, two quantities of great interest are the outer radius of disks and the surface density radial profile Σ(R) (see Williams & Cieza 2011; Armitage 2011, 2015; Morbidelli & Raymond 2016, for recent reviews). Both of them are of fundamental importance because they are directly related to the amount of mass present in the disk, which is obviously linked to the planet formation potential of a system. The surface density profile also provides one of the main parameters for any planet formation model. Moreover, in the commonly assumed disk evolution framework (the so-called α-disk scenario, Lüst 1952; Shakura & Sunyaev 1973; Lynden-Bell & Pringle 1974), the gas outer radius determines the global evolution of the gaseous component, since the viscous timescale is directly related to the radial extent of the disk. In star forming regions, it can also imprint the respective importance of environmental effects, such as star-disk encounters (e.g. Clarke & Pringle 1993; Olczak et al. 2006; Rosotti et al. 2014; Dai et al. 2015; Vincke et al. 2015) and external photoevaporation (e.g. Hollenbach et al. 1994; Johnstone et al. 1998; Adams et al. 2004; Facchini et al. 2016; Haworth et al. 2016a), which can both truncate the disks radially.

Interestingly, observations have shown that the outer radius of the gaseous component of a disk, as probed by the bright 12CO emission, is generally larger than that of the dust component, probed by (sub)mm continuum emission. Historically, this difference has been considered to be due to optical depth effects, with the gas lines (in particular the 12CO line) more optically thick than the continuum emission at similar wavelengths (e.g. Dutrey et al. 1998; Guilloteau & Dutrey 1998; Panić et al. 2009). This assumption motivated observers to try fitting both the gas tracers and continuum emission with the same surface density profiles, leading Hughes et al. (2008) and Andrews et al. (2009) to propose a tapered power law profile, theoretically justified by the self-similar solutions by Lynden-Bell & Pringle (1974). However, recent observations with higher sensitivity and angular resolution have shown that the apparent friction between the dust (continuum) and gas (molecular) intensity profiles in disks cannot be reconciled even with such a surface density profile. In well resolved systems, in particular those imaged with ALMA, the dust outer edge decreases too sharply with radius, and cannot be reproduced by a tapered outer disk (e.g. Andrews et al. 2012, 2016; de Gregorio-Monsalvo et al. 2013; Piétu et al. 2014; Cleeves et al. 2016).

The natural explanation that has been proposed to interpret such observations is a combination of grain growth and consequent radial drift of large particles from the outer to the inner disk. The fact that dust particles can grow to at least mm/cm sizes in protoplanetary disks has now been proven by many observations at different radio wavelengths (e.g. Testi et al. 2003; Natta et al. 2004; Lommen et al. 2007, 2010; Andrews & Williams 2005, 2007; Ricci et al. 2010a,b). Recent reviews on the topic are Testi et al. (2014), Andrews (2015), Birnstiel et al. (2016). Moreover, there is increasing evidence that the maximum grain size attained in a disk is a function of distance from the star, as expected from grain growth models (e.g. Guilloteau et al. 2011; Pérez et al. 2012, 2015; Menu et al. 2014; Tazzari et al. 2016). These observational constraints have been accompanied by further theoretical modelling of both grain growth and radial drift processes (e.g. Birnstiel et al. 2010, 2012). In particular, Birnstiel & Andrews (2014) have shown that the sharp edge observed in the (sub-)mm continuum is predicted by dust evolution models.

In order to investigate this situation further, simultaneous modelling of both dust and gas is needed. The size distribution evolution of the dust particles does affect the chemistry occurring in protoplanetary disks, in particular the CO abundance and excitation, via multiple physical-chemical processes some of which are directly linked to the total dust surface area available at any point of the disk. A first important effect is that grain growth, vertical settling, and radial drift affect the penetration depth of the UV (ultraviolet) photons into the disk; these photons can originate both from the central star or from the surrounding environment. Some studies have started addressing this topic in static disks, for example Jonkheid et al. (2004, 2007), Nomura et al. (2007), Vasyunin et al. (2011), Fogel et al. (2011), Woitke et al. (2016), Cleeves (2016). Enhanced UV penetration affects both the chemistry and the gas temperature, which both determine the molecular (and atomic) emission. Moreover, the relation between UV penetration depth and grain growth can also have an important dynamical effect, making the disk more or less prone to photoevaporate (Gorti et al. 2015; Facchini et al. 2016). Another important connection between dust evolution and chemistry is that the amount of dust surface area available determines the level of thermal coupling between the gas and solid phases. Finally, the dust surface area affects the balance between freeze-out and desorption, since they both depend on the total surface area. The non linear combination of all these effects make the modelling of CO emission lines (and of other molecules) less straightforward than simply assumed.

In this paper, we aim to explore how the properties of physically justified dust grain size distributions affect the thermal and chemical structure of a protoplanetary disk. In particular, we focus on how dust evolution (grain growth, radial drift, and vertical settling) affects the observability of both dust and gas outer radius and surface densities, as probed by (sub-)mm continuum observations and low-J rotational lines of different CO isotopologues. To do so, we combine the dust evolution models by Birnstiel et al. (2015) with the thermo-chemical code DALI (Dust And LInes, Bruderer et al. 2012; Bruderer 2013). The method used in the paper is detailed in Sect. 2 and the setup parameters for the models investigated are given in Sect. 3. In Sect. 4 we describe our findings, which are then discussed in Sect. 6. In Sect. 7 we summarise the results and draw our conclusions.

2. Method

thumbnail Fig. 1

Diagram of the method used in this paper. The blue box indicates the input parameters and the green box indicates the output observables. Boxes with orange borders show new modules, compared to the standard DALI code. Orange text underlines modifications in previous DALI modules. The red text lists input parameters needed to compute the dust distribution.

Open with DEXTER

In this section we present the method used in this paper to couple dust evolution with the thermo-chemical processes in a protoplanetary disk. The method consists in merging the thermo-chemical code DALI (Bruderer et al. 2009, 2012; Bruderer 2013), which has been widely used to model gas emission in protoplanetary disks (e.g. Miotello et al. 2014, 2016; Bruderer et al. 2014, 2015; van der Marel et al. 2015, 2016; Kama et al. 2016; Fedele et al. 2016), with the grain size distribution reconstruction routine by Birnstiel et al. (2015). More precisely, given an initial gas density distribution, the radial dependence of the grain size distribution is computed with the grain size reconstruction routine. Subsequently, the vertical distribution of the dust is calculated solving the advection-diffusion equation of vertical settling in steady-state at every radius for every dust bin considered. The local grain size distribution is then used to compute the opacities at every grid point of the disk (radius and height), which are then used in the continuum radiative transfer module. Moreover, the local particle distribution is used to compute the total dust surface area available for the thermo-chemical processes, in particular gas-grain collisions, H2 formation rate, freeze-out, thermal and non-thermal desorption, and hydrogenation. With these new rates, both the gas temperature and the chemical abundances are computed self-consistently. A schematic of the method is illustrated in Fig. 1.

2.1. Gas density distribution and stellar parameters

The main input physical quantity for the used method is the gas density structure. We use the following gas surface density Σgas dependence on cylindrical radius R (Andrews et al. 2009; Bruderer 2013): (1)where Rc is the characteristic radius determining the radial length scale of the disk, and Σc determines the total mass of the disk. Assuming hydrostatic equilibrium, isothermality in the vertical direction, and zR (where z is the vertical coordinate), we obtain the gas mass density ρgas: (2)where h = hc(R/Rc)ψ, ψ is the flaring index, and hc is the aspect ratio of the disk at the characteristic radius. The expression in Eq. (1) is justified by the self-similar solutions by Lynden-Bell & Pringle (1974). This surface density profile would require the system to be older than a viscous timescale, which is not always the case. Thus the outer regions of protoplanetary disks might have surface densities that do not follow such profiles and instead reflect the initial conditions. However, we choose to use it throughout this paper since it can describe the surface densities with very few parameters, and to compare our models with those present in the literature exploiting the same parametrisation.

The other input parameters describe the stellar properties, in particular the stellar mass M and radiation parameters. The stellar spectrum is modelled as a black body, determined by its total luminosity L and effective temperature Teff. Accretion onto the star is also considered, which can be important in setting the total high energy radiation. The accretion luminosity is modelled as a black body emitting at 10 000 K, with total luminosity equalling the total potential energy per time unit emitted on the stellar surface by a mass accretion rate (Kama et al. 2016). This prescription is used to estimate the UV excess emitted in the accretion shock. We note, however, that the disk evolution is not modelled self-consistently with this mass accretion rate.

2.2. Radial grain size distribution

In order to determine the grain size distribution at every radius, given the gas density structure defined as above, the semi-analytical prescription described in Birnstiel et al. (2015) is used, which has proven to be a good representation of the more complete numerical models by Birnstiel et al. (2010). Here we briefly summarise the main physical mechanisms determining the radial distribution of particles sizes. The analytical details are given in Birnstiel et al. (2015).

Grain growth in disks occurs due to the mutual collisions of dust particles. Depending on the relative velocity Δv, such collisions lead either to sticking of the particles or to fragmentation or erosion (or various other outcomes that are not modelled here). The boundary between these two regimes is set by the fragmentation velocity vfrag, where we assume that perfect sticking occurs when Δv<vfrag. The relative motions of dust particles and their transport within the disks are driven by four processes: turbulent mixing, Brownian motions, gas drag, and radial drift (Weidenschilling 1977). The combination of such processes can lead to two different regimes: 1) dust grains grow until the relative velocities are high enough that fragmentation occurs (fragmentation barrier); 2) the drifting timescale becomes shorter than the collisional timescale, and the maximum particle size is determined by radial drift (drift barrier).

The physical conditions of protoplanetary disks are such that the fragmentation barrier is unlikely to dominate in the entire disk. It is possible to assign a fragmentation radius Rfrag, inside which the maximum grain size is set by fragmentation. In the outer regions of the disk, the maximum grain size is determined by radial drift. The actual boundary between the two regimes is smoothed by diffusion. As in Birnstiel et al. (2015), in this paper the grain size distribution inside Rfrag is determined by the analytical fits of coagulation and fragmentation equilibrium by Birnstiel et al. (2011). Outside the fragmentation radius, semi-analytical treatments including diffusion and radial drift are exploited (Pavlyuchenkov & Dullemond 2007; Birnstiel et al. 2015). In this paper, a pd parameter equalling 2.5 is assumed, where pd is the power law index of the radial dependence of a given grain size bin in the outer disk (see Sect. 2.3 in Birnstiel et al. 2015, for details).

The final radial distribution of grain sizes will depend on many input parameters, in particular the gas surface density of the disk, the dust-to-gas mass ratio Δdg, the mass of the star, and the dust temperature profile, which is implicitly set by the flaring angle ψ. Besides these quantities, which are all implemented as initial input parameters (see Sect. 2.1), the output of the grain size reconstruction routine depends on the fragmentation velocity, the disk turbulence, and the grain mass density ρgr. Laboratory experiments have shown that fragmentation velocity can assume values of a few 1−10 m s-1 (e.g. Blum & Wurm 2000; Gundlach & Blum 2015). Throughout this paper we assume vfrag = 10 m s-1. We parametrise turbulence via the dimensionless parameter α (Shakura & Sunyaev 1973), and we will vary this parameter in the simulations presented below. Finally, we assume compact dust grains, with ρgr = 2.5 g cm-3, and assume the dust-to-gas ratio to be uniform in the whole disk, with Σdust(R) / Σgas(R) = Δdg = 0.01. In future work, we will address how more complex models, with a radially varying dust-to-gas ratio, will affect the final results. A combination of some of these parameters can easily lead to an analytical estimate of the fragmentation radius (Birnstiel et al. 2015): (3)Given the input parameters defined in Sect. 2.1 and the additional ones listed in this section, a radially dependent grain size distribution Σdust(R,a) can be retrieved, which is defined as the dust surface density at a given radius R, with particle sizes between a and a + da (in radius). For numerical reasons, a discrete grain size grid of 250 size bins logarithmically spaced is used, with a minimum grain size amin = 50 Å, and a maximum grain size amax = 1 m.

thumbnail Fig. 2

Gas density structure used in all models.

Open with DEXTER

2.3. Dust vertical settling

Once the radial dependence of the grain size distribution is computed, the next step is to determine the vertical distribution of dust particles. This can be done by taking into account the physical processes responsible for dust settling, with the assumption that the turbulence acting in the vertical direction is the same as that responsible for the angular momentum transport and for the turbulent motions leading to particle collisions, thus with a gas diffusion coefficient: (4)where ν is the kinematic viscosity (Shakura & Sunyaev 1973), and .

Following Birnstiel et al. (2010), we use the results derived by Youdin & Lithwick (2007) and we relate the dust diffusivity to the gas diffusivity through the Schmidt number: (5)where the Stokes number St in the Epstein regime in the disk mid-plane can be computed as (e.g. Birnstiel et al. 2011): (6)The equation regulating the vertical motions of dust particles can be written as (Dubrulle et al. 1995; Dullemond & Dominik 2004; Fromang & Papaloizou 2006; Armitage 2010): (7)where the friction timescale of a particle of size a in the Epstein regime is: (8)and the thermal velocity of the gas vth is: (9)The temperature T (or better, the sound speed cs) used in Eq. (9) is derived from the flaring angle of the disk, which is an input parameter of the models. In particular, cs = hR/ Ω (from the ideal gas law).

In this work, we assume that the dust structure has reached a steady-state in the vertical direction, that is when the gravitational force acting along the vertical direction and the aerodynamical drag caused by turbulent motions balance out (Fromang & Nelson 2009). We can therefore neglect the term on the left hand side of Eq. (7), and solve the advection-diffusion equation for ρd(R,z,a) with a 1D integration in z. We normalise each bin by requiring that: (10)We note that the “canonical” grain size distribution f(R,z,a) is related to ρdust(R,z,a) by: (11)where f(R,z,a)da is the number of grains per unit volume, with size between a and a + da. At this stage, we have obtained the dust density distribution for every size bin in every grid point of the disk.

A limitation of our model is the underlying vertical gas density distribution (see Eq. (2)), which implicitly assumes isothermality along the z-direction. As shown in this paper, and in all the literature computing the thermal balance of the gas phase material, the gas temperature is an increasing function of z, with the upper layers being more heavily irradiated by the central star. Thus, the disk upper layers are expected to “puff up” and be less dense than what we assume here, in order to reach hydrostatic equilibrium in the vertical direction. Moreover, smaller particles are expected to be stirred up more easily than what we assume, since Eq. (9) considers the temperature derived from the flaring index (i.e. assuming vertical isothermality again), and does not take into account a vertical thermal gradient. In order to obtain a completely consistent density profile, one should iterate between steps 1 and 6 of the diagram shown in Fig. 1 to obtain a vertical hydrostatic solution for both the gas and the dust density structure. At the moment, this is computationally too expensive for these complex models. Moreover, such iterations may result in disks that are too strongly flared and/or unstable in the very inner regions (e.g. Woitke et al. 2009).

2.4. Average grain size

The grain size distribution affects both the opacities and the thermo-chemistry of the disk. To reduce the computational cost of such calculations, it is possible to define average quantities of the grain size distribution, where an ensemble of dust defined by these average properties has the same mass and total surface area as the original ensemble. In particular, following Vasyunin et al. (2011), the second and third moment of the grain size distribution can be defined as: The average grain size is then given by: (14)This approach employs an average size that yields the same mass and total surface area of the original ensemble (Vasyunin et al. 2011).

Dividing the total dust mass per cm-3ρdust(R,z) by the mean mass of one dust particle we can obtain the number of grains per volume ngr. Defining the mean geometrical cross section as , the total grain geometrical cross section per volume at an (R,z) location in the disk is given by: (15)where (16)

2.5. Dust opacities and PAHs abundance

The grain size distribution is used to compute the dust opacities at every position in the disk. We have first computed a library of dust opacities containing the opacities κν(ai) for every bin size ai. These opacities are computed using a standard ISM (interstellar medium) dust composition following Weingartner & Draine (2001), with an MNR (Mathis et al. 1977) grain size distribution between ai and ai + Δai. The mass extinction coefficients are calculated using Mie theory with the miex code (Wolf & Voshchinnikov 2004) and optical constants by Draine (2003) for graphite and Weingartner & Draine (2001) for silicates. This opacity library is then used to obtain the dust opacities in every grid point of the disk by mass averaging: (17)These opacities are implemented in the Monte-Carlo continuum radiative transfer module of DALI, as described in the Appendix of Bruderer et al. (2012), Bruderer (2013). A first stage of the Monte-Carlo method computes the dust temperature Tdust, whereas a second stage computes the mean intensities over the entire spectrum, from UV to radio frequencies.

In the radiative transfer, we also consider the opacity of PAHs (Polycyclic Aromatic Hydrocarbons). PAHs are also taken into account in the thermo-chemical modelling, where they can be important heating sources via photoelectric effect. In this work, the PAHs abundance is assumed to be 0.1 of the typical ISM abundance in the whole disk (see Bruderer et al. 2012), following observations suggesting a PAHs deficit in protoplanetary disks (e.g. Geers et al. 2006; Oliveira et al. 2010).

2.6. Thermo-chemistry

The thermo-chemical module of DALI is used to determine the chemical abundances, molecular and atomic excitations, and gas temperature Tgas with non-LTE (local thermodynamical equilibrium) calculations. The details are described in Bruderer et al. (2012), Bruderer (2013). Since some reaction rates in the thermal-chemical module do depend on the total dust surface area available, we apply small changes that take into account such dependence. The details of such rates, and their dependence on the total dust surface area, are explained below. After this final step, the synthetic emission maps of both continuum, and of specific molecular lines, can be obtained by using the ray-tracer described in Bruderer et al. (2012), Bruderer (2013).

2.6.1. Gas-grain collisions

Gas-grain collisions are important in setting the thermal coupling between dust and gas. We follow the prescription by Young et al. (2004, see their Appendix), where the dust-gas energy transfer is expressed as: (18)where Sd is the dust geometrical cross section per baryon, that is:(19)The direct proportionality to the local dust-to-gas ratio δdg = ρdust/ρgas shows that the total dust cross section is higher for a larger amount of dust, whereas the inverse dependence on indicates that a smaller average grain size increases the surface-to-mass ratio of the grain size distribution. Thus, a more top-heavy grain size distribution leads to a less effective thermal coupling between the gas and dust components. We note that usually DALI assumes Sd = δdg × 6.09 × 10-20 cm2 (Young et al. 2004).

2.6.2. H2 formation rate

The H2 formation rate on the surface of dust grains depends on the total dust surface area available. Following the prescription by Cazaux & Tielens (2002b,a), we use a recombination rate for H2 equal to (see also Hollenbach et al. 1971): (20)where nH is the abundance of H atoms in the gas phase, ϵH2 is the H2 recombination efficiency, and SH(Tdust) is the sticking coefficient of H atoms on the grain surface, which depends on the dust temperature.

2.6.3. Freeze-out, desorption, and hydrogenation

The other chemical mechanisms that depend on the dust surface area in the chemical network used in this paper are freeze-out (or adsorption), evaporation (or desorption), and grain-surface hydrogenation reactions. For freeze-out, we follow Charnley et al. (2001) and Visser et al. (2009), where the adsorption rate coefficient for a molecule X can be expressed as: (21)where M(X) is the molecular mass of species X.

Thermal desorption is treated similarly. Following Visser et al. (2011), the thermal evaporation rate is prescribed as: (22)where Eb(X) is the binding energy of species X. The most relevant molecule in this work is CO. The binding energy for pure CO ice is ~855 K (Sandford & Allamandola 1993; Collings et al. 2003; Öberg et al. 2005; Bisschop et al. 2006), but it can be as high as ~1100−1500 K on H2O ice or mixed with CO2 ice (Martín-Doménech et al. 2014; Fayolle et al. 2016). We assume a binding energy for pure CO ice, 855 K, unless stated otherwise. The values of the pre-exponential factor A(X) can be found in Visser et al. (2011). The factor f(X) sets the boundary between zeroth to first order desorption, where only one monolayer of ice is considered active in the evaporation process: (23)where the number of binding sites per unit grain surface Nss = 8 × 1014 cm-2 (e.g. Hollenbach et al. 2009) and nice is the sum of number densities of all individual ice species.

Photodesorption is modelled following again Visser et al. (2011, see their Eq. (6)): (24)where FUV is the local value of UV flux (obtained from the radiative transfer module, see Sect. 2.5) and Y(X) is the number of photodesorbed molecules per grain per incident UV photon. The values used in this paper are the same as in Visser et al. (2011), and are taken from laboratory experiments (Öberg et al. 2007, 2009a,b). More recent works on the subject can be found in Paardekooper et al. (2016, and references therein). To mimic photodesorption induced by cosmic rays in the dense optically thick regions of the disk, we assume a floor value for the local UV flux of 104cm-2 s-1 (Shen et al. 2004). We do not include cosmic ray spot heating.

Finally, DALI includes a set of simplified grain-surface hydrogenation reactions, in particular to turn C, N, and O into CH4, NH3, and H2O. The simple prescription implemented is as follows: (25)where again the reactions are allowed only in the top monolayer and nhydro represents the sum of the number densities of all ice species that can be hydrogenated; in this paper, C, CH, CH2, CH3, N, NH, NH2, O, and OH. In freeze-out, evaporation, and hydrogenation rates, DALI usually assumes a typical grain size of 0.1 μm and a grain number density ngr = 10-12nH. Instead, in this paper ngr and are computed directly from the local grain size distribution.

With these new rates, the molecular abundances and excitations and gas temperature are computed. Using the ray-tracing module, the synthetic emission maps in both continuum and molecular and atomic lines can be obtained.

2.7. Standard DALI models

Since the implementation of a more realistic grain size distribution into DALI requires a few changes, we will compare the results of this new methodology with similar results from the standard DALI code to better estimate the effects of grain growth and radial drift onto the chemistry and emission of a few molecular lines. To mimic grain growth and settling, DALI usually considers two populations of dust grains, with sizes ranging between 50 Å and 1 μm for small grains and ranging between 50 Å and 1 mm for large grains (as in D’Alessio et al. 2006). The scale height for the small population is the same as for the gas, whereas the scale height for the large population is reduced to χh, where χ< 1. The mass ratio between the large and small populations is defined as f.

To be more specific, these standard parametric models have a gas density distribution as defined in Eq. (2), but a dust density distribution of:

We will compare these standard models with the more complete ones including dust evolution.

3. Models setup

As our representative model, a model using a gas density distribution that resembles early ALMA observations of the HD 163296 system is considered. We stress here that this paper is not aimed to model specifically the HD 163296 system. However, we consider it as a good representative case, since de Gregorio-Monsalvo et al. (2013) have shown that the dust and gas emission cannot be modelled simultaneously by a tapered outer disk surface density profile. Moreover, the disk is bright and large enough that resolved intensity profiles have been obtained for both continuum and CO isotopologues. We thus take disk parameters that have been used by de Gregorio-Monsalvo et al. (2013) to fit the Science Verification ALMA data. Other parametrisations for the disk structure have been used (e.g. Rosenfeld et al. 2013; Williams & McPartland 2016), but we will discuss them more thoroughly in Sect. 5. In particular, we use Rc = 125 AU, γ = 0.9, Σc = 7.03 g cm-2, hc = 0.11, and ψ = 0.25. These parameters yield a disk gas mass Mdisk ~ 7 × 10-2M. The gas density structure used in all models is shown in Fig. 2. The total dust-to-gas mass ratio Δdg is set to 0.01. We use a stellar mass of 2.47 M, with a bolometric luminosity L = 38 L, and an effective temperature of 104K (Tilling et al. 2012). The X-ray luminosity is set to LX = 6 × 1029erg s-1 (Günther & Schmitt 2009). Finally, we use an accretion rate of 4.5 × 10-7M yr-1 (Mendigutía et al. 2013).

All disk models are represented with a 2D (R,z) grid, with 150 grid points in the radial direction and 80 points in the vertical direction. The first 50 radial points are logarithmically sampled between Rin and 25 AU. We set an inner disk radius Rin = 0.42 AU, since the dust sublimation radius is ~, if a dust sublimation temperature of ~1500 K (Dullemond et al. 2001) is assumed. The grid is then sampled linearly between 25 and 800 AU. The vertical grid is linear and samples the disk between 0 and 6 local scale heights.

As standard reference models using standard DALI, we define two models, which we label STN (which stands for “standard”) and STN-SM, where the main difference between the two is that the second one does not have grains larger than 1 μm in the disk. For model STN, we consider f = 0.85 and a settling parameter χ = 0.2 (see Sect. 2.7). Instead, for model STN-SM, we set f = 0, that is, all the dust is in the population of small grains. In these reference models, we compute the dust opacities as described in Sect. 2.5, with just two large size bins. The rather unrealistic STN-SM model is used mostly to determine the thermal structure of a hypothetical disk that did not witness grain growth yet, in order to compare chemical timescales with grain growth timescales in Sect. 6.3.

For the models including dust evolution, the three new parameters defining the dust distribution are α, vfrag, and ρgr. As noted above, vfrag = 10 m s-1 and ρgr = 2.5 g cm-3. Here we focus on exploring the dependence of the dust and gas properties on turbulence, by assigning α the values of 10-2, 10-3 and 10-4, that is, ranging between high and low values of disk viscosity (e.g. Hartmann et al. 1998; Turner et al. 2014, where the latter is a recent theoretical review). We note that the first tentative measurement of disk turbulence by Hughes et al. (2011), Flaherty et al. (2015) and Teague et al. (2016) seem to indicate that turbulence is quite low (α ≲ 10-3), at least in the outer disks.

In all the models we assume other parameters that are relevant for the calculations performed for this paper. In particular, in the radiative transfer module in the first stage we use 3 × 107 photon packages for both the star and the environment radiation, which is assumed to be 1 G0, where G0 ~ 2.7 × 10-3erg s-1 cm-2 is the UV interstellar radiation field between 911 Å and 2067 Å (Draine 1978). The number of photon packages used in stage 2 is 3 × 106 in every wavelength bin (see Bruderer 2013, for more details). The number of photon packages in both stages has been chosen such that the dust temperature and the intensity field are smooth functions in the spatial grid specified above.

In the thermo-chemical module, we assume initial ISM-like abundances, with [C] / [H] = 1.35 × 10-4 and [O] / [H] = 2.88 × 10-4, where notation [X] indicates element X in all its volatile forms (i.e. not locked up in refractory dust). We do not consider CO isotope selective photodissociation (e.g. van Dishoeck & Black 1988), since the mass of the disk and the stellar mass used in these models are high enough that this effect should not be too significant (Miotello et al. 2014, 2016). The assumed cosmic ray ionisation rate is ζCR = 5 × 10-17s-1. The calculations are performed in time-dependent mode, for a timescale of 1 Myr. The upper CO emitting layers of disks reach chemical equilibrium in <1 Myr. The same timescale is used in the dust evolution calculations, where, however, the grain size distribution reaches a steady-state in ~2 × 105yr.

Finally, the ray tracing assumes a distance that is comparable to that of the HD 163296 system, 122 pc (van den Ancker et al. 1998). The synthetic observables are ray traced assuming a disk inclination of 45°, with a beam convolution of 0.52″ × 0.38″, where the position angle of the beam is 82° (de Gregorio-Monsalvo et al. 2013). The position angle of the major axis of the disk is taken as 137°.

thumbnail Fig. 3

Grain size distribution for α = 10-2, 10-3 and 10-4, from left to right, respectively. The two rows show the same plots in two different radial scales. The dashed horizontal line is the total dust-to-gas ratio Δdg = 0.01, whereas the dashed-dotted vertical line indicates the fragmentation radius Rfrag.

Open with DEXTER

thumbnail Fig. 4

From top to bottom: STN, α = 10-2, 10-3 and 10-4 models. From left to right: local dust-to-gas ratio δdg, average grain size and total dust surface area per volume ngrσ used in the thermo-chemical module. Contour levels are indicated in the colour bar.

Open with DEXTER

4. Results

In this section we show the results of the models, focusing in particular on the dust properties and on the effects that these have on the chemical abundances and excitation of the main CO isotopologues, 12CO, 13CO, and C18O.

4.1. Dust density structure and average grain size

The three values of turbulence induce important differences in the radial dependence of the grain size distribution (see Fig. 3). In the inner regions, lower viscosities lead to larger grain sizes since turbulent velocity decreases with α. Moreover, for low viscosities, most of the grain size distribution is limited by radial drift and not by fragmentation. This is apparent by looking at the location of the fragmentation radius Rfrag for the three different cases, where Rfrag varies between ~10500 AU with α ranging between 10-4 and 10-2 (in fact Rfrag is expected to be roughly linear with turbulence parameter, see Eq. (3)). In the outer regions, where the gas surface density exponentially decays, higher viscosities lead to grain size distributions which are less bottom-heavy. The reason is that in the outer exponential tail of the surface density profile, even the smallest grains have a Stokes number 1, and thus all grains reside in the smallest size bin, being the maximum grain size set by the radial drift limit. However, as viscosity gets higher, some larger grains from smaller radii are diffused outwards, leading to a more top-heavy grain size distribution in the outer regions of the disk when compared to the one resulting from lower viscosities. As we will show, this is important in setting the penetration depth of UV photons into the outer disk. All three models show a quite sharp transition between mm-size particles and micron-size grains at ~200 AU. Such a transition becomes sharper for lower viscosities, since the maximum grain size attained at every radius becomes a steeper function of the distance from the star. Interestingly, in all models there is a significant depletion of small grains just outside the fragmentation radius. As explained in Birnstiel et al. (2015), this occurs because the relative motions between dust particles are not high enough to induce fragmentation. Thus, particles grow to large enough sizes that they radially drift inwards and the small particles are not replenished efficiently. This is in fact a possible origin of gaps in scattered light images, as detailed in Birnstiel et al. (2015). Such a gap in scattered light should be correlated to a possibly small enhancement of the opacity at (sub-)mm wavelengths.

The results of the combination of grain growth, radial drift, and vertical settling in the dust density distribution are apparent in the left panels of Fig. 4, where the grain size evolution models are compared to the STN one. First, the dust density structure in the most viscous case is similar to that of the STN model. As turbulence gets lower, the vertical settling becomes more severe, with differences in the dust density in the upper layers of the disk of more than four orders of magnitude between the α = 10-2 and 10-4 cases. This is due to a combination of grain growth and vertical settling: for low viscosities, larger size bins are more populated because the fragmentation limit is very close to the central star, and the low turbulence is very inefficient in maintaining the massive dust grains at significant altitudes in the disk.

The average grain size and the total dust surface area per volume ngrσ are two other important quantities that are directly related to the dust density structure (see Sect. 2.4). From the combination of grain growth and vertical settling, lower viscosities lead to average grain sizes that are a steep function of both radius R and height z, as can be seen from the central panels of Fig. 4. In the α = 10-4 case, can be as high as >10 μm in the disk mid-plane. However, in this case steeply decreases to average grain sizes of the order of a few nm at relatively low scale heights. The dust-to-gas ratio and the average grain size jointly determine the total dust surface area per unit volume (see Eqs. (15) and (19)), which at any given location clearly decreases with turbulence. The difference in the upper layers of the disk can be higher than four orders of magnitude for α varying between 10-4 and 10-2. Interestingly, the total dust surface area that the STN model assumes is quite similar to the α = 10-2 case, indicating that the parameters chosen for the STN model, with a settling parameter χ = 0.2, are representative of a quite turbulent disk.

4.2. The effects of grain growth, radial drift, and dust settling on gas properties

The distribution of the dust particles in the disk has multiple effects, which are addressed below. In particular, their concurrence can significantly affect the radial distribution of the CO emission and it is thus important to take them into account when modelling CO emission in a protoplanetary disk. The focus of this section is on the outer disk, which usually dominates the emission both in continuum and in the low-lying CO rotational lines.

4.2.1. Dust and gas thermal structures

The opacities in a disk do depend on the size distribution of dust particles within the disk (see Sect. 2.5). The UV wavelength range is particularly important because it provides one of the most important heating sources for the gas phase via photoelectric heating of both PAHs and small particles (e.g. Weingartner & Draine 2001; Croxall et al. 2012). At low viscosities, the smallest grains become the dominant mass carriers in the outer disk (see discussion in Sect. 4.1) and are generally more abundant than in the STN model, or in the α = 10-2 case (which again look very similar). This has the important effect that for low viscosities, the small grains are very effective in screening the FUV (far-UV) photons (see Appendix A for more details). In these cases, the transition between optically thick and optically thin disk regions at FUV wavelengths becomes very sharp in the vertical direction. However, at longer wavelengths, the exact opposite happens. At lower viscosities, the opacities in the IR (infrared) at intermediate scale heights are much lower than in the more turbulent cases, simply due to the fact that the vertical settling for low α values is more severe. Thus, the IR photons, which are the most important in determining the dust temperature, can penetrate further into the disk, heating the dust to higher temperatures. At intermediate scale heights in the outer disk, close to the 12CO emitting layer, the dust temperatures can be as high as ~100 K in the α = 10-4 case, whereas they reach at most ~50 K in the α = 10-2 case.

thumbnail Fig. 5

Mid-plane dust temperatures of the models using the gas density parameters by de Gregorio-Monsalvo et al. (2013). The grey lines show the sublimation temperature for CO in the α = 10-3 model; the dashed line uses a binding energy of 855 K, the dotted-dashed line of 1100 K.

Open with DEXTER

We now focus on the disk mid-plane, since this is where most of the mass lies, and where important gas-grain chemistry dependent on the dust temperature occurs. The dust temperatures in the disk mid-plane are shown in Fig. 5. While the STN and STN-SM models show a monotonically decreasing temperature (which is usually assumed in most models), grain growth and vertical settling concur in yielding a non monotonic dust temperature. The effect becomes more severe for lower values of turbulence. The reason can be understood as follows. At intermediate radii (100−300 AU), the average size of the dust particles has grown considerably, in particular for low turbulence (see Fig. 3). For such large particles, the aerodynamic drag of the turbulent eddies is quite inefficient and dust particles can easily settle close to the disk mid-plane. However, in the inner regions the gas densities are high enough that the dust particles are stirred to high altitudes (see left panels of Fig. 4). The net effect is that due to settling there is an intermediate region where the disk is self-shadowing (Dullemond & Dominik 2004; Birnstiel et al. 2015). At larger radii the average grain size decreases again, due to radial drift. The abundant small particles are again stirred up easily, making them intercept significant stellar radiation. Thus, the dust temperature increases at large radii (R ≳ 300 AU). The result is that there is a region in the disk mid-plane where, for low turbulent models, the dust temperatures can be lower by 20% with respect to the STN model (Fig. 5). Another way to understand this mechanism is by looking at the radial gradient of the vertically averaged grain size, which in intermediate regions decreases more steeply than the gas surface density profile (Fig. 6). Thus, Eq. (6) implies that there is a region where the average Stokes number decreases with radius, with less dust particles being stirred up towards the upper layers. We note that the disk analysed here is a relatively warm disk, with most of the disk mid-plane too warm for CO freeze-out. Interestingly, the difference between the models is just around the CO freeze-out temperature (see Fig. 5), which is computed as the dust temperature where kads(CO) = kthdes(CO).

The gas temperature has somehow the opposite dependence on disk turbulence. The total dust surface area per volume at intermediate scale heights is much lower for lower turbulence values than in the STN and α = 10-2 models, and this reduces the efficiency of photoelectric heating (Jonkheid et al. 2004, 2007). The low dust surface area in the low turbulence models also implies that the gas-grain collisions are very ineffective in yielding thermal coupling between the dust and gas components. The net effect is that the gas can get significantly colder in low viscosity disks, with gas temperatures <15−20 K in large regions of the disk, even though in those same regions the dust temperature is warmer. In the left panels of Fig. 7, it is apparent that the region of thermal coupling close to the disk mid-plane (highlighted in white in the left panels of Fig. 7) becomes vertically thicker as viscosity increases. More plots on the resulting disk thermal structure are shown in Appendix A.

By comparing the left panels of Fig. 7 with Fig. A.1, it emerges that the region where the gas temperature is lower than the dust temperatures occurs just below the τFUV = 1 layer, where the heating via photoelectric effect is quenched. Moreover, line cooling becomes very significant in the same region. This is apparent if we look at the abundances plots of atomic carbon C and CO, shown in Fig. 7. The region where Tgas/Tdust < 1 in the low turbulence cases (highlighted in blue in the left panels of Fig. 7) coincides with the transition between atomic C and CO, that is when the cooling becomes significant via the lines of these two species.

The same mechanisms can be observed in the vertical cuts in the outer disk (at ~550 AU) shown in Fig. A.3. In these regions of the disk, while the mid-plane temperature is the same for all the models, Tdust gets warmer at z > 0 in the low turbulence cases, as described above. At the same time, the gas temperature becomes considerably colder, due to the lack of photoelectric heating, to the poor thermal coupling, and to the effective cooling in the region of the C+ – C – CO transition. Finally, we confirm the general results of the simpler models by Jonkheid et al. (2004, 2007) and Vasyunin et al. (2011), where they suggest that settling causes the maximum of the abundances in the vertical direction to shift closer to the disk mid-plane. In particular, we note that the difference of this work from Vasyunin et al. (2011) is that in the latter the dust evolution considered the fragmentation limit only, and more importantly the gas thermal balance was not computed.

thumbnail Fig. 6

Vertically averaged grain size () as a function of radius for the turbulent models, compared to the gas surface density profile rescaled to some value Σ0. For both α = 10-3 and α = 10-4 there is an intermediate region where the radial gradient of is steeper than the radial gradient of Σgas, thus leading to disk self-shadowing.

Open with DEXTER

thumbnail Fig. 7

From top to bottom: STN, α = 10-2, 10-3 and 10-4 models. From left to right: ratio Tgas/Tdust, C abundance, and gas phase CO abundance. The yellow and red solid lines in the left and right panels indicate the vertical τ = 1 layer for the 12CO J = 3–2 transition.

Open with DEXTER

thumbnail Fig. 8

Left and central panel: column density of CO gas and CO ice (JCO) of all models. We note the different scale on the y-axis in the two panels: N(JCO) is always much lower than N(CO) in the gas phase at all radii, even for α = 10-3. Right panel: abundance of CO (solid lines) and CO ice (dashed lines) along the disk mid-plane.

Open with DEXTER

thumbnail Fig. 9

From top left to bottom right: CO ice abundance in the STN, α = 10-2, 10-3, and 10-4 models. We note the different vertical scale from other similar plots.

Open with DEXTER

thumbnail Fig. 10

From top left to bottom right: vertical cuts at ~300 AU of the STN, α = 10-2, 10-3, and 10-4 models. Legend: solid lines in black, red, green, and blue represent abundances of CO, C, C+, and CO ice, respectively. Dashed lines in orange and brown are Tgas and Tdust, respectively. We note that Tgas falls below Tdust at low values of α and that this thermal de-coupling starts at the C–CO transition.

Open with DEXTER

4.2.2. CO abundance and snowline

The vertical column density of CO in the gas phase is almost the same in all models, as shown in the left panel of Fig. 8. The only differences appear at very large radii, where the STN and STN-SM models have higher photodissociation rates due to the deeper UV penetration in the disk (see Fig. A.1), which leads to slightly lower column densities. Settling does not affect the gas CO abundances significantly, as already suggested by Jonkheid et al. (2007).

The thermal structure of both dust and gas phases in the disks, and the total dust surface area available for freeze-out and desorption, have a significant impact on the amount of CO ice (see column density of CO ice in the central panel of Fig. 8) and the location of the CO snow surface. The snow surface is defined as the location (in both R and z) where 50% of a species is in the gas phase and 50% is frozen out onto dust grains. Similarly, the snowline is defined as the radius where the same happens along the disk mid-plane. The CO ice abundance of most models is shown in Fig. 9. There are at least two effects that can be identified. The first one is that in the STN model, some CO ice is present out to very large radii (>600 AU), since the dust is colder in the outer regions than in other models. The second one is that the CO ice abundance gets to rather high values (10-4) for α ≤ 10-3. This is mainly due to vertical settling self-shadowing these regions of the disk, as explained in the previous section. In fact there is a clear correlation between the column density of CO ice (Fig. 8) and the dust temperature in the disk mid-plane (Fig. 5). The dust temperatures in these models fall to ~17−20 K in the disk mid-plane, at R ~ 200−350 AU, and the large amount of dust surface area leads to significant CO freeze-out. This is shown more quantitatively in Fig. 10, where the abundances of the main C carriers and the gas and dust temperature are shown in a vertical cut at ~300 AU. The amount of CO ice in the STN and α = 10-2 case are again rather similar. Moreover, the abundance of CO ice clearly reflects the vertical settling of the dust particles, with the abundance being a steep function of z for the low turbulence cases. Finally, we note that for the warm disk considered in this work, the column density of the gas phase CO is very similar in all models, since N(CO) N(JCO) even in the lowest turbulence case (Fig. 8), since the gas phase CO column densities are dominated by the layers above the CO snow surface in all cases.

The low turbulence models predict that CO freezes out in the disk mid-plane only in a very specific radial range. At larger radii, where the growth of dust particles has been limited by radial drift, the higher dust temperatures in the mid-plane (see Fig. 5) enhance thermal desorption, leading to very low abundances of CO ice. This suggests that radial drift, grain growth, and settling can concur in having a second desorption front of CO at large radii, where the average grain size decreases substantially (Cleeves 2016). We thus confirm that radial drift can indeed lead to a thermal inversion in the outer disk, thus leading to a second CO desorption front. Differently from the models by Cleeves (2016), where the second CO desorption front is due to photodesorption, due to a low dust-to-gas ratio in the outer disk, in our models it is caused by thermal desorption, since in the outer regions Tdust becomes higher than the CO sublimation temperature. Interestingly, a few observations suggest that such a second CO desorption front could indeed be present in IM Lup, AS 209, and TW Hya, just outside the submm radius (Öberg et al. 2015; Huang et al. 2016; Schwarz et al. 2016, respectively). Whether these observations can be explained by the models presented here, or whether non-thermal desorption is important in these outer regions, is difficult to say at this stage, in particular given the very different properties of both disk and stellar mass of these three objects from the parameters used in this work.

4.3. Emission

Hitherto we have discussed the differences that grain evolution has on the density and temperature structure, and on CO abundances, in a representative model of a protoplanetary disk. We now investigate how these properties affect observable quantities, in particular the intensity profiles of both (sub-)mm continuum and CO isotopologues.

4.3.1. Continuum

The peak-normalised continuum emission at 850 μm of all models is shown in Fig. 11. The STN-SM model is clearly different from all the others, since it lacks grains larger than 1 μm. Its mm emission is thus the most compact in radial extent. The radial extent of the STN, α = 10-2, and 10-3 cases is quite similar. The STN model does not have a radially varying grain size distribution, that is, radial drift is not considered in this parametric model. However, the sharpness of the continuum emission in the very outer regions does depend on the radial gradient of the grain size distribution (or on the radial gradient of the dust-to-gas ratio, which we did not consider in this paper, see Birnstiel & Andrews 2014). As the turbulence gets lower (and thus the steepness of the radial grain size distribution becomes more significant), the outer edge of the continuum emission becomes sharper. Thus, this sharp structure can be naturally explained by a combination of grain growth and radial drift. Interestingly, these models predict that the sharpness of the edge is correlated to a large scale substructure in the continuum emission that is clearly visible in the α = 10-4 case in Fig. 11, but also present in the α = 10-3 model. We recall that this substructure is due to the accumulation of large grains (and depletion of small particles) outside the fragmentation radius (Birnstiel et al. 2015). For this radial substructure to arise, there is no need of any planet or hydrodynamical instability, it is a pure consequence of radial drift and fragmentation of dust particles.

From our dust evolution models, we also recover the well-known result that the size of the emitting surface in continuum depends on the wavelength, as predicted by grain models considering radial drift. As an example, the right panel of Fig. 11 compares the normalised synthetic emission profiles at 850 μm and 8 mm. Apart from the STN-SM model, which does not consider either grain growth or radial drift, all other cases show more compact emission at longer wavelengths.

thumbnail Fig. 11

Peak-normalised continuum intensity profiles of all models. In all panels the solid lines show the profiles at 850 μm. The difference between the two left panels is the radial scale. The dashed grey line in the two left panels indicates the normalised input surface density profile. In the right panel, the coloured dashed lines show the peak-normalised continuum intensity profiles at 8 mm, which look more compact than the 850 μm profiles, apart from the STN-SM model.

Open with DEXTER

thumbnail Fig. 12

From left to right: intensity profiles of the 12CO, 13CO, and C18O J = 3–2 line. The grey dashed line depicts the input surface density, convolved by the same beam size as the synthetic line intensity profiles. The horizontal dashed-dotted line represents an arbitrary intensity cut, with a dynamic range of 30 for 12CO, and a dynamic range of 10 for 13CO and C18O.

Open with DEXTER

4.3.2. CO versus continuum radial profiles

The emission of CO isotopologues is often used to determine the temperature and density structure of protoplanetary disks. The most abundant 12CO is used to constrain the thermal profile, whereas the more optically thin 13CO and C18O are studied to determine the gas surface density profile (e.g. de Gregorio-Monsalvo et al. 2013; van der Marel et al. 2016; Huang et al. 2016; Cleeves et al. 2016; Schwarz et al. 2016).

The radial intensity profiles of the CO J = 3–2 line are shown in Fig. 12. While the STN and α = 10-2 models are very similar as usual, the α = 10-3 and 10-4 profiles become dimmer in the outer regions of the disks, in particular for 12CO. This can be well understood from what has been outlined in Sect. 4.2.1. For lower turbulence values, the gas becomes colder, in particular in the regions where 12CO is emitting most, not far from the CCO transition. The lower temperature implies that the J = 3 state is less populated, leading to lower emission in the J = 3–2 line. The τ = 1 surface of the J = 3–2 12CO line crosses the region of thermal de-coupling of dust and gas in the low turbulent cases (see Fig. 7), where the gas can be much colder than the dust. This is confirmed by looking at the line ratios, in particular at the ratio of the J = 6–5/J = 2–1 lines, the J = 3–2/J = 2–1 lines, and the J = 3–2/J = 1–0 lines, for both 12CO and 13CO (Fig. 13). The 12CO line ratios clearly show that there is a correlation between the line ratio being steeper and the steeper intensity profile in the J = 3–2 line. This confirms that the difference seen between the models in the profiles of 12CO is due to a different thermal structure, in particular due to the vertical settling (and weak thermal dust-gas coupling) leading to lower temperatures in the outer disk. Such temperature effects are still seen in the 13CO line ratios, even though the trend with turbulence is less clear, due to the lower optical thickness of the lines. This confirms that the discrepancy is not caused by a difference in the column density, which is almost the same in all models (see Fig. 8), but by an excitation effect.

Figure 12 shows another significant result. Given the same gas surface density profile structure (the input gas density structure is the same for all models), the radial intensity profile of all CO isotopologues can be quite different depending on the dust properties of the disk, that is, depending on the radial grain size distribution and vertical settling. In particular, observations performed with ALMA at high sensitivity usually lead to a dynamic range in 12CO of the order of ~30, whereas in 13CO and C18O this is of the order of 10. In all panels of Fig. 12 the horizontal dashed-dotted line shows an arbitrary intensity cut with such dynamic ranges. These plots indicate that an observed disk with the same gas surface density profile would lead to very different estimates of the disk gas outer radius (), since this would be observationally determined by the sensitivity cut. In particular, from the 12CO line, the outer radius would vary between ~350−550 AU for α ranging between 10-2 and 10-4, whereas the same estimate based on the 13CO line would yield disk outer radii varying between ~250−550 AU for the same range of turbulence. Such an effect is less prominent for the C18O isotope, since its lower optical depth makes it depend more strongly on the CO column density. For optically thinner isotopologues, the line emission would look smaller in size, as observed by for example. Isella et al. (2016). In general, for the same gas surface density profile, disks with lower turbulence would look smaller in CO intensity maps.

thumbnail Fig. 13

Line ratios of J = 6–5/J = 2–1 (left column), J = 3–2/J = 2–1 (middle column), and J = 3–2/J = 1–0 (right column) for 12CO and 13CO (top and bottom row, respectively). The legend is the same as in Fig. 12.

Open with DEXTER

Both in the turbulent models and in the STN model, the 12CO emission is more extended than the continuum (compare Figs. 11, 12), in most cases by a factor two, considering a dynamic range of 30 for both intensity profiles. This confirms the early suggestion by Dutrey et al. (1998) and Guilloteau & Dutrey (1998) that the difference in radial extents between gas and dust is (mostly) due to optical depth effects. However, the ratio of and (taken as the radius where the intensity profiles reach a dynamic range of 30 at 850 μm), depends strongly on turbulence. In particular, for the turbulent models shown in Figs. 11, 12, the ratio scales from ~1.4 for α = 10-4 to ~4 for α = 10-2, indicating that the ratio increases with turbulence.

Finally, we note that in this paper we have not considered how the CO emission is affected by potential C and O depletion in the disk, which has been suggested by thermo-chemical models of recent observations (e.g. Kama et al. 2016; Bergin et al. 2016; McClure et al. 2016; van’t Hoff et al. 2017; Miotello et al. 2017).

5. A closer view of HD 163296

As a representative case for our models, we consider the HD 163296 system. As mentioned in the Introduction, this is one of the first systems where the need for grain growth and radial drift has been advocated to explain the different radial intensity profiles in (sub)mm continuum and CO rotational lines, and in particular to explain the steep radial decline of the 850 μm emission. This system has been analysed in continuum and CO emission by a few different groups (de Gregorio-Monsalvo et al. 2013; Rosenfeld et al. 2013; Flaherty et al. 2015; Guidi et al. 2016; Boneberg et al. 2016; Williams & McPartland 2016), who have all exploited the ALMA Science Verification data, program 2011.0.00010.SV. In particular, de Gregorio-Monsalvo et al. (2013) and Rosenfeld et al. (2013) derived a density structure of the disk from the continuum and 12CO data, whereas more recently Williams & McPartland (2016) fitted the more optically thin 13CO and C18O to obtain directly the gas surface density profile. They all used a surface density parametrisation given by Eq. (1), with very similar stellar mass, stellar luminosity, and flaring angle. We thus run models with the best fit parameters determining the gas surface density structure from these three papers. The actual values are reported in Table 1. We run both STN and STN-SM models of the three parametrisations, together with models with turbulence α = 10-2, 10-3, and 10-4. All parameters not specified in Table 1 are kept fixed. All models presented in the previous sections have the gas surface density parameters obtained by de Gregorio-Monsalvo et al. (2013, see Sect. 3. We stress that we do not aim to find a best fit model for HD 163296. The aim is to show how the properties of dust grains affect the emission of CO isotopologues, which are the most used probe for the gas surface density, and whether we can recover the continuum and lines intensity profiles of actual observations within a single framework.

Table 1

Main parameters of HD 163296 used to produce the models in Fig. 16, as derived in the following references: [A] de Gregorio-Monsalvo et al. (2013); [B] Rosenfeld et al. (2013); [C] Williams & McPartland (2016).

First, we check that the models do recover the integrated continuum fluxes observed at (sub)mm wavelengths. In Fig. 14 we show a comparison of the continuum fluxes at 0.85, 1.3, 2.8, 8, 10, and 52 mm of the models using the parametrisation by de Gregorio-Monsalvo et al. (2013) with the continuum fluxes reported in Isella et al. (2007, at 2.8 and Guidi et al. (2016). The STN-SM model under-predicts the flux at all wavelengths, confirming that grain growth is needed to explain the observed emission. The other models are all roughly consistent with the data. All models under-predict the flux at 5.2 cm, but free-free emission is likely dominating at these long wavelengths and probably contributes also at 1 cm (Guidi et al. 2016).

Qi et al. (2015) determined a radial location for the CO snowline at 90 ± 10 AU by modelling the observed N2H+ and C18O emission (see also the earlier work by Mathews et al. 2013, on DCO+). Even though the recent work by van’t Hoff et al. (2017) has shown that the interpretation of N2H+ emission as a CO snowline tracer is more subtle than previously assumed, the C18O emission still suggests that the snowline is at about 90−100 AU. In the models shown in Fig. 8, we obtain a CO snowline around 250 AU, larger by a factor of two in the most optimistic case (α = 10-3). However, all those models assume a CO binding energy of 855 K, typical for pure CO ice (see Sect. 2.6.3 for a more detailed discussion). By running a turbulent model (α = 10-3) with a CO binding energy more appropriate for ice mixtures (1100 K), we obtain a location of the snowline of ~115 AU, which is more consistent with the data (see Fig. 15), as is also suggested by Qi et al. (2015). The location of the snowline corresponds to the radius at which the dust temperature reaches the CO sublimation temperature (see Fig. 5). The column density of gaseous CO is not significantly affected. Interestingly, the least turbulent models predict a second thermal desorption front in the outer regions of the disk, around 500 AU. Deeper ALMA observations should be able to determine whether this prediction is correct, for example by looking for DCO+ emission in these outer regions.

thumbnail Fig. 14

Comparison of continuum fluxes predicted by all models using the gas density parameters from de Gregorio-Monsalvo et al. (2013) with the data from Isella et al. (2007, at 2.8 and Guidi et al. (2016).

Open with DEXTER

Figure 16 shows the observed and modelled peak-normalised intensity profiles of both continuum at 850 μm and of the 12CO and 13CO J = 2–1 lines. The 12CO and 13CO synthetic images are obtained with a beam convolution of 0.68″ × 0.55″ and 0.72″ × 0.57″, respectively, with a PA = 72°. The vertical scale is used to mimic a peak value of ~300σ. With the parameters used by Williams & McPartland (2016), we always obtain a dust disk outer radius that is too large compared to the ALMA data. More interestingly, the continuum observations of the STN models all fail to reproduce the sharp edge in the outer disk, whereas models including dust evolution are more consistent with the data. For low turbulence values, the sharp edge is associated with substructure in the disk continuum emission, with the level of substructure increasing with decreasing turbulence (i.e. with Rfrag shifting towards the inner regions of the disk). The 12CO data are well recovered by many models, in particular with α = 10-2−10-3. When α = 10-4, the temperature gradient along the CO emitting layer becomes too steep, leading to a fast decline in CO emission, which is not observed in the ALMA data.

6. Discussion

6.1. Gas versus dust radial extent

So far, the quantity that has been mostly derived for protoplanetary disks is disk mass, since the observations have been limited by angular resolution and sensitivity. This fact has led the community to mostly look at models of disk integrated CO emission, in particular at the fluxes of different CO isotopologues which can be used to constrain the gaseous disk mass (e.g. Williams & Best 2014; Miotello et al. 2016). However, recent and upcoming ALMA observations allow us to investigate spatially resolved chemical and physical properties of the gaseous component of disks, in particular the surface density profile and the disk outer radius (e.g. Williams & McPartland 2016).

Our models including physically motivated grain properties confirm the early suggestion by Dutrey et al. (1998) and Guilloteau & Dutrey (1998) that the observed different radial extents in gas and dust can be largely explained by the difference in the optical depth of gas versus dust. This feature is retrieved by all models mimicking grain growth and vertical settling (i.e. the STN models), without the need to invoke radial drift. However, a simple quantity such as the disk gas outer radius as derived from 12CO and 13CO radial intensity profiles can depend significantly on the properties of the dust grains. In particular, the radial extent of the (sub)mm continuum does not depend strongly on turbulence (considering all the assumptions in our models), whereas the radial extent of the CO emission does. In our models, the ratio of the CO and mm continuum radius is directly related to the amount of turbulence in the disk, with the ratio increasing with α (Fig. 17). The exact values of the ratio will depend on the dynamic range of the CO emission used to determine . The value of the ratio will also depend on the wavelength considered to extract , since the size of the dust disk can depend on wavelength (Fig. 11), with disk size slowly decreasing with increasing wavelength. In this preliminary study we did not explore the dependence on other important parameters, such as stellar mass (thus stellar spectrum), disk mass, and surface density profile, although qualitatively the trend should be robust.

thumbnail Fig. 15

Abundance of CO (solid lines) and CO ice (dashed lines) along the disk mid-plane for the α = 10-3 model, with CO binding energy of 855 K (blue lines) and 1100 K (brown lines). The vertical dashed-dotted black line shows the location of the CO snowline as determined by Qi et al. (2015). The parametrisation of the model is by de Gregorio-Monsalvo et al. (2013). To reproduce the observed snowline location, the CO binding energy needs to be higher than that of pure CO ice, with Eb = 1100 K leading to a much better agreement than Eb = 855 K.

Open with DEXTER

thumbnail Fig. 16

Comparison between the Science Verification ALMA data and models of HD 163296. The dashed lines represent the ALMA data, the solid lines the models. From top to bottom: STN, α = 10-2, 10-3, and 10-4 models. From left to right: model parameters taken from de Gregorio-Monsalvo et al. (2013), Rosenfeld et al. (2013), and Williams & McPartland (2016). The α = 10-2 case in the middle column is the best representation of the data among our models.

Open with DEXTER

In order to interpret measurements of the disk outer radius and line emission, a proper treatment of the dust properties is needed. In the future, observations at multiple (sub)mm frequencies will allow us to determine some fundamental properties of the large dust grains as a function of radius, in particular meaningful constraints on the radial dependence of the grain size distribution (e.g. Guilloteau et al. 2011; Pérez et al. 2012, 2015; Menu et al. 2014; Tazzari et al. 2016). Such inferred dust properties could then be used as a parametric input for the thermo-chemical models, in order to interpret the emission of important molecules (as CO) more consistently with the continuum observations.

In this work we have mostly looked at the large scale structure of disks, which are assigned a smooth gas surface density profile. We have not explored the potential effects of small substructures in the surface density (as observed by ALMA Partnership et al. 2015; Andrews et al. 2016; Tsukagoshi et al. 2016; Isella et al. 2016) and we have not looked at how substructures in the dust properties can affect the gas emission. However, from the results shown in this paper, it is reasonable to assume that radial substructures in dust properties will affect the CO emission, by changing the gas thermal variations and possibly abundances on the same scales (e.g. Bruderer 2013; Yen et al. 2016b; Isella et al. 2016; Teague et al. 2017).

6.2. Effects of grain growth, vertical settling, and radial drift on chemistry and line emission

The combination of dust evolution processes, in particular grain growth, vertical settling, and radial drift, largely determine the dust density and thermal structure in protoplanetary disks. Moreover, they also affect the continuum mean intensity throughout the disk, which is fundamental in the thermo-chemical processes of the gaseous material. In particular, severe vertical settling and substantial grain growth decrease the efficiency of physical and chemical processes dependent on the dust surface area. In these cases, the gas can be largely colder than the dust component at intermediate scale heights due to the low rate of gas-grain collisions. Optically thick molecular lines emitting from the disk intermediate layers, as low J CO transitions, can show the imprint of such cold gas temperatures by being less excited than expected from simple models. Thus, the largest effect of dust evolution and radial drift is on the gas temperature structure; enhanced UV penetration heating the gas and dissociating or desorbing CO has a minor effect. As a result, the radial CO emission profile can be significantly steeper than the actual gas surface density profile. In this paper we have focused on CO emission lines, but the same can indeed occur for other molecular lines emitted from the same regions of disks.

The dust properties also affect the abundances of ices in the disk mid-plane. The dust density structure and grain size distribution determine the dust temperature in the disk mid-plane and the grain size distribution sets the total dust surface area available for adsorption and desorption processes. We have shown how these two mechanisms determine the amount of CO ice in the disk mid-plane for the disk parametrisation analysed in this paper, which is typical of a Herbig star. The combination of radial drift and vertical settling in particular can induce a thermal inversion in the disk mid-plane in the outer regions of the disk, without the need to reduce the dust-to-gas ratio in these outer regions. Since the ice abundances are so sensitive to the thermal structure and dust total surface area available, including more realistic dust properties in the thermo-chemical models is important to follow the chemical evolution of some species, in particular for processes depending on the dust surface area, such as grain surface reactions in the disk mid-plane. The formation of some precursors of complex molecules, such as CH3OH and H2CO (the latter has been detected in HD 163296 by Yen et al. 2016a), is indeed coupled to the ice composition and abundances on the mantles of dust grains.

thumbnail Fig. 17

Ratio derived for the models using the disk parametrisation of de Gregorio-Monsalvo et al. (2013), with a 0.52″ synthetic beam. The quantity is derived assuming a dynamic range of 30 for 12CO and a dynamic range of 10 for 13CO and C18O. For the continuum at 850 μm, the assumed dynamic range is 30. For all isotopologues, the ratio is a steep increasing function of disk turbulence.

Open with DEXTER

6.3. Dust growth versus freeze-out timescales

In the results shown above, the main implicit assumption is that the growth timescale of dust particles is shorter than the main chemical timescales for the reactions considered in our chemical network, such that the thermal balance can be computed on static dust properties. This assumption is made only for numerical reasons, that is, there is still no code capable of coupling dust evolution and thermo-chemistry in a 2D disk for a viscous timescale. We can however check a posteriori whether the methodology presented here should be improved in the near future (see Haworth et al. 2016b, for a more generic discussion). We can do so by comparing the growth timescale of dust particles of a given size with the CO freeze-out timescale on particles of the same size, in order to focus in particular on the CO snow surface. The model that is best suited to make such comparison is the STN-SM one, where the grains are still small (m). In the comparison, the gas thermal structure computed by the model is used.

The growth timescale tgrowth of a dust particle of size a due to turbulent motions can be written as tgrowth = a/ȧ, where ȧ is (Brauer et al. 2008): (28)The relative velocities due to turbulence can be expressed as (Ormel & Cuzzi 2007): (29)where the local Stokes number can be written as: (30)We thus obtain that: (31)The quantity cs is obtained from the STN-SM model, where , and the contribution of chemical species to both P and ρgas is taken into account.

thumbnail Fig. 18

Ratio of the growth timescale of dust particles and the freeze-out timescale of 12CO. For both timescales, particles with sizes of 0.1 μm have been considered. The thermal structure is taken from the STN-SM model. A turbulent α = 10-3 has been used in the calculations (Eq. (32)).

Open with DEXTER

Similarly, the freeze-out timescale can be expressed as tads ~ 1 /kads (see Eq. (21)), where kads is the adsorption rate. If we consider an ensemble of grains with size a for the freeze-out timescale, the ratio of the two timescales can be written as: (32)For the most abundant 12CO, M(CO) = 28.

Figure 18 shows the ratio of the two timescales for particles of size a = 0.1 μm, an assumed turbulent α = 10-3, and ρgr = 2.5 g cm-3. All the other quantities are taken from the STN-SM model. In the whole disk, tgrowth is larger than tads, and in particular in the mid-plane of the inner disk, tgrowthtads. This implies that grains should be growing once the CO gas has already frozen onto dust grains, where the physical conditions of the disk allow this to happen. This fact strongly supports the hypothesis that in order to determine the abundance of ices in a protoplanetary disk, more sophisticated models, with coupled dynamical grain growth and thermo-chemistry, are needed. Some recent papers are indeed going in this direction, coupling different aspects of dust evolution (e.g. radial drift, grain growth, or vertical settling) with the chemical evolution of a disk (e.g. Piso et al. 2015, 2016; Krijt et al. 2016; Bergin et al. 2016; Kama et al. 2016).

7. Conclusions

In this work we have coupled dust evolution models to thermo-chemical models of protoplanetary disks, focusing in particular on the radial properties of continuum and CO rotational lines emission. Detailed results are presented for a model representative of the HD 163296 disk and compared with observations of that disk. Our main conclusions can be summarised as follows:

  • 1.

    Differences in gas and dust radial extents of protoplanetarydisks are readily found in our models and can be largely explainedby the difference in optical depth of gas versus dust, without theneed to invoke radial drift (the STN model leads to a similar size incontinuum as the turbulent models).

  • 2.

    Different turbulence values can dramatically affect the estimate of the disk gas outer radius. The gas outer radius probed by 12CO emission can differ by a factor of ~ three for a turbulent α ranging between 10-4 and 10-2, with the ratio of the CO and mm radius increasing with turbulence (Fig. 17).

  • 3.

    For the massive disk considered in this paper, 13CO exhibits the same trend as 12CO, with the outer radius increasing with turbulence. The intensity profile of the more optically thin C18O depends more shallowly on turbulence.

  • 4.

    The effect of radial drift is primarily visible in the sharp outer edge of the (sub)mm continuum intensity profile, in particular for low values of turbulence, when the maximum attained dust grain size is a steep function of radius. This may also lead to radial substructure in the continuum emission at large radii.

  • 5.

    Proper treatment of dust evolution via grain growth, radial drift, and vertical settling leads to gas outer disks that are significantly colder than in models with simpler parametrised dust treatments. This is caused mostly by settled large grains coupling poorly with the gas, with the gas temperature even falling below the dust temperature at intermediate disk heights. This low gas temperature should be reflected in the CO line intensity ratios. Enhanced penetration of UV radiation has a smaller effect on CO intensities.

  • 6.

    Lower turbulence allows particles to grow to larger sizes. Together with low vertical stirring, this leads to severe settling. As a consequence, steeper CO intensity profiles than the actual surface density profile are obtained for low α, due to a steeper gas temperature gradient in the intermediate disk layers.

  • 7.

    Radial drift and dust settling concur in causing a thermal inversion of the dust, potentially leading to a second mid-plane CO thermal (rather than photo) desorption front at large radii.

  • 8.

    Using HD 163296 as a test case, we are able to obtain good agreement between our models and the (sub)mm continuum and 12CO – 13CO lines intensity profiles for fairly high values of α.

  • 9.

    The CO snowline location of HD 163296 at ~90 AU, as determined by Qi et al. (2015), can be retrieved by using a binding energy for CO typical of ice mixtures (Eb ~ 1100 K). Lower binding energies, in particular for pure CO ice (Eb ~ 855 K) fail to reproduce the CO snowline location.

In summary, models mimicking grain growth and vertical settling with a parametrised dust treatment but no radial drift reproduce several of the observed features well, but fail to explain sharp dust edges and result in different radial and vertical gas temperature structures. Using such models to infer the underlying gas surface density structure from observed radial profiles could therefore also lead to incorrect conclusions.

Acknowledgments

We thank Ruud Visser, Paola Pinilla, Arthur Bosman, Merel van’t Hoff, and Paolo Cazzoletti for useful discussions, and Michiel Hogerheijde and the Allegro team, and especially Nico Salinas, for providing the cleaned ALMA data. We thank the anonymous referee for providing comments that helped improve the clarity of the paper. This paper makes use of the following ALMA data: ADS/JAO.ALMA#2011.0.00010.SV. ALMA is a partnership of ESO (representing its member states), NSF (USA), and NINS (Japan), together with NRC (Canada) and NSC and ASIAA (Taiwan), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO, and NAOJ. T.B. acknowledges funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme under grant agreement No. 714769. Astrochemistry in Leiden is supported by the European Union A-ERC grant 291141 CHEMPLAN, by the Netherlands Research School for Astronomy (NOVA), and by a Royal Netherlands Academy of Arts and Sciences (KNAW) professor prize. All the figures were generated with the python-based package matplotlib (Hunter 2007).

References

Appendix A: Disk thermal structure

Figure A.1 shows the attenuated UV field in STN and turbulent models. Figure A.2 portrays the dust and gas thermal structure of the same models. Finally, Fig. A.3 shows the gas and dust temperatures, and the abundances of CO, C, and C+, in a vertical cut at 550 AU.

thumbnail Fig. A.1

From top left to bottom right: STN, α = 10-2, 10-3, and 10-4 model. Colours show the UV field in G0 units, where G0 ~ 2.7 × 10-3erg s-1 cm-2 is the UV interstellar radiation field (ISRF) between 911 Å and 2067 Å (Draine 1978). The cyan dashed lines indicate the τFUV = 1 surface.

Open with DEXTER

thumbnail Fig. A.2

From top to bottom: STN, α = 10-2, 10-3, and 10-4 models. From left to right: dust temperature Tdust, gas temperature Tgas, and ratio of the two.

Open with DEXTER

thumbnail Fig. A.3

From top left to bottom right: vertical cuts at ~550 AU of the STN, α = 10-2, 10-3, and 10-4 models. Legend: solid lines in black, red, and green represent abundances of CO, C, and C+, respectively. Dashed lines in orange and brown are Tgas and Tdust, respectively.

Open with DEXTER

Appendix B: Line profiles of CO isotopologues

The spatially unresolved line profiles of three CO isotopologues discussed and shown in Sect. 4.3.2 already contain important information. The profiles of the J = 3–2 transition are reported in Fig. B.1. The first result is that the total fluxes do not change significantly between different models, that is, settling and grain growth do not affect the total CO flux (e.g. Miotello et al. 2014; Woitke et al. 2016). Moreover, the optically thinner isotopologues (in particular C18O) show almost perfect agreement between all models. For optically thicker CO isotopes, the differences between the models become more prominent. This suggests that the discrepancy is not caused by a difference in the column density, which is almost the same in all models (see Fig. 8), but by an excitation effect, that is, by a difference in the gas temperature. In particular, the emission from the outer regions (thus at low velocities) is clearly correlated to the disk turbulence, with the peaks of the double-horned line profile being brighter for higher viscosity. Moreover, the 12CO line profile shows brighter high velocity wings for lower viscosities, suggesting higher temperatures from the very inner regions.

thumbnail Fig. B.1

From left to right: integrated line profile of the 12CO, 13CO, and C18O J = 3−2 transition of the different models shown in Fig. 12.

Open with DEXTER

Appendix C: High resolution synthetic images

All the synthetic intensity profiles shown in the paper use a convolution beam of 0.52″ × 0.38″ (Sect. 3). However, ALMA has and is going to provide higher angular resolution observations of protoplanetary disks. We thus show the peak-normalised synthetic intensity profiles of both continuum and CO lines of our main models convolved with a smaller circular beam (0.1″ resolution) in Figs. C.1, C.2. With higher angular resolution, the substructure outside the fragmentation radius is more prominent. Moreover, the difference of the sharpness of the outer edge between the STN model and models including dust evolution is even more significant. The models with α ≤ 10-3 show a rather sharp edge in the 850 μm intensity profiles.

The intensity profiles of the line emission are similar to the low angular resolution ones. The main difference appears in the inner disk (see top panels of Fig. C.2), where part of the line flux is lost due to high optical thickness of the continuum. This is apparent in the α = 10-4 case, where the inner disk is heavily populated by large grains, due to low turbulent velocities (see Fig. 3).

In Fig. C.3 we show the peak normalised moment 0 maps of the 12CO J = 3–2 line of the models of HD 163296, where the parameters by de Gregorio-Monsalvo et al. (2013) are used for the gas density structure. While the dust radial extent seems similar for all models, expect the STN-SM one, the gas radial extent depends significantly on α, with the disk appearing smaller for lower turbulent parameters. The gas extents are derived for a fixed dynamic range of 100 in the line emission.

thumbnail Fig. C.1

As in Fig. 11, with a 0.1″ beam convolution.

Open with DEXTER

thumbnail Fig. C.2

As in Fig. 12, with a 0.1″ beam convolution. Top panels show the same plots as in the bottom panels, but with a different radial scale.

Open with DEXTER

thumbnail Fig. C.3

Moment 0 map of the 12CO J = 3–2 line of HD 163296 parametrised with the prescription by de Gregorio-Monsalvo et al. (2013). The moment 0 map is shown in log scale, normalised to peak value, with dynamic range equalling 100 in the emission. Contours show peak normalised continuum levels at 850 μm, at peak value over 2, 4, 8, 16, 32, 64, and 128. The images are shown with an angular resolution 0.1″.

Open with DEXTER

All Tables

Table 1

Main parameters of HD 163296 used to produce the models in Fig. 16, as derived in the following references: [A] de Gregorio-Monsalvo et al. (2013); [B] Rosenfeld et al. (2013); [C] Williams & McPartland (2016).

All Figures

thumbnail Fig. 1

Diagram of the method used in this paper. The blue box indicates the input parameters and the green box indicates the output observables. Boxes with orange borders show new modules, compared to the standard DALI code. Orange text underlines modifications in previous DALI modules. The red text lists input parameters needed to compute the dust distribution.

Open with DEXTER
In the text
thumbnail Fig. 2

Gas density structure used in all models.

Open with DEXTER
In the text
thumbnail Fig. 3

Grain size distribution for α = 10-2, 10-3 and 10-4, from left to right, respectively. The two rows show the same plots in two different radial scales. The dashed horizontal line is the total dust-to-gas ratio Δdg = 0.01, whereas the dashed-dotted vertical line indicates the fragmentation radius Rfrag.

Open with DEXTER
In the text
thumbnail Fig. 4

From top to bottom: STN, α = 10-2, 10-3 and 10-4 models. From left to right: local dust-to-gas ratio δdg, average grain size and total dust surface area per volume ngrσ used in the thermo-chemical module. Contour levels are indicated in the colour bar.

Open with DEXTER
In the text
thumbnail Fig. 5

Mid-plane dust temperatures of the models using the gas density parameters by de Gregorio-Monsalvo et al. (2013). The grey lines show the sublimation temperature for CO in the α = 10-3 model; the dashed line uses a binding energy of 855 K, the dotted-dashed line of 1100 K.

Open with DEXTER
In the text
thumbnail Fig. 6

Vertically averaged grain size () as a function of radius for the turbulent models, compared to the gas surface density profile rescaled to some value Σ0. For both α = 10-3 and α = 10-4 there is an intermediate region where the radial gradient of is steeper than the radial gradient of Σgas, thus leading to disk self-shadowing.

Open with DEXTER
In the text
thumbnail Fig. 7

From top to bottom: STN, α = 10-2, 10-3 and 10-4 models. From left to right: ratio Tgas/Tdust, C abundance, and gas phase CO abundance. The yellow and red solid lines in the left and right panels indicate the vertical τ = 1 layer for the 12CO J = 3–2 transition.

Open with DEXTER
In the text
thumbnail Fig. 8

Left and central panel: column density of CO gas and CO ice (JCO) of all models. We note the different scale on the y-axis in the two panels: N(JCO) is always much lower than N(CO) in the gas phase at all radii, even for α = 10-3. Right panel: abundance of CO (solid lines) and CO ice (dashed lines) along the disk mid-plane.

Open with DEXTER
In the text
thumbnail Fig. 9

From top left to bottom right: CO ice abundance in the STN, α = 10-2, 10-3, and 10-4 models. We note the different vertical scale from other similar plots.

Open with DEXTER
In the text
thumbnail Fig. 10

From top left to bottom right: vertical cuts at ~300 AU of the STN, α = 10-2, 10-3, and 10-4 models. Legend: solid lines in black, red, green, and blue represent abundances of CO, C, C+, and CO ice, respectively. Dashed lines in orange and brown are Tgas and Tdust, respectively. We note that Tgas falls below Tdust at low values of α and that this thermal de-coupling starts at the C–CO transition.

Open with DEXTER
In the text
thumbnail Fig. 11

Peak-normalised continuum intensity profiles of all models. In all panels the solid lines show the profiles at 850 μm. The difference between the two left panels is the radial scale. The dashed grey line in the two left panels indicates the normalised input surface density profile. In the right panel, the coloured dashed lines show the peak-normalised continuum intensity profiles at 8 mm, which look more compact than the 850 μm profiles, apart from the STN-SM model.

Open with DEXTER
In the text
thumbnail Fig. 12

From left to right: intensity profiles of the 12CO, 13CO, and C18O J = 3–2 line. The grey dashed line depicts the input surface density, convolved by the same beam size as the synthetic line intensity profiles. The horizontal dashed-dotted line represents an arbitrary intensity cut, with a dynamic range of 30 for 12CO, and a dynamic range of 10 for 13CO and C18O.

Open with DEXTER
In the text
thumbnail Fig. 13

Line ratios of J = 6–5/J = 2–1 (left column), J = 3–2/J = 2–1 (middle column), and J = 3–2/J = 1–0 (right column) for 12CO and 13CO (top and bottom row, respectively). The legend is the same as in Fig. 12.

Open with DEXTER
In the text
thumbnail Fig. 14

Comparison of continuum fluxes predicted by all models using the gas density parameters from de Gregorio-Monsalvo et al. (2013) with the data from Isella et al. (2007, at 2.8 and Guidi et al. (2016).

Open with DEXTER
In the text
thumbnail Fig. 15

Abundance of CO (solid lines) and CO ice (dashed lines) along the disk mid-plane for the α = 10-3 model, with CO binding energy of 855 K (blue lines) and 1100 K (brown lines). The vertical dashed-dotted black line shows the location of the CO snowline as determined by Qi et al. (2015). The parametrisation of the model is by de Gregorio-Monsalvo et al. (2013). To reproduce the observed snowline location, the CO binding energy needs to be higher than that of pure CO ice, with Eb = 1100 K leading to a much better agreement than Eb = 855 K.

Open with DEXTER
In the text
thumbnail Fig. 16

Comparison between the Science Verification ALMA data and models of HD 163296. The dashed lines represent the ALMA data, the solid lines the models. From top to bottom: STN, α = 10-2, 10-3, and 10-4 models. From left to right: model parameters taken from de Gregorio-Monsalvo et al. (2013), Rosenfeld et al. (2013), and Williams & McPartland (2016). The α = 10-2 case in the middle column is the best representation of the data among our models.

Open with DEXTER
In the text
thumbnail Fig. 17

Ratio derived for the models using the disk parametrisation of de Gregorio-Monsalvo et al. (2013), with a 0.52″ synthetic beam. The quantity is derived assuming a dynamic range of 30 for 12CO and a dynamic range of 10 for 13CO and C18O. For the continuum at 850 μm, the assumed dynamic range is 30. For all isotopologues, the ratio is a steep increasing function of disk turbulence.

Open with DEXTER
In the text
thumbnail Fig. 18

Ratio of the growth timescale of dust particles and the freeze-out timescale of 12CO. For both timescales, particles with sizes of 0.1 μm have been considered. The thermal structure is taken from the STN-SM model. A turbulent α = 10-3 has been used in the calculations (Eq. (32)).

Open with DEXTER
In the text
thumbnail Fig. A.1

From top left to bottom right: STN, α = 10-2, 10-3, and 10-4 model. Colours show the UV field in G0 units, where G0 ~ 2.7 × 10-3erg s-1 cm-2 is the UV interstellar radiation field (ISRF) between 911 Å and 2067 Å (Draine 1978). The cyan dashed lines indicate the τFUV = 1 surface.

Open with DEXTER
In the text
thumbnail Fig. A.2

From top to bottom: STN, α = 10-2, 10-3, and 10-4 models. From left to right: dust temperature Tdust, gas temperature Tgas, and ratio of the two.

Open with DEXTER
In the text
thumbnail Fig. A.3

From top left to bottom right: vertical cuts at ~550 AU of the STN, α = 10-2, 10-3, and 10-4 models. Legend: solid lines in black, red, and green represent abundances of CO, C, and C+, respectively. Dashed lines in orange and brown are Tgas and Tdust, respectively.

Open with DEXTER
In the text
thumbnail Fig. B.1

From left to right: integrated line profile of the 12CO, 13CO, and C18O J = 3−2 transition of the different models shown in Fig. 12.

Open with DEXTER
In the text
thumbnail Fig. C.1

As in Fig. 11, with a 0.1″ beam convolution.

Open with DEXTER
In the text
thumbnail Fig. C.2

As in Fig. 12, with a 0.1″ beam convolution. Top panels show the same plots as in the bottom panels, but with a different radial scale.

Open with DEXTER
In the text
thumbnail Fig. C.3

Moment 0 map of the 12CO J = 3–2 line of HD 163296 parametrised with the prescription by de Gregorio-Monsalvo et al. (2013). The moment 0 map is shown in log scale, normalised to peak value, with dynamic range equalling 100 in the emission. Contours show peak normalised continuum levels at 850 μm, at peak value over 2, 4, 8, 16, 32, 64, and 128. The images are shown with an angular resolution 0.1″.

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.