Different dust and gas radial extents in protoplanetary disks: consistent models of grain growth and CO emission

TL;DR

This study models the dust and gas distributions in protoplanetary disks, showing that differences in their radial extents are mainly due to optical depth effects, with implications for interpreting ALMA observations.

Contribution

It introduces a comprehensive model combining dust evolution and radiative transfer to explain dust and gas size discrepancies without solely relying on radial drift.

Findings

Gas and dust sizes differ mainly due to optical depth effects.

Grain growth and settling thermally decouple gas and dust components.

Model successfully reproduces CO snowline in HD 163296 using ice mixture binding energy.

Abstract

ALMA observations of protoplanetary disks confirm earlier indications that there is a clear difference between the dust and gas radial extents. The origin of this difference is still debated, with both radial drift of the dust and optical depth effects suggested in the literature. In this work, the feedback of realistic dust particle distributions onto the gas chemistry and molecular emissivity is investigated, with a particular focus on CO isotopologues. The radial dust grain size distribution is determined using dust evolution models that include growth, fragmentation and radial drift. A new version of the code DALI is used to take into account how dust surface area and density influence the disk thermal structure, molecular abundances and excitation. The difference of dust and gas radial sizes is largely due to differences in the optical depth of CO lines and millimeter continuum, without the need to invoke radial drift. The effect of radial drift is primarily visible in the sharp outer edge of the continuum intensity profile. The gas outer radius probed by $^{12}$CO emission can easily differ by a factor of $\sim 2$ between the models for a turbulent $\alpha$ ranging between typical values. Grain growth and settling concur in thermally decoupling the gas and dust components, due to the low collision rate with large grains. As a result, the gas can be much colder than the dust at intermediate heights, reducing the CO excitation and emission, especially for low turbulence values. Also, due to disk mid-plane shadowing, a second CO thermal desorption (rather than photodesorption) front can occur in the warmer outer mid-plane disk. The models are compared to ALMA observations of HD 163296 as a test case. In order to reproduce the observed CO snowline of the system, a binding energy for CO typical of ice mixtures needs to be used rather than the lower pure CO value.

Figures

Tables

Table 1. Table 1: Main parameters of HD163296 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 ) .

$M_{\mathrm{disk}}\left(M_{\odot }\right)$	$R_c$ (AU)	$\gamma$	Reference
$7\times {10}^{-2}$	$125$	$0.9$	A
$9\times {10}^{-2}$	$115$	$0.8$	B
$4.8\times {10}^{-2}$	$213$	$0.39$	C

Equations

$${\Sigma }_{\mathrm{gas}}(R)={\Sigma }_{\mathrm{c}}{\left(\frac{R}{R_{\mathrm{c}}}\right)}^{-\gamma }\exp \left[-{\left(\frac{R}{R_{\mathrm{c}}}\right)}^{2-\gamma }\right],$$

$${\rho }_{\mathrm{gas}}(R,z)=\frac{{\Sigma }_{\mathrm{gas}}(R)}{\sqrt{2\pi }Rh}\exp \left[-\frac{1}{2}{\left(\frac{z}{Rh}\right)}^2\right],$$

$$\frac{R_{\rm frag}}{100\,{\rm AU}}\sim\frac{M_{*}}{M_{\odot}}\frac{\alpha}{10^{-3}}\frac{\Delta_{\rm dg}}{0.01}\Bigg{(}\frac{|\gamma|}{2.75}\Bigg{)}^{-1}\Bigg{(}\frac{v_{\rm frag}}{10\,{\rm m\,s}^{-1}}\Bigg{)}.$$

$$D_{\mathrm{gas}}\sim \nu =\frac{\alpha c_{\mathrm{s}}^2}{\Omega },$$

$$\mathrm{Sc}\equiv \frac{D_{\mathrm{gas}}}{D_{\mathrm{dust}}}=1+{\mathrm{St}}^2,$$

$$\mathrm{St}=\frac{a{\rho }_{\mathrm{gr}}}{{\Sigma }_{\mathrm{gas}}}\frac{\pi }{2}.$$

