© ESO 2018

1 Introduction

Protoplanetary disk evolution adds to the construction of the host young star and simultaneously provides the material for the formation of planets and large rocky bodies. Mass accretion onto the central star and planet formation are therefore connected processes where the gaseous and solid disk components play important and complementary roles (see Armitage 2015, for a review). It is generally assumed that gas is the main disk component accounting for 99% of its mass, while only 1% is in the form of dust. Furthermore, it is also the gaseous component that drives disk evolution and thus the disk structure at different stages. What the relative importance of different mechanisms transporting the material onto the star and dispersing the disk within few Myr is, remains an open question. Disk evolution leads to mass accretion onto the central star (Hartmann et al. 2016), but part of the disk content can be dissipated by high-energy radiation driven winds or magnetic torques (Alexander et al. 2014; Gorti et al. 2016; Bai 2016) and external processes such as encounters (e.g., Clarke & Pringle 1993; Pfalzner et al. 2005) and external photoevaporation (e.g., Clarke 2007; Anderson et al. 2013; Facchini et al. 2016). Depending on which process is dominant, the disk structure will be considerably different. In particular this will impact the distribution of the gas as a function of the distance from the star, which is described by the surface density distribution, Σ_gas . Morbidelli & Raymond (2016) discuss the difference in the resulting disk structures between disk wind models and the standard viscously evolving α-disk model. Suzuki et al. (2016) find that the inner disk can be heavily depleted as winds efficiently remove angular momentum from the disk andcause large accretion onto the star. On the other hand, a different disk wind model by Bai (2016) predicts a surface density profile that scales inversely with the radius in the innermost 10 au, similarly to what is expected in the α-disk case. Reliable observational measurements of Σ_gas as function of radius would therefore be key to understand disk evolution and the relative importance of different processes, as well as how planet formation occurs (Morbidelli & Raymond 2016).

Successful attempts have been made to constrain the disk dust surface density Σ_dust from high S/N spatially resolved millimeter-continuum observations. The thermal emission of millimeter-sized grains is largely optically thin in protoplanetary disks and allows tracing directly the column density of the bulk of the dust (e.g., Andrews 2015; Tazzari et al. 2017). Dust and gas surface densities, however, are not necessarily the same. Even though theyare strongly coupled, different mechanisms are thought to affect their evolution, such as radial drift of millimeter-sized grains modifying the dust surface density on short timescales. This is the case, for example, in HD 163296 where the carbon monoxide (CO) emission extends much further than the continuum dust emission (de Gregorio-Monsalvo et al. 2013; Isella et al. 2016), although the much higher optical depth of the CO lines versus millimeter continuum also contributes significantly to this difference in appearance (Dutrey et al. 1998; Facchini et al. 2017). Another example comes from the Herbig disk HD 169142, where ALMA observations show that the dust appears to be concentrated in two rings between 20–35 au and 56–83 au, whereas gas seen in CO isotopologues is still present inside the dust cavity and gap (Fedele et al. 2017). More extreme examples are transitional disks with large cavities severely depleted in dust content, but less in gas, with gas cavities smaller than dust cavities (van der Marel et al. 2013, 2016a; Bruderer et al. 2014; Zhang et al. 2014). Taken together, it is clear that gas and millimeter dust emission do not follow each other and that a more direct tracer of the gas surface density profile is needed.

The main gaseous component is molecular hydrogen (H₂), whose emission lines at near- and mid-infrared wavelengths are very weak and are superposed on strong continuum emission (e.g., Thi et al. 2001; Pascucci et al. 2006). In contrast, CO pure rotational transitions at millimeter wavelengths are readily detected (e.g., Dutrey et al. 1996; Thi et al. 2001; Dent et al. 2005; Panić et al. 2008). CO is the second-most abundant molecule, with a chemistry that is very well studied and implemented in physical-chemical codes. Owing to their optically thin emission lines, less abundant CO isotopologues are considered better tracers of the gas column than ¹² CO (Miotello et al. 2014, 2016; Schwarz et al. 2016). This paper investigates whether spatially resolved observations of such lines may be good probes of the gas surface density distribution, once freeze-out and isotope-selective photodissociation are taken into account.

Thanks to the unique sensitivity of ALMA, large surveys of disks in different star-forming regions are being carried out to trace gas and dust simultaneously (Ansdell et al. 2016; Barenfeld et al. 2016; Eisner et al. 2016; Pascucci et al. 2016). In particular, CO isotopologues pure rotational lines have been targeted but lines are weaker than expected and at the short integration time of typically 1 min, the line data have low S/N. This may be interpreted as lack of gas due to fast disk dispersal, or as lack of volatile carbon that leads to faint CO lines (Favre et al. 2013; Ansdell et al. 2016; Schwarz et al. 2016; Miotello et al. 2017). As a consequence, ¹³ CO lines are optically thin toward a larger part of the disk and become much better tracers of the gas columns than C¹⁸ O lines, which were thought to be preferred pre-ALMA, but are often undetected. Furthermore, ¹³ CO is less affected by isotope-selective photodissociation than C¹⁸O (Miotello et al. 2014). The focus of this paper is thus on ¹³CO, rather than on other isotopologues. Recently, Zhang et al. (2017) have proposed the optically thin lines of ¹³ C¹⁸O, successfully detected in the TW Hya disk, to be a good probe of Σ_gas (R) at small radii. We briefly investigate this possibility from our modeling perspective using another low abundance isotopologue C¹⁷ O.

Williams & McPartland (2016) have presented a method to extract disk gas surface density profiles from ¹³ CO ALMA observations. Assuming a self-similar density distribution as given by viscous evolution theory, the authors create a large gallery of parametrized models, which they use to extract the disk parameters with a MCMC analysis. Although we share the same scientific question, this paper follows a different approach, i.e., we fit a power-law profile to our simulated ¹³ CO radial profiles as if they were observations, in order to recover the power-law index γ. Our analysis shows that it may be difficult to infer both γ and the critical radius independently for small disks.

In Sect. 2, the physical-chemical modeling is presented. The results of the modeling investigation are presented in Sect. 3 and their implications are discussed in Sect. 4.

2 Modeling

2.1 DALI

The physical-chemical code DALI is used for the disk modeling (Dust And LInes; Bruderer et al. 2012; Bruderer 2013). DALI calculates the dust temperature, T_dust, and local continuum radiation field from UV to mm wavelengths with a 2D Monte Carlo method, given an input disk physical structure and a stellarspectrum as source of irradiation. Then, a time-dependent (1 Myr) chemical network simulation is run. Subsequently, the gas temperature, T_gas, is calculated through an iterative balance between heating and cooling processes until a self-consistent solution is found and the non-LTE excitation of the molecules is computed. In this paper, the final outputs are spectral image cubes of CO isotopologues computed by a ray-tracer module. Throughout this paper, the J = 3–2 line is produced and analyzed, but results are similar for J = 2–1. Line intensities are provided in K km s⁻¹, rather than Jy beam⁻¹ km s⁻¹, in order to not depend on the beam size. The line intensity can be converted to Jy beam⁻¹ km s⁻¹ using the relation in Appendix B of Bruderer (2013). For the ray-tracing, a distance of 150 pc is assumed.

DALI has been extensively tested and benchmarked against observations (see e.g., Fedele et al. 2013) and a detailed description of the code is presented by Bruderer et al. (2012) and Bruderer (2013). For this work, a complete treatment of CO freeze-out and isotope-selective processes is included (for more details, see Miotello et al. 2014, 2016). The volatile carbon abundance is assumed to be high, [C]/[H] = 1.35 × 10⁻⁴, and constant across the disk.

2.2 Grid of models

The disk surface density distribution is often parametrized by an exponentially tapered power-law function, following the prescription proposed by Andrews et al. (2011). Physically, this represents a viscously evolving disk, where the viscosity is expressed by ν ∝ R^γ (Lynden-Bell & Pringle 1974; Hartmann et al. 1998). The self-similar surface profile is expressed by $${\text{Σ}}_{\text{gas}}={\text{Σ}}_c{\left(\frac{R}{R_c}\right)}^{-\text{γ}}\exp \left[-{\left(\frac{R}{R_c}\right)}^{2-\text{γ}}\right],$$(1)