$$\frac{\partial {\rho }_{\mathrm{dust}}(z,a)}{\partial t}=D_{\mathrm{dust}}\left[{\rho }_{\mathrm{gas}}\frac{\partial }{\partial z}\left(\frac{{\rho }_{\mathrm{dust}}}{{\rho }_{\mathrm{gas}}}\right)\right]+\frac{\partial }{\partial z}({\Omega }^2{t}_{\mathrm{fric}}{\rho }_{\mathrm{dust}}z),$$

$$t_{\mathrm{fric}}=\frac{{\rho }_{\mathrm{gr}}}{{\rho }_{\mathrm{gas}}}\frac{a}{v_{\mathrm{th}}},$$

$$v_{\mathrm{th}}=\sqrt{\frac{8k_{\mathrm{B}}T}{\pi \mu m_{\mathrm{H}}}}=\sqrt{\frac{8}{\pi }}{c}_{\mathrm{s}}.$$

$${\int }_{\infty }^{\infty }{\rho }_{\mathrm{dust}}(R,z,a)dz={\Sigma }_{\mathrm{dust}}(R,a).$$

$${\rho }_{\mathrm{dust}}(R,z,a)=\frac{4}{3}\pi {\rho }_{\mathrm{gr}}{a}^{3}f(R,z,a)da,$$

$$⟨a^2⟩(R,z)=\int a^{2}f(R,z,a)da;$$

$$⟨a^3⟩(R,z)=\int a^{3}f(R,z,a)da.$$

$$\bar{a}(R,z)=\frac{⟨a^3⟩(R,z)}{⟨a^2⟩(R,z)}.$$

$$n_{\mathrm{gr}}⟨\sigma ⟩(R,z)=\frac{{\rho }_{\mathrm{dust}}(R,z)}{{\rho }_{\mathrm{gr}}}\frac{3}{4\bar{a}(R,z)},$$

$${\rho }_{\mathrm{dust}}(R,z)={\int }_{a_{\min }}^{a_{\max }}\frac{4}{3}\pi {\rho }_{\mathrm{gr}}{a}^{3}f(R,z,a)da.$$

$${\kappa }_{\nu }(R,z)=\frac{\sum _i{\kappa }_{\nu }(a_i){\rho }_{\mathrm{dust}}(R,z,a_i)}{\sum _i{\rho }_{\mathrm{dust}}(R,z,a_i)}.$$

$${\Lambda }_{\mathrm{gd}}=$$

$$S_{\mathrm{d}}=\frac{n_{\mathrm{gr}}⟨\sigma ⟩}{n_{\mathrm{gas}}}\propto \frac{{\delta }_{\mathrm{dg}}}{\bar{a}}.$$

$$k({\mathrm{H}}_2)=\frac{1}{2}{n}_{\mathrm{H}}^{-1}{v}_{\mathrm{th}}{n}_{\mathrm{gr}}⟨\sigma ⟩{ϵ}_{{\mathrm{H}}_2}{S}_{\mathrm{H}}(T_{\mathrm{dust}}),$$

$$k_{\mathrm{ads}}(X)=n_{\mathrm{gr}}⟨\sigma ⟩\frac{v_{\mathrm{th}}}{\sqrt{M(X)}},$$

$$k_{\rm thdes}=4n_{\rm gr}\langle\sigma\rangle A(X)f(X)\exp{\Bigg{[}-\frac{E_{\rm b}(X)}{k_{\rm B}T_{\rm dust}}\Bigg{]}},$$

$$f(X)=\frac{1}{\max (n_{\mathrm{ice}},N_{\mathrm{ss}}{n}_{\mathrm{gr}}⟨\sigma ⟩)},$$

$$k_{\mathrm{phdes}}=n_{\mathrm{gr}}⟨\sigma ⟩f(X)Y(X)F_{\mathrm{UV}},$$

$$k_{\mathrm{hydro}}(X)=n_{\mathrm{gr}}⟨\sigma ⟩\frac{v_{\mathrm{th}}}{\max (n_{\mathrm{hydro}},N_{\mathrm{ss}}{n}_{\mathrm{gr}}⟨\sigma ⟩)},$$

$${\rho }_{\mathrm{dust},\mathrm{small}}(R,z)=\frac{(1-f){\Delta }_{\mathrm{dg}}{\Sigma }_{\mathrm{gas}}(R)}{\sqrt{2\pi }Rh}\exp \left[-\frac{1}{2}{\left(\frac{z}{Rh}\right)}^2\right];$$