where R_c is the critical radius and Σ_c is the surface density at the critical radius. This parametrization has often been employed to model disks with DALI (see e.g., Miotello et al. 2014, 2016). Adopting an exponential taper to the power-law profile of the surface density distribution (see Eq. (1)), as suggested by viscous evolution theory, has the inconvenience that the profile slope depends on two free parameters, R_c and γ, which can be difficult to disentangle with the current low S/N of the data in the outer regions of disks. A simplification comes from the assumption that Σ_gas has a pure power-law dependence with radius. In this case the power-law index γ is left as a single free parameter: $${\text{Σ}}_{\text{gas}}=\{\begin{array}{ll}{\text{Σ}}_c{\left(\frac{R}{R_c}\right)}^{-\text{γ}},\hfill & R\le R_{\text{out}},\hfill \\ 0,\hfill & R>R_{\text{out}}.\hfill \end{array}$$(2)

In this work both parametrizations have been employed to design the input density structures.

First, some of the model results presented in Miotello et al. (2014) using the self similar profile have been analyzed. The free parameters are then R_c, γ, and M_disk, whose values are reported in Col. 2 of Table 1. A value for R_out is also reported, but this is simply needed for the simulation and does not have any effect on the disk structure, as the exponential cutoff truncates the disk at smaller radii. Parameters defining the disk vertical structure are reported in Table 1 of Miotello et al. (2014) and simulate a medium-flared disk with large dust grains settled. From here on, such disk models where the input surface density distribution is set by Eq. (1) will be called self-similar disk models.

Additional models have then been run with the simple power-law surface density structure (see Eq. (2)). In this case the free parameters are R_out, γ, and M_disk as shown in Col. 3 of Table 1. From here on, disk models where the input surface density distribution is set by Eq. (2) will be called simple power-law disk models.

In orderto investigate different disk mass regimes and then compare with recent observations of CO isotopologues in protoplanetary disks (Ansdell et al. 2016, 2018), a range of models from less massive to high mass disks have been run (see Table 1). These observations also show moderately extended emission, therefore the outer radius in the models has been set to 200 au. The raytracing of all models presented here is carried out assuming the disk to be at a distance of 150 pc, representative of that of the nearby star-forming regions. Furthermore, the output synthetic images have either been left unconvolved or convolved with a moderate resolution beam of 0.2′′ . In the first case, the resolution is that of the radial grid assumed for the simulation. The calculation is carried out on 75 cells in the radial direction (logarithmically spaced up to 30 au and then linearly spaced for a radial resolution of ~ 10 au) and 60 in the vertical direction models (Miotello et al. 2014). In the new models with the simple power-law surface density structure, a denser grid with 95 cells in the radial direction (improved radial resolution in the outer disk of ~ 3 au) and 80 in the vertical direction was chosen.

All models presented in this work are T Tauri-like disk models. The stellar spectrum is assumed to be a blackbody at a temperature T_eff = 4000 K with excess UV due to accretion and the stellar luminosity is L_bol = 1 L_⊙ (see Miotello et al. 2014, 2016). In order to explore the effect of different stellar properties on the results of this work, we have compared the self similar disk models with some of the Herbig-like disk models presented in Miotello et al. (2016), where T_eff = 10 000 K and L_bol = 10 L_⊙. We find that stellar temperature and luminosity do not quantitatively affect the results found with the T-Tauri like disk models, except for pushing the preferred fitting region further out in the disk for massive disk models (see Sect. 3.4). Such massive Herbig disk models (M_disk = 10⁻³ M_⊙) show radial intensity profiles which are a factor of 1.6 brighter than those found for the respective T-Tauri models.

3 Results

The first test is to compare the input disk surface density profiles for the DALI models with the simulated ¹³ CO profiles, as if they were observed at infinite resolution, i.e., not convolved with any beam. In practice, the image resolution is the physical grid resolution in our models (see Sect. 2.2). The next step is then to convolve the simulated ¹³ CO profiles with a typical ALMA observation beam. The results obtained in the two cases, where Σ_gas is parametrized by a simple power-law or following the self-similar solution are presented below.

3.1 Simple power-law

The input disk surface density profiles, assumed to be a power-law function of the radius (see Eq. (2)), are compared with the simulated ¹³ CO profiles for three representative models (M_disk = 10⁻⁴ M_⊙, R_out = 200 au, γ = 0.8, 1, 1.5) and are shown in panel a of Fig. 1. The surface density profiles (in g cm⁻²) have been rescaled by a constant factor in order to visually match the line intensity profiles. The disk can be divided into three regions: an inner part in which the ¹³CO radial profile underestimates Σ_gas, a central zone where the two coincide, and an outer region where the line emission drops and deviates from the surface density distribution (see case with γ = 1.5, for clarity). At small radii the divergence iscaused by the fact that ¹³CO emission becomes optically thick as the columns are very high. At large radii the ¹³ CO column densities become too low, lower than 10¹⁵ cm⁻², and ¹³ CO is not able to efficiently self-shield (van Dishoeck & Black 1988; Visser et al. 2009).

There are two ways to further investigate these cases. The column densities of gas-phase ¹³ CO calculated from the surface to the midplane are shown as function of the disk radius in panel c of Fig. 1 in blue, light blue, and purple for models with γ = 0.8, 1, 1.5, respectively.The column density can be compared with the line luminosity profiles presented in panel a of Fig. 1. The radius where the line luminosity profiles decrease and deviate from the surface density distribution is similar to that at which the ¹³ CO column densities drop below 10¹⁵ cm⁻². This illustrates that inefficient CO self-shielding affects CO isotopologues intensity profiles and needs to be considered in the disk outer regions.

As reported in Table 1, less and more massive models have been run for the simple power-law case. We observe similar trends, but the radial location of the three regions described above are radially shifted to smaller or larger radii, depending on disk mass (see Fig. A.2). Compared with the representative model with mass M_disk = 10⁻⁴ M_⊙, for M_disk = 10⁻⁵ M_⊙ the simulated radial intensity profiles better follow Σ_gas at smaller radii as the line is less optically thick. On the other hand, ¹³ CO column densities smaller than N = 10¹⁵cm⁻² are reached at smaller radii. The opposite happens for M_disk = 10⁻² M_⊙, where the line emission is optically thick up to almost 100 au, but CO self-shielding becomes inefficient much further out, around 200 au.

Fig. 1 Panel a: 13CO line intensity radial profiles (solid lines) obtained with three representative disk models with input surface density distribution Σgas (dashed lines) chosen to be a simple power-law (see Eq. (2)). The model parameters are Rout = 200 au, Mdisk = 10−4 M⊙, and γ = 0.8, 1, 1.5 shown in dark blue, light blue, and purple, respectively (top, middle, and bottom panels). No beam convolution is applied to the simulated images and a distance of 150 pc is assumed. Panel b: 13 CO line intensity radial profiles compared with input surface density profiles as in panel a, but convolved with a 0.2′′ beam. Panel c: column densities of gas-phase 13CO calculated from the surface to the midplane shown as function of the disk radius for the three representative models with simple power-law surface density. The dotted black line indicates the column density at which CO self-shielding becomes inefficient (N = 1015 cm−2). Ice-phase CO column densities are not shown in the plot as they are below the values shown here.

3.1.1 Fitting of the surface density power-law index γ

We now address the question of whether it is possible to retrieve the underlying surface density power-law index γ by fitting the line intensity profiles. We refer to the power-law index used as input in DALI as γ_input and we label the fitted value γ_fit.

The fit of a power-law profile to the simulated intensity profiles is performed as a linear fit in the log I- logR space. As it is clear from Fig. 1 that ¹³CO radial intensity profiles follow Σ_gas(R) only in a particular region, the fit is not carried out over the whole extent of the disk, but over a radial range that spans from R_start to R_cut. The starting point R_start is kept fixed and chosen to be just outside the inner region where the emission is optically thin, typically 30–40 au. The range of radii over which the fitting is performed is then varied by changing the value of R_cut.

The results of the fitting are presented in Fig. 2 for the single power-law case. The dashed lines illustrate the value of γ_input, while the colored dots show the value of γ_fit if different R_cut are chosen for the fitting. The starting radius for the fitting procedure R_start is indicated by the dotted vertical line.

For low mass disk models (M_disk = 10⁻⁴ M_⊙, panel a of Fig. 2), the power-law index γ_input can be retrieved within 20% of uncertainty in a range of radii that goes from 30–40 au to ~100 au, the typical scales probed with ALMA. Furthermore, the retrieved γ_fit usually overestimates the input value, except for the γ = 1.5 case. Fixing R_start to either 30or 40 au does not change qualitatively the outcome of the fit. For R_cut larger than 100 au, the fitted values deviate significantly from γ_input as one enters the region where self-shielding becomes inefficient and the intensity profiles deviate from Σ_gas. For higher mass disk models (M_disk = 10⁻² M_⊙, panel b of Fig. 2), one needs to chose a larger R_start of 100 au, as ¹³CO is optically thick much further out than for the low mass disk case. Furthermore, the uncertainties on γ_fit are larger.

In order to choose the correct radial range over which to perform the fit, one would need to know in which mass regime the observed disk is. More precisely, one would need to constrain the total ¹³ CO gas mass. For faint observed ¹³CO fluxes, Miotello et al. (2016) have found a linear relation between line emission and total disk mass. This could be used to constrain M_disk. For brighter observed ¹³CO fluxes, this relation is less reliable and one could employ the dust masses obtained from continuum emission, multiplied by the canonical factor of 100. This would provide an upper limit to the disk mass (see Sect. 4 for discussion on carbon depletion). Another possibility is to combine ¹³ CO and C¹⁸ O line fluxes (Miotello et al. 2016), if the latter are available from the same observation.

Fig. 2 Results of the power-law fitting of the simulated 13CO line intensity profiles obtained with the simple power-law models (Mdisk = 10−4 M⊙ and Mdisk = 10−2 M⊙ in panels a and b, respectively). The values of the fitted power-law index γfit are shown by the filled dots as a function of the radial range over which the fitting was carried out for the high-resolution images (0.01′′ beam). The empty symbols in panel a show the fitted power-law index γfit when the simulated images are convolved with a 0.2′′ beam. The model input power-law indexes are shown by the dashed lines. The dotted line shows the starting point of the fitting Rstart.

3.1.2 Convolution with a beam

Finally, the simulated line intensity radial profiles need to be convolved with a synthesized beam to check how the trends are affected by this operation. Observations taken from recent ALMA disk surveys have moderate angular resolution, between 0.2′′ and 0.3′′, corresponding to 30–40 au diameter at 150 pc (Ansdell et al. 2016). To mimic such observations, we convolve our simulated profiles with a 0.2′′ beam. The convolved ¹³CO line intensity radial profiles are shown in panel b of Fig. 1.

As in Sect. 3.1.1, a power-law fitting of the convolved radial intensity profiles has been carried out. The starting point R_start has been fixed to 30 au and the results are shown in panel a of Fig. 2. The values of γ_fit do not qualitatively differ from those found with the unconvolved radial profiles, shown in panel a of Fig. 2. Even with a beam convolution of 0.2′′ the power-law index can be retrieved to an accuracy of ~20%, if the fit is carried out in a region that goes from ~30 au to ~100 au.

3.2 Self-similar disk models

Similar trends as those presented in Sect. 3.1 have been found for the self-similar case. The input disk surface density profiles compared with the simulated ¹³ CO profiles for three representative self-similar models (M_disk = 10⁻³ M_⊙, R_c = 200 au, γ = 0.8, 1, 1.5) are presented in panel a of Fig. 3. The surface density profiles have been rescaled by a constant factor in the plots to visually match the line intensity profiles. As for the simple power-law case it is possible to identify three regions: the inner disk in which ¹³ CO lines are optically thick and the emission profile is flatter than the Σ_gas profile, a central region in which the two profiles match, and the outer disk where the line intensity underestimates the surface density distribution.

In the self-similar models, a secondary effect adds to inefficient CO self-shielding to reduce the emission in the outer disk. As shown in Fig. A.1, the simple power-law representative model (panel b) for a less massive disk (M_disk = 10⁻⁴ M_⊙) is warmer than the self-similar representative disk model (panel a, M_disk = 10⁻³M_⊙). In this second case, the dust temperature T_dust drops below20 K in the midplane at radii larger than 50 au and a significant fraction of ¹³ CO is frozen out, thus not contributing to the line emission (see panel a of Fig. A.1). Moreover, the inner region where the emission lines are optically thick is smaller compared to the power-law case, and the radial range over which the line intensity profile follows Σ_gas is shifted toward smaller radii.

The column density, reported in panel b of Fig. 3, can be compared with the line intensity profiles presented in panel a. The radius where the line intensity profiles decrease and deviate from the surface density distribution is similar to that at which ¹³ CO column densities drop below 10¹⁵ cm⁻² (at ~ 300 au for γ = 0.8, ~ 250 au for γ = 1.0, and ~ 200 au for γ = 1.5).

Similar to the simple power-law case, simulated ¹³CO radial intensity profiles have been fitted with an exponentially tapered power-law (see Eq. (1)) to try and retrieve the power-law index γ_input and the critical radius R_c,input. This procedure is applied to models where M_disk = 10⁻⁴ M_⊙, in order to directly compare with the results found in the simple power-law case, instead of M_disk = 10⁻³ M_⊙, as discussed above. The intensity profiles for the lower mass (M_disk = 10⁻⁴ M_⊙) self-similar disk models are reported in Fig. A.3. Fit parameters are much more ambiguous when fitting the intensity profiles obtained with the self-similar disk models. As shown in Fig. 4, the retrieved γ_fit and R_c,fit can significantly deviate from the original γ_input. If R_c,input is 200 au, it is still possible to distinguish γ_fit between models with different γ_input (panel a of Fig. 4). This holds in particular if R_cut > 100 au, in which case γ_input is retrieved within 20% and R_c,input within 12%, consistent with Williams & McPartland (2016) large HD163296 disk. If instead R_c is smaller, i.e., 30 au or 60 au (panels b and c of Fig. 4, respectively), this is not true as γ_input cannot be reliably retrieved. It is still possible to obtain a good estimate for R_c, which can be retrieved within 13% and 30% if R_c,input = 60 and 30 au, respectively. The complication is given by the fact that for smaller values of R_c the line intensity radial profile follows the Σ_gas profile over a very limited radial range and the fitting does not perform well enough. The effects of either optical thickness or inefficient self-shielding dominate throughout most of the disk’s extent and the radial line intensity profile does not resemble an exponentially tapered power-law. In summary, the obtained γ_fit has a poor relation with γ_input if R_c,input is small (i.e., 30 and 60 au). On the other end, it is always possible to retrieve R_c with moderate uncertainty.

Recent ALMA disk surveys (e.g., Ansdell et al. 2016; Pascucci et al. 2016) have shown that typical protoplanetary disks are usually fainter and less extended than previously thought. Therefore, the typical sensitivity of these observations does not allow to detect the regions where the exponential taper would dominate, if Σ_gas were described by Eq. (1). Moreover, our results show that in such external regions the inefficiency of CO self-shielding prevents CO isotopologues to be used to constrain the surface density distribution. Therefore, we consider the power-law prescription as a simplification of the self-similar disk model. Both behave in the same way, as a simple power-law disk model, in the region where CO isotopologue can be reliably employed to trace Σ_gas.

Fig. 3 Panel a: 13CO line intensity radial profiles (solid lines) obtained with three representative disk models with input surface density distribution Σgas (dashed lines) given by the self-similar disk model (see Eq. (1)). The model parameters are Rc = 200 au, Mdisk = 10−3 M⊙, and γ = 0.8, 1, 1.5 shown in dark blue, light blue, and purple, respectively (top, middle, and bottom panels). Panel b: column densities of gas-phase (solid lines) and ice-phase (dashed lines) 13 CO calculated from the surface to the midplane shown as function of the disk radius for the three representative self-similar disk models. The dotted black line indicates the column density at which CO self-shielding becomes inefficient (N = 1015 cm−2).

Fig. 4 Results of the power-law fitting of the simulated 13CO line intensityprofiles obtained with the self-similar disk models with Mdisk = 10−4 M⊙. The values of the fitted power-law index γfit are shown by the filled dots as a function of the radial range in which the fitting was carried out. The model input power-law indexes are shown by the dashed lines. The dotted line shows the starting point of the fitting Rstart .

3.3 Inner disk surface density profile from C¹⁷O line intensity radial profiles

Our models show that ¹³CO is not a good tracer of Σ_gas in the inner 30–40 au. Recently, Zhang et al. (2017) have used the rarer ¹³ C¹⁸O isotopologue to try determining the surface density profile in the inner regions of the well-studied TW Hya disk, inside the CO snowline. Another rather low abundance CO isotopologue is C¹⁷ O, expected to be around 24 times more abundant than ¹³ C¹⁸O based on isotope ratios. We thus analyze the intensity profiles obtained for C¹⁷ O to check how they can be used to infer the surface density profile in the disk inner regions.

The simulated C¹⁷O intensity profiles and column density profiles for six power-law disk models (R_out = 200 au, M_disk = 10⁻⁴, 10⁻²M_⊙ and γ = 0.8, 1, 1.5) are presented in Fig. 7. For the 10⁻⁴M_⊙ mass disk models, the intensity profiles follow the shape of Σ_gas in the very inner disk (R < 10 au, see panel a of Fig. 7), apart for the steeper γ = 1.5 model where optical thickness starts to play a role for R < 5 au. On the other hand, C¹⁷O column densities drop below 10¹⁵cm⁻² very early, at R < 10 au. Thus, any information about Σ_gas from C¹⁷O intensity profiles is lost for the outer disk regions.

The picture changes if the results of the more massive (M_disk = 10⁻²M_⊙) disk model are analyzed. C¹⁷ O emission is optically thick out to 50 au, except for the model with γ = 0.5, where Σ_gas is less steepin the inner disk (see panel b in Fig. 7). Therefore, C¹⁷ O is not a good tracer of the surface density distribution in the inner disk in massive disks. An even rarer CO isotopologue, as ¹³ C¹⁸O, should be used in these cases, as was done for TW Hya (Zhang et al. 2017). A possible complication comes from the fact that its faint emission may be shielded by optically thick continuum emission, and this would prevent ¹³ C¹⁸O from being a good tracer of the gas column density.

As for ¹³CO, it is possible to fit the simulated intensity profiles with a power-law and check for which conditions γ_fit resembles γ_input well enough.

3.4 A simple prescription to fit γ: the slope–pivot region

The fitting method of the surface density power-law index presented in Sect. 3.1 is not straightforward to be applied to observations, as this approach involves the choice of two parameters, R_start and R_cut, which would be difficult to determine. On the other hand, there is a relation between the “slope–pivot-point” (see Fig. 5), where ¹³ CO radial intensity profile starts to deviate from the surface density profile due to optical depth, and the column density. The slope–pivot point will always be close to where τ = 1, thus around the same intensity modulo temperature effects, which however are not significant for low-J CO lines. One can then perform the fit of γ between the slope–pivot point and the point where the ¹³CO column density is larger than 10¹⁵ cm⁻². Such radial region will always be around the same narrow range of line intensity.

More specifically, at the radius where ${\tau }_{{13}_{\text{CO}}}\text{=}\text{1}$, the ¹³CO intensity is roughly 6 K km s⁻¹ (see Fig. 5). This results from the line intensity, which is controlled by a combination of the line opacity and the temperature. The intensity at the frequency ν of a gas column with line opacity τ_ν and uniform temperature T is $$I_{\text{ν}}=B_{\text{ν}}(T)\cdot (1-e^{-{\tau }_{\text{ν}}}),$$(3)

where B_ν(T) is the Planck function. Thus, for τ_ν = 1, I_ν depends only on the temperature. In the Rayleigh–Jeans approximation, I_ν ∝ T. Protoplanetary disks have strong gradients in physical conditions both in radial and vertical directions and thus cannot be taken as at uniform temperature. However, the bulk of the disks mass is close to the midplane, where the physical conditions change less than in the warm molecular layer of the upper disk (van Zadelhoff et al. 2001; Bruderer et al. 2012). Mass tracers such as ¹³ CO trace mostly regions close to the midplane. In the midplane, the temperature is a weak function of the radius with T ~r^−q with 0 < q < 1. In the simple power-law case, the pivot point is at radii between 60 and 100 au (see Fig. 5) and the temperature at the radius of the pivot point varies by less than a factor of 2. Thus, the intensity also changes by less than a factor of 2. This explains the similar intensity at the pivot point. The explanation may also apply to other molecules with emission mostly from regions close to the midplane (e.g., ¹²CO or C ¹⁸ O).

We have compared our results in the simple power-law case and self-similar disk case, to test that one can always retrieve γ safely around the same range or line intensity. In Fig. 5 the ¹³CO emission radial profiles are compared with the total line optical depth τ and the ¹³ CO column density. The line emission can be used to trace the gas distribution if τ < 1 (shaded region in panels a and d), i.e., for radii larger than those shown by the dotted lines. On the other hand, the ¹³ CO column density needs to be larger than 10¹⁵ cm⁻² where the self-shielding is effective (shaded region in panels b and e), i.e., for radii smaller than those shown by the dashed lines. Combining these two requirements one is restricted to a radial range, the “pivot region”, where the ¹³ CO intensity is always around 6 K km s⁻¹, for each value of γ and no matter if theunderlying Σ_gas is given by a simple power-law or a self-similar disk model. For all models analyzed in this work, the minimum and maximum flux values relative to the slope–pivot regions are 3 and 8 K km s⁻¹. For the Herbig models with M_disk = 10⁻³ M_⊙ in particular, if one computes the line optical depth and the column density, as shown in Fig. 5, the pivot point shifts toward larger radii. The slope–pivot region is however still located around the value of 6 K km s⁻¹ as for the T Tauri models. This means that, in practice, one would only need to measure the slope of the ¹³ CO radial intensity over a narrow range of radii close to this contour to retrieve the surface density slope (see shaded region in panels c and f). The pivot region will always be very small, as the line intensity that corresponds to N(¹³CO) = 10¹⁵ cm⁻² is very close – within a factor of a few – to that where the line becomes optically thin. The models presented in Fig. 5 have a total disk mass of 10⁻⁴M_⊙ and therefore ¹³CO emission lines are still optically thin in the interested region of the disk. If the disk was more massive, for example, with M_disk = 10⁻²M_⊙, the same argument described above would be valid at larger radii or for less abundant CO isotopologues, for example, for C¹⁷ O. Also for C¹⁷ O, one could retrieve γ by fitting the radial profile at a given contour at 6 K km s⁻¹. In fact, in principle one could put together the surface density profile over the entire radial range by putting together piecemeal the slopes at the pivot points of each of the CO isotopologues, from ¹³ C¹⁸O in the innermost part of the disk to C¹⁷O, C¹⁸ O, and ¹³ CO moving outward in the disk, and finally ¹²CO itself in the outermost region. If the observed disk is inclined by more than about 30°, the data should be deprojected in order to apply this method.

This finding is promising as it relates the observable line intensity radial profile directly with the column density, with no need to know a priori the radial range for the fit of γ or disk mass. On the other hand, the observer would need good angular resolution and sensitivity to perform this fit at a given contour. Typically, one would need a resolution of ~ 0.1′′–0.2′′ to determine the slope of the surface density distribution in disks. The peak flux at the slope–pivot point would lie around 10 mJy beam⁻¹ km s⁻¹ for the ¹³CO J = 3−2 line (and around 3.5 mJy beam⁻¹ km s⁻¹ for the ¹³CO J = 2−1 transition) at a resolution of 0.1′′, and 40 mJy beam⁻¹ km s⁻¹ (and around 14 mJy beam⁻¹ km s⁻¹ for the ¹³CO J = 2−1 transition) at a resolution of 0.2′′ for a disk at 150 pc. Requesting a resolution of 0.1′′ and a S/N of approximatively 5 would make the integration time very long even with ALMA, but increasing the beam size to 0.2′′ results in sufficient S/N in less than 1 h observation, as the data do not need to be spectrally resolved.

Fig. 5 Total 13CO line optical depth τ (panels a and d), 13CO column density (panels b and e), and 13CO intensity radial profiles for simple power-law disk models (top panels) and self-similar disk models (bottom panels). The total disk mass is set to Mdisk = 10−4M⊙ and γ = 1. The shaded regions in panels a and d show the range where τ < 1, which holds for radii larger than those shown by the dotted lines. The shaded region in panels b and e shows the region where 13CO column densities are larger than 1015 cm−2, i.e., for radii smaller than those shown by the dashed lines. The shaded regions in panels c and f show the range of 13 CO line intensityfor which the conditions τ < 1 and N13CO > 1015 cm−2 hold. The fit of γ is performed over the “slope–pivot region” between the dotted and dashed lines, shown by the red arrow in panel c.

3.5 Case study: TW Hya

We test the slope–pivot region fitting on the TW Hya disk observations (i ~6°, d = 60.1 pc as found by Bailer-Jones et al. 2018) published by Schwarz et al. (2016). More precisely, the C¹⁸ O J = 3−2 transition has been considered. According to the model presented by Schwarz et al. (2016), the ¹³ CO J = 3−2 line is optically thick throughout the entire disk, as expected since TW Hya is a much brighter and more massive disk than those considered in our work. Therefore, C¹⁸ O is needed in order to constrain Σ_gas, as this is found to be optically thin everywhere. The observation beam size is ~ 0 . ′′5 × 0 . ′′ 3 and, at such resolution, the slope–pivot point is expected to be at 0.15 Jy beam⁻¹ km s⁻¹.

Figure 6 shows the observed C¹⁸O (J = 3−2) radial profile of TW Hya in black, with the 3Σ level of uncertainty shown by the shaded region. If a power-law is fitted to the observed radial profile between ~ 0.1 and 0.2 Jy beam⁻¹ km s⁻¹, corresponding to radii between ~15 and 30 au (see red line in Fig. 6), the fitted power-law index is γ = 0.85, which is within 15% of the model value γ = 1 assumed by Schwarz et al. (2016) and Kama et al. (2016) for the TW Hya disk. With a S/N > 20 in the slope–pivot region, the error on gamma is negligible, less than ± 0.05. The slight mismatch between the data and the power-law profile on the left side of the slope–pivot region is due to the optical thickness of the emission that becomes more important. On the right side of the slope–pivot region, instead, the mismatch is probably due to the fact that C¹⁸O is being selectively photodissociated. For less massive disks, one should prefer ¹³ CO over C¹⁸ O to trace Σ_gas, as the first is less affected by isotope-selective photodissociation compared with the latter (as shown by Miotello et al. 2014).

High angular resolution and sensitivity are however needed to apply this method, but ALMA can achieve them. In the particular case of TW Hya, where our method could be applied, the rms was 12 mJy beam⁻¹ in a 0.1 km s⁻¹ channel.

Fig. 6 Radial intensity profile of the C18O (3–2) line observed with ALMA (Schwarz et al. 2016) is shown in black and the 3Σ level of uncertainty is presented by the shaded region. The red line shows the power-law function used to fit the data at the slope–pivot region.

4 Discussion

The results presented in Sect. 3 indicate that constraining Σ_gas from spatially resolved observationsof ¹³CO is a challenging task. However, if the underlying surface density distribution is a power-law, it is possible to constrainits steepness by fitting the emission coming from the correct portion of the disk, in particular at large-enough distances from the star that the observed line is optically thin, but at small-enough radii that photodissociation and freeze-out are not too significant. Moreover, in these intermediate regions, the radial gradient of gas temperature is small, and thus does not affect the gradient of the radial intensity profile of the ¹³ CO line. From the models shown in this paper, it is possible to retrieve the dependence of the surface density with radius with good accuracy (20% on the power-law exponent) especially if one has anestimate of the ¹³CO total mass or if one has high-enough resolution data to fit around the pivot-point region at 6 K km s⁻¹.

However, in the approach presented in Sects. 3.1.1 and 3.2, there is the very strong assumption that the surface density is either a power-law function of the radius or it is given by the self-similar model. However, the actual underlying surface density in real protoplanetary disks is of course not known. Moreover, the main difficulty in the fitting is due to the relatively narrow radial range where ¹³ CO traces well the underlying surface density, with this radial range limited by optical depth effects in the inner disk and inefficient self-shielding in the outer regions. Determining additional free parameters, such as the characteristic radius of the self-similar profile, leads to a large degeneracy between the power-law index and the new parameters. To overcome the issue of ¹³ CO being optically thick at small radii, one possibility is to probe Σ_gas in the inner disk by observing very optically thin isotopologues, such as C¹⁷ O and ¹³ C¹⁸O. The latter was successfully detected with ALMA in the closest and best studied protoplanetary disk, TW Hya, by Zhang et al. (2017). The inconveniences of this emission line are its faintness combined with the possibility that the continuum becomes optically thick shielding ¹³ C¹⁸O emission at very small radii. Moreover, the gas temperature structure has a strong impact on the line emission in the inner disk, therefore possibly independent constraints of the thermal structure would be needed.

A much simpler approach is to fit the power-law index γ at the slope–pivot region as explained in Sect. 3.4. This allows a reliable estimate of γ both in the case of a simple power-law and self similar disk, by only fitting the intensity profile at the right intensity contour. The case study presented in Sect. 3.5 shows the applicability of this method to resolved observations of disks with ALMA.

An additional caveat is volatile carbon depletion, a process that may be happening in a large fraction of protoplanetary disks (Favre et al. 2013; Kama et al. 2016; Schwarz et al. 2016; Miotello et al. 2017). In our models, the volatile carbon abundance is assumed to be high and constant throughout the disk. If carbon depletion takes place, but it is constant throughout the disk especially around the pivot region, this should not add major complications to retrieving the surface-density profile. However, the nature of this depletion process is not yet known, but it is possible that carbon is sequestered from CO into CO₂ and more complex ice species in the outer disk, and then drifted inward following the large dust grains. If this is true, one would expect a radially dependent decrease of CO abundances in the outer disk, on top of the reduction due to inefficient self-shielding and freeze-out. As the ices reach the inner disk, the carbon should be liberated into gas phase and quickly return into CO, which would then present an increased emission (Du et al. 2015). This would introduce a new degree of degeneration in the employment of CO isotopologues as tracers of the disk surface density distribution. Interestingly, the expected enhancement of ¹³ C¹⁸O in the inner disk of TW Hya has not been found by Zhang et al. (2017), showing that much is still to be understood about volatile depletion in protoplanetary disks.

5 Summary and conclusion

In this work we have addressed the issue of determining the gas surface density distribution in protoplanetary disks by using resolved high signal-to-noise observationsof ¹³CO. Simulated ¹³CO intensity radial profiles have been produced using the physical-chemical code DALI (Bruderer et al. 2012; Bruderer 2013), with two different input surface density profiles: a simple power-law with radius, and the self-similar solution given by viscously evolving disk theory. By comparing the input Σ_gas profiles with the output intensity profiles, one always finds the following:

• ¹³CO emission follows the slope of Σ_gas(R), but only in a specific disk region. For very small radii the low-J ¹³CO lines become optically thick and underestimate the surface density, while in the outer disk this happens because the ¹³CO column density drops below 10¹⁵ cm⁻² and self-shielding becomes inefficient.

• When fitting a power-law to the simulated intensity profiles, it is possible to recover the input power-law index γ_input only by performing the fit over the appropriate radial range. This holds in particular for simple power-law disk models, where γ_input can be retrieved within 20% uncertainty between 30 and 100 au, even when the profiles are convolved with a 0.2′′ beam.

• In the self-similar case, it is not always possible to reliably retrieve γ_input by fitting a self-similar model to the intensity profile. R_c can instead be always retrieved within 30%.

• Fitting the power-law index γ in a narrow range around the slope–pivot region of the intensity profile allows a reliable estimate of γ in the case of both a simple power-law and self-similar disk. The slope–pivot region is always located around 6 K km s⁻¹. Application of such a method is shown in the case study of the TW Hya disk.

• If carbon depletion were constant throughout the disk, this would not introduce an additional uncertainty in the employment of CO isotopologues as tracers of the disk surface density distribution.

¹³ C¹⁸O may be a better tracer of Σ_gas(R) in the inner regions for massive disks, circumventing the ¹³CO optical depth issue, as suggested by Zhang et al. (2017). For lower mass disks, C¹⁷ O and C¹⁸ O may be more appropriate. However, both dust optical depth and gas temperature may limit the analysis. Thus, combining observations of optically thin tracers with ¹³CO may be the best option after all.

Fig. 7 Panels a and b: C17O line intensity radial profiles (solid lines) obtained with six disk models with input surface density distribution Σgas (dashed lines) chosen to be a simple power-law (see Eq. (2)). The model parameters are Rout = 200 au, Mdisk = 10−4M⊙ (left panel), 10−2M⊙ (right panel) and γ = 0.8, 1, 1.5 shown in dark blue, light blue, and purple, respectively (top, middle, and bottom panels). Panels c and d: column densities of gas-phase (solid lines) and ice-phase (dashed lines) of C17 O calculated from the surface to the midplane shown as function of the disk radius for six models with simple power-law surface density (Rout = 200 au, Mdisk = 10−4, 10−2M⊙). Top, middle, and bottom panels: models with γ = 0.8, 1, 1.5, respectively. The dotted black lines indicate the column density at which CO self-shielding becomes inefficient (N = 1015cm−2).

Acknowledgements

The authors thank the referee J.P. Williams, S. Andrews, L. Testi, I. Pascucci, and I. Kamp for the comments which helped to improve the paper, and E. Bergin, and K. Schwarz for sharing their data. Astrochemistry in Leiden is supported by the Netherlands Research School for Astronomy (NOVA), by a Royal Netherlands Academy of Arts and Sciences (KNAW) professor prize, and by the European Union A-ERC grant 291141 CHEMPLAN. A.M. acknowledges an ESO Fellowship.

Appendix A: Additional figures

Here some ancillary figures are reported.

Fig. A.1 Abundance of 13CO in one quadrant the disk for two representative models. Panel a: self-similar disk model (Rc = 200 au, Mdisk = 10−3 M⊙, and γ = 1); panel b: simple power-law disk model (Rout = 200 au, Mdisk = 10−4 M⊙, and γ = 1). The white contours indicates the Tdust = 20 K surface below which CO freeze-out becomes important.

Fig. A.2 13CO line intensity radial profiles (solid lines) obtained with input surface density distribution Σgas (dashed lines) given by the simple power-law disk models (see Eq. (1)). The model parameters are Mdisk = 10−5 M⊙ (panel a) and Mdisk = 10−2 M⊙ (panel b), γ = 0.8, 1, 1.5 shown in dark blue, light blue, and purple, respectively (top, middle, and bottom panels).

Fig. A.3 13CO line intensity radial profiles (solid lines) obtained with input surface density distribution Σgas (dashed lines) given by the self-similar disk model (see Eq. (1)). The model parameters are Mdisk = 10−4 M⊙ and γ = 0.8, 1, 1.5 shown in dark blue, light blue, and purple, respectively (top, middle, and bottom panels). Model with Rc = 30, 60, 200 au are presented in panels a, b, and c, respectively.

References

All Tables

All Figures

Fig. 1 Panel a: 13CO line intensity radial profiles (solid lines) obtained with three representative disk models with input surface density distribution Σgas (dashed lines) chosen to be a simple power-law (see Eq. (2)). The model parameters are Rout = 200 au, Mdisk = 10−4 M⊙, and γ = 0.8, 1, 1.5 shown in dark blue, light blue, and purple, respectively (top, middle, and bottom panels). No beam convolution is applied to the simulated images and a distance of 150 pc is assumed. Panel b: 13 CO line intensity radial profiles compared with input surface density profiles as in panel a, but convolved with a 0.2′′ beam. Panel c: column densities of gas-phase 13CO calculated from the surface to the midplane shown as function of the disk radius for the three representative models with simple power-law surface density. The dotted black line indicates the column density at which CO self-shielding becomes inefficient (N = 1015 cm−2). Ice-phase CO column densities are not shown in the plot as they are below the values shown here.

In the text

Fig. 2 Results of the power-law fitting of the simulated 13CO line intensity profiles obtained with the simple power-law models (Mdisk = 10−4 M⊙ and Mdisk = 10−2 M⊙ in panels a and b, respectively). The values of the fitted power-law index γfit are shown by the filled dots as a function of the radial range over which the fitting was carried out for the high-resolution images (0.01′′ beam). The empty symbols in panel a show the fitted power-law index γfit when the simulated images are convolved with a 0.2′′ beam. The model input power-law indexes are shown by the dashed lines. The dotted line shows the starting point of the fitting Rstart.

In the text

Fig. 3 Panel a: 13CO line intensity radial profiles (solid lines) obtained with three representative disk models with input surface density distribution Σgas (dashed lines) given by the self-similar disk model (see Eq. (1)). The model parameters are Rc = 200 au, Mdisk = 10−3 M⊙, and γ = 0.8, 1, 1.5 shown in dark blue, light blue, and purple, respectively (top, middle, and bottom panels). Panel b: column densities of gas-phase (solid lines) and ice-phase (dashed lines) 13 CO calculated from the surface to the midplane shown as function of the disk radius for the three representative self-similar disk models. The dotted black line indicates the column density at which CO self-shielding becomes inefficient (N = 1015 cm−2).

In the text

Fig. 4 Results of the power-law fitting of the simulated 13CO line intensityprofiles obtained with the self-similar disk models with Mdisk = 10−4 M⊙. The values of the fitted power-law index γfit are shown by the filled dots as a function of the radial range in which the fitting was carried out. The model input power-law indexes are shown by the dashed lines. The dotted line shows the starting point of the fitting Rstart .

In the text

Fig. 5 Total 13CO line optical depth τ (panels a and d), 13CO column density (panels b and e), and 13CO intensity radial profiles for simple power-law disk models (top panels) and self-similar disk models (bottom panels). The total disk mass is set to Mdisk = 10−4M⊙ and γ = 1. The shaded regions in panels a and d show the range where τ < 1, which holds for radii larger than those shown by the dotted lines. The shaded region in panels b and e shows the region where 13CO column densities are larger than 1015 cm−2, i.e., for radii smaller than those shown by the dashed lines. The shaded regions in panels c and f show the range of 13 CO line intensityfor which the conditions τ < 1 and N13CO > 1015 cm−2 hold. The fit of γ is performed over the “slope–pivot region” between the dotted and dashed lines, shown by the red arrow in panel c.

In the text

	Fig. 6 Radial intensity profile of the C18O (3–2) line observed with ALMA (Schwarz et al. 2016) is shown in black and the 3Σ level of uncertainty is presented by the shaded region. The red line shows the power-law function used to fit the data at the slope–pivot region.
In the text

Fig. 7 Panels a and b: C17O line intensity radial profiles (solid lines) obtained with six disk models with input surface density distribution Σgas (dashed lines) chosen to be a simple power-law (see Eq. (2)). The model parameters are Rout = 200 au, Mdisk = 10−4M⊙ (left panel), 10−2M⊙ (right panel) and γ = 0.8, 1, 1.5 shown in dark blue, light blue, and purple, respectively (top, middle, and bottom panels). Panels c and d: column densities of gas-phase (solid lines) and ice-phase (dashed lines) of C17 O calculated from the surface to the midplane shown as function of the disk radius for six models with simple power-law surface density (Rout = 200 au, Mdisk = 10−4, 10−2M⊙). Top, middle, and bottom panels: models with γ = 0.8, 1, 1.5, respectively. The dotted black lines indicate the column density at which CO self-shielding becomes inefficient (N = 1015cm−2).

In the text

	Fig. A.1 Abundance of 13CO in one quadrant the disk for two representative models. Panel a: self-similar disk model (Rc = 200 au, Mdisk = 10−3 M⊙, and γ = 1); panel b: simple power-law disk model (Rout = 200 au, Mdisk = 10−4 M⊙, and γ = 1). The white contours indicates the Tdust = 20 K surface below which CO freeze-out becomes important.
In the text

	Fig. A.2 13CO line intensity radial profiles (solid lines) obtained with input surface density distribution Σgas (dashed lines) given by the simple power-law disk models (see Eq. (1)). The model parameters are Mdisk = 10−5 M⊙ (panel a) and Mdisk = 10−2 M⊙ (panel b), γ = 0.8, 1, 1.5 shown in dark blue, light blue, and purple, respectively (top, middle, and bottom panels).
In the text

Fig. A.3 13CO line intensity radial profiles (solid lines) obtained with input surface density distribution Σgas (dashed lines) given by the self-similar disk model (see Eq. (1)). The model parameters are Mdisk = 10−4 M⊙ and γ = 0.8, 1, 1.5 shown in dark blue, light blue, and purple, respectively (top, middle, and bottom panels). Model with Rc = 30, 60, 200 au are presented in panels a, b, and c, respectively.

In the text