$${\rho }_{\mathrm{dust},\mathrm{large}}(R,z)=\frac{f{\Delta }_{\mathrm{dg}}{\Sigma }_{\mathrm{gas}}(R)}{\sqrt{2\pi }R\chi h}\exp \left[-\frac{1}{2}{\left(\frac{z}{R\chi h}\right)}^2\right].$$

$$\dot{a}=\frac{{\rho }_{\mathrm{dust}}}{{\rho }_{\mathrm{gr}}}\Delta v.$$

$$\Delta v=c_{\mathrm{s}}\sqrt{\frac{3}{2}\alpha \mathrm{St}},$$

$$\mathrm{St}=\frac{a{\rho }_{\mathrm{gr}}}{{\rho }_{\mathrm{dust}}{v}_{\mathrm{th}}}\Omega .$$

Full text

*Aims. *In thermo-chemical models, the dust properties are usually prescribed by simple parametrisations. In this work, the feedback of more realistic dust particle distributions onto the gas chemistry and molecular emissivity is investigated, with a particular focus on CO isotopologues.

*Methods. *The radial dust grain size distribution is determined using dust evolution models that include growth, fragmentation, and radial drift for a given static gas density structure. The vertical settling of dust particles is computed in steady-state. A new version of the code DALI is used to take into account how dust surface area and density influence the disk thermal structure, molecular abundances, and excitation. Synthetic images of both continuum thermal emission and low $J$ CO isotopologues lines are produced.

*Results. *The difference of dust and gas radial sizes is largely due to differences in the optical depth of CO lines and millimeter continuum, without the need to invoke radial drift. The effect of radial drift is primarily visible in the sharp outer edge of the continuum intensity profile. The gas outer radius probed by 12CO emission can easily differ by a factor of $\sim$ two between the models for a turbulent $\alpha$ ranging between $10^{-4}$ and $10^{-2}$, with the ratio of the CO and mm radius $R_{\rm CO}^{\rm out}/R_{\rm mm}^{\rm out}$ increasing with turbulence. Grain growth and settling concur in thermally decoupling the gas and dust components, due to the low collision rate with large grains. As a result, the gas can be much colder than the dust at intermediate heights, reducing the CO excitation and emission, especially for low turbulence values. Also, due to disk mid-plane shadowing, a second CO thermal desorption (rather than photodesorption) front can occur in the warmer outer mid-plane disk. The models are compared to ALMA observations of HD 163296 as a test case. In order to reproduce the observed CO snowline of the system, a binding energy for CO typical of ice mixtures, with $E_{\rm b}\geq 1100$ K, needs to be used rather than the lower pure CO value.

*Conclusions. *The difference between observed gas and dust extent is largely due to optical depth effects, but radial drift and grain size evolution also affect the gas and dust emission in subtle ways. In order to properly infer fundamental quantities of the gaseous component of disks, such as disk outer radii and gas surface density profiles, simultaneous modelling of both dust and gas observations including dust evolution is needed.

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 $\Sigma(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 $\alpha$-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 Sec. 2 and the setup parameters for the models investigated are given in Sec. 3. In Sec. 4 we describe our findings, which are then discussed in Sec. 6. In Sec. 7 we summarise the results and draw our conclusions.

2 Method

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 $\Sigma_{\rm gas}$ dependence on cylindrical radius $R$ (Andrews et al., 2009; Bruderer, 2013):

[TABLE]

where $R_{\rm c}$ is the characteristic radius determining the radial length scale of the disk, and $\Sigma_{\rm c}$ determines the total mass of the disk. Assuming hydrostatic equilibrium, isothermality in the vertical direction, and $z\ll R$ (where $z$ is the vertical coordinate), we obtain the gas mass density $\rho_{\rm gas}$:

[TABLE]

where $h=h_{\rm c}(R/R_{\rm c})^{\psi}$, $\psi$ is the flaring index, and $h_{\rm c}$ 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 $T_{\rm eff}$. 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 $10000$ K, with total luminosity equalling the total potential energy per time unit emitted on the stellar surface by a mass accretion rate $\dot{M}$ (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 $\Delta 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 $v_{\rm frag}$, where we assume that perfect sticking occurs when $\Delta v<v_{\rm frag}$. 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 $R_{\rm frag}$, 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 $R_{\rm frag}$ 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 $p_{\rm d}$ parameter equalling $2.5$ is assumed, where $p_{\rm d}$ is the power law index of the radial dependence of a given grain size bin in the outer disk (see Sec. 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 $\Delta_{\rm dg}$, the mass of the star, and the dust temperature profile, which is implicitly set by the flaring angle $\psi$. Besides these quantities, which are all implemented as initial input parameters (see Sec. 2.1), the output of the grain size reconstruction routine depends on the fragmentation velocity, the disk turbulence, and the grain mass density $\rho_{\rm 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 $v_{\rm frag}=10$ m s*-1*. We parametrise turbulence via the dimensionless parameter $\alpha$ (Shakura & Sunyaev, 1973), and we will vary this parameter in the simulations presented below. Finally, we assume compact dust grains, with $\rho_{\rm gr}=2.5\,$g cm*-3*, and assume the dust-to-gas ratio to be uniform in the whole disk, with $\Sigma_{\rm dust}(R)/\Sigma_{\rm gas}(R)=\Delta_{\rm 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):

[TABLE]

Given the input parameters defined in Sec. 2.1 and the additional ones listed in this Section, a radially dependent grain size distribution $\Sigma_{\rm 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 $a_{\rm min}=50\,\AA$, and a maximum grain size $a_{\rm max}=1$ m.

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:

[TABLE]

where $\nu$ is the kinematic viscosity (Shakura & Sunyaev, 1973), and $\Omega=\sqrt{GM_{*}/R^{3}}$.

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:

[TABLE]

where the Stokes number ${\rm St}$ in the Epstein regime in the disk mid-plane can be computed as (e.g. Birnstiel et al., 2011):

[TABLE]

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):

[TABLE]

where the friction timescale of a particle of size $a$ in the Epstein regime is:

[TABLE]

and the thermal velocity of the gas $v_{\rm th}$ is:

[TABLE]

The temperature $T$ (or better, the sound speed $c_{\rm s}$) used in Eq. 9 is derived from the flaring angle of the disk, which is an input parameter of the models. In particular, $c_{\rm s}=hR/\Omega$ (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 Equation 7, and solve the advection-diffusion equation for $\rho_{\rm d}(R,z,a)$ with a 1D integration in $z$. We normalise each bin by requiring that:

[TABLE]

We note that the “canonical” grain size distribution $f(R,z,a)$ is related to $\rho_{\rm dust}(R,z,a)$ by:

[TABLE]

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:

[TABLE]

[TABLE]

The average grain size is then given by:

[TABLE]

This approach employs an average size $\bar{a}$ 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* $\rho_{\rm dust}(R,z)$ by the mean mass of one dust particle $(\bar{m}_{\rm gr}=\rho_{\rm gr}4\pi/3\times\bar{a}^{3}),$ we can obtain the number of grains per volume $n_{\rm gr}$. Defining the mean geometrical cross section as $\langle\sigma\rangle=\pi\bar{a}^{2}$, the total grain geometrical cross section per volume at an $(R,z)$ location in the disk is given by:

[TABLE]

where

[TABLE]

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 $\kappa_{\nu}(a_{i})$ for every bin size $a_{i}$. 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 $a_{i}$ and $a_{i}+\Delta a_{i}$. 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:

[TABLE]

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 $T_{\rm dust}$, 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 $T_{\rm gas}$ 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:

[TABLE]

where $S_{\rm d}$ is the dust geometrical cross section per baryon, that is:

[TABLE]

The direct proportionality to the local dust-to-gas ratio $\delta_{\rm dg}=\rho_{\rm dust}/\rho_{\rm gas}$ shows that the total dust cross section is higher for a larger amount of dust, whereas the inverse dependence on $\bar{a}$ 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 $S_{\rm d}=\delta_{\rm dg}\times 6.09\times 10^{-20}\,{\rm cm}^{2}$ (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):

[TABLE]

where $n_{\rm H}$ is the abundance of H atoms in the gas phase, $\epsilon_{{\rm H}_{2}}$ is the H2 recombination efficiency, and $S_{\rm H}(T_{\rm dust})$ 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:

[TABLE]

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:

[TABLE]

where $E_{\rm b}(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 $\sim 855\,$K (Sandford & Allamandola, 1993; Collings et al., 2003; Öberg et al., 2005; Bisschop et al., 2006), but it can be as high as $\sim 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:

[TABLE]

where the number of binding sites per unit grain surface $N_{\rm ss}=8\times 10^{14}$ cm*-2* (e.g. Hollenbach et al., 2009) and $n_{\rm ice}$ is the sum of number densities of all individual ice species.

Photodesorption is modelled following again Visser et al. (2011, see their eq. 6):

[TABLE]

where $F_{\rm UV}$ is the local value of UV flux (obtained from the radiative transfer module, see Sec. 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, 2009b). 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 $10^{4}\,$cm*-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:

[TABLE]

where again the reactions are allowed only in the top monolayer and $n_{\rm hydro}$ 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\,\mu$m and a grain number density $n_{\rm gr}=10^{-12}n_{\rm H}$. Instead, in this paper $n_{\rm gr}$ and $\bar{a}$ 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\,\AA$ and $1\,\mu$m for small grains and ranging between $50\,\AA$ 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 $\chi h$, where $\chi<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:

[TABLE]

[TABLE]

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 Sec. 5. In particular, we use $R_{\rm c}=125\,$AU, $\gamma=0.9$, $\Sigma_{\rm c}=7.03\,$g cm*-2*, $h_{\rm c}=0.11$, and $\psi=0.25$. These parameters yield a disk gas mass $M_{\rm disk}\sim 7\times 10^{-2}\,M_{\odot}$. The gas density structure used in all models is shown in Fig. 2. The total dust-to-gas mass ratio $\Delta_{\rm dg}$ is set to $0.01$. We use a stellar mass of $2.47\,M_{\odot}$, with a bolometric luminosity $L_{*}=38\,L_{\odot}$, and an effective temperature of $10^{4}\,$K (Tilling et al., 2012). The X-ray luminosity is set to $L_{X}=6\times 10^{29}\,$erg s*-1* (Günther & Schmitt, 2009). Finally, we use an accretion rate of $4.5\times 10^{-7}\,M_{\odot}$ 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 $R_{\rm in}$ and $25$ AU. We set an inner disk radius $R_{\rm in}=0.42\,$AU, since the dust sublimation radius is $\sim 0.07\sqrt{L_{*}/L_{\odot}}$, if a dust sublimation temperature of $\sim 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 [math] 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\,\mu$m in the disk. For model STN, we consider $f=0.85$ and a settling parameter $\chi=0.2$ (see Section 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 Sec 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 Section 6.3.

For the models including dust evolution, the three new parameters defining the dust distribution are $\alpha$, $v_{\rm frag}$ , and $\rho_{\rm gr}$. As noted above, $v_{\rm frag}=10\,$m s*-1* and $\rho_{\rm gr}=2.5\,$g cm*-3*. Here we focus on exploring the dependence of the dust and gas properties on turbulence, by assigning $\alpha$ 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 ($\alpha\lesssim 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\times 10^{7}$ photon packages for both the star and the environment radiation, which is assumed to be $1\,G_{0}$, where $G_{0}\sim 2.7\times 10^{-3}\,$erg s*-1* cm*-2* is the UV interstellar radiation field between $911\,\AA$ and $2067\,\AA$ (Draine, 1978). The number of photon packages used in stage 2 is $3\times 10^{6}$ 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 $[{\rm C}]/[{\rm H}]=1.35\times 10^{-4}$ and $[{\rm O}]/[{\rm H}]=2.88\times 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 $\zeta_{\rm CR}=5\times 10^{-17}\,$s*-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 $\sim 2\times 10^{5}\,$yr.

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^{\circ}$, with a beam convolution of $0.52\arcsec\times 0.38\arcsec$, where the position angle of the beam is $82^{\circ}$ (de Gregorio-Monsalvo et al., 2013). The position angle of the major axis of the disk is taken as $137^{\circ}$.

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 $\alpha$. 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 $R_{\rm frag}$ for the three different cases, where $R_{\rm frag}$ varies between $\sim 10$ – $500$ AU with $\alpha$ ranging between $10^{-4}$ and $10^{-2}$ (in fact $R_{\rm frag}$ 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 $\lesssim 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 $\sim 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 $\alpha=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 $\bar{a}$ and the total dust surface area per volume $n_{\rm gr}\langle\sigma\rangle$ are two other important quantities that are directly related to the dust density structure (see Sec. 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 $\alpha=10^{-4}$ case, $\bar{a}$ can be as high as $>10\,\mu$m in the disk mid-plane. However, in this case $\bar{a}$ 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 Eq. 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 $\alpha$ varying between $10^{-4}$ and $10^{-2}$. Interestingly, the total dust surface area that the STN model assumes is quite similar to the $\alpha=10^{-2}$ case, indicating that the parameters chosen for the STN model, with a settling parameter $\chi=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 Sec. 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 Section 4.1) and are generally more abundant than in the STN model, or in the $\alpha=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 $\alpha$ 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 $\sim 100\,$K in the $\alpha=10^{-4}$ case, whereas they reach at most $\sim 50\,$K in the $\alpha=10^{-2}$ case.

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\gtrsim 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 $\gtrsim 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 $k_{\rm ads}({\rm CO})=k_{\rm thdes}({\rm 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 $\alpha=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. 19, it emerges that the region where the gas temperature is lower than the dust temperatures occurs just below the $\tau_{\rm 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 $T_{\rm gas}/T_{\rm dust}<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 $\sim 550\,$AU) shown in Fig. 21. In these regions of the disk, while the mid-plane temperature is the same for all the models, $T_{\rm dust}$ 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.

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. 19), 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 ($\gtrsim 10^{-4}$) for $\alpha\leq 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 $\sim 17-20\,$K in the disk mid-plane, at $R\sim 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 $\sim 300\,$AU. The amount of CO ice in the STN and $\alpha=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)$\gg 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 $T_{\rm dust}$ 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\,\mu$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\,\mu$m. Its mm emission is thus the most compact in radial extent. The radial extent of the STN, $\alpha=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 $\alpha=10^{-4}$ case in Fig. 11, but also present in the $\alpha=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\,\mu$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.

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 $\alpha=10^{-2}$ models are very similar as usual, the $\alpha=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 Sec. 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 C – CO transition. The lower temperature implies that the $J$=3 state is less populated, leading to lower emission in the $J$=3-2 line. The $\tau=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 $\sim 30$, whereas in 13CO and C18O this is of the order of $\lesssim 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 ($R_{\rm CO}^{\rm out}$), since this would be observationally determined by the sensitivity cut. In particular, from the 12CO line, the outer radius would vary between $\sim 350-550\,$AU for $\alpha$ ranging between $10^{-2}$ and $10^{-4}$, whereas the same estimate based on the 13CO line would yield disk outer radii varying between $\sim 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.

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 $\gtrsim$ 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 $R_{\rm CO}^{\rm out}$ and $R_{\rm mm}^{\rm out}$ (taken as the radius where the intensity profiles reach a dynamic range of $30$ at $850\,\mu$m), depends strongly on turbulence. In particular, for the turbulent models shown in Figs. 11-12, the ratio $R_{\rm CO}^{\rm out}/R_{\rm mm}^{\rm out}$ scales from $\sim 1.4$ for $\alpha=10^{-4}$ to $\sim 4$ for $\alpha=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\,\mu$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 $\alpha=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 Sec. 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.

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\,$mm) 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\pm 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 ($\alpha=10^{-3}$). However, all those models assume a CO binding energy of $855\,$K, typical for pure CO ice (see Sec. 2.6.3 for a more detailed discussion). By running a turbulent model ($\alpha=10^{-3}$) with a CO binding energy more appropriate for ice mixtures ($1100\,$K), we obtain a location of the snowline of $\sim 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 $\gtrsim 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.

Figure 16 shows the observed and modelled peak-normalised intensity profiles of both continuum at $850\,\mu$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\arcsec\times 0.55\arcsec$ and $0.72\arcsec\times 0.57\arcsec$, respectively, with a PA$=72^{\circ}$. The vertical scale is used to mimic a peak value of $\sim 300\sigma$. 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 $R_{\rm frag}$ shifting towards the inner regions of the disk). The 12CO data are well recovered by many models, in particular with $\alpha=10^{-2}-10^{-3}$. When $\alpha=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 $R_{\rm CO}^{\rm out}/R_{\rm mm}^{\rm out}$ is directly related to the amount of turbulence in the disk, with the ratio increasing with $\alpha$ (Fig. 17). The exact values of the ratio will depend on the dynamic range of the CO emission used to determine $R_{\rm CO}^{\rm out}$. The value of the ratio will also depend on the wavelength considered to extract $R_{\rm mm}^{\rm out}$, 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.

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 HD163296 by Yen et al., 2016a), is indeed coupled to the ice composition and abundances on the mantles of dust grains.

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 ($\bar{a}\lesssim 0.1\,\mu$m). In the comparison, the gas thermal structure computed by the model is used.

The growth timescale $t_{\rm growth}$ of a dust particle of size $a$ due to turbulent motions can be written as $t_{\rm growth}=a/\dot{a}$, where $\dot{a}$ is (Brauer et al., 2008):

[TABLE]

The relative velocities due to turbulence can be expressed as (Ormel & Cuzzi, 2007):

[TABLE]

where the local Stokes number can be written as:

[TABLE]

We thus obtain that:

[TABLE]

The quantity $c_{\rm s}$ is obtained from the STN-SM model, where $c_{\rm s}=\sqrt{P/\rho_{\rm gas}}$, and the contribution of chemical species to both $P$ and $\rho_{\rm gas}$ is taken into account.

Similarly, the freeze-out timescale can be expressed as $t_{\rm ads}\sim 1/k_{\rm ads}$ (see Eq. 21), where $k_{\rm ads}$ 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:

[TABLE]

For the most abundant 12CO, $M({\rm CO})=28$.

Figure 18 shows the ratio of the two timescales for particles of size $a=0.1\,\mu$m, an assumed turbulent $\alpha=10^{-3}$, and $\rho_{\rm gr}=2.5\,$g cm*-3*. All the other quantities are taken from the STN-SM model. In the whole disk, $t_{\rm growth}$ is larger than $t_{\rm ads}$, and in particular in the mid-plane of the inner disk, $t_{\rm growth}\gg t_{\rm ads}$. 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:

Differences in gas and dust radial extents of protoplanetary disks are readily found in our models and can be largely explained by the difference in optical depth of gas versus dust, without the need to invoke radial drift (the STN model leads to a similar size in continuum as the turbulent models). 2. 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 $\sim$ three for a turbulent $\alpha$ ranging between $10^{-4}$ and $10^{-2}$, with the ratio of the CO and mm radius $R_{\rm CO}^{\rm out}/R_{\rm mm}^{\rm out}$ increasing with turbulence (Fig. 17). 3. 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. 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. 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. 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 $\alpha$, due to a steeper gas temperature gradient in the intermediate disk layers. 7. 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. 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 $\alpha$. 9. 9.

The CO snowline location of HD 163296 at $\sim 90\,$AU, as determined by Qi et al. (2015), can be retrieved by using a binding energy for CO typical of ice mixtures ($E_{\rm b}\sim 1100\,$K). Lower binding energies, in particular for pure CO ice ($E_{\rm b}\sim 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.

Acknowledgements.

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).

Appendix A Disk thermal structure

Figure 19 shows the attenuated UV field in STN and turbulent models. Figure 20 portrays the dust and gas thermal structure of the same models. Finally, Fig. 21 shows the gas and dust temperatures, and the abundances of CO, C, and C+, in a vertical cut at $550\,$AU.

Appendix B Line profiles of CO isotopologues

The spatially unresolved line profiles of three CO isotopologues discussed and shown in Section 4.3.2 already contain important information. The profiles of the $J$=3-2 transition are reported in Fig. 22. 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.

Appendix C High resolution synthetic images

All the synthetic intensity profiles shown in the paper use a convolution beam of $0.52\arcsec\times 0.38\arcsec$ (Section 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\arcsec$ resolution) in Figs. 23- 24. 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 $\alpha\leq 10^{-3}$ show a rather sharp edge in the $850\,\mu$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. 24), where part of the line flux is lost due to high optical thickness of the continuum. This is apparent in the $\alpha=10^{-4}$ case, where the inner disk is heavily populated by large grains, due to low turbulent velocities (see Fig. 3).

In Fig. 25 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 $\alpha$, 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.