# C18O emission as an effective measure of gas masses of protoplanetary disks

Maxime Ruaud  NASA Ames Research Center, Moffett Field, CA, USA  Carl Sagan
Center, SETI Institute, Mountain View, CA, USA  Uma Gorti  NASA Ames Research
Center, Moffett Field, CA, USA  Carl Sagan Center, SETI Institute, Mountain
View, CA, USA  David J. Hollenbach  Carl Sagan Center, SETI Institute,
Mountain View, CA, USA

(Received January 1, 2021; Revised January 7, 2021; Accepted November 8, 2021)

###### Abstract

Many astrochemical models of observed CO isotopologue line emission, earlier
considered a good proxy measure of H2 and hence disk gas mass, favor large
deviations in the carbon and oxygen gas phase abundances and argue that severe
gas phase CO depletion makes it a poor mass tracer. Here, we show that C18O
line emission is an effective measure of the gas mass, and despite its complex
chemistry, a possibly better tracer than HD. Our models are able to reproduce
C18O emission from recent ALMA surveys and the TW Hya disk to within a factor
of $\sim 2-3$ using carbon and oxygen abundances characteristic of the
interstellar medium (C/H$=1.4\times 10^{-4}$; O/H$=3.2\times 10^{-4}$) without
having to invoke unusual chemical processing. Our gas and dust disk structure
calculations considering hydrostatic pressure equilibrium and our treatment of
the CO conversion on grains are primarily responsible for the very different
conclusions on disk masses and CO depletion. As did previous studies, we find
that a gas phase C/O of $\sim 1-2$ can explain observed hydrocarbon emission
from the TW Hya disk; but significantly, we find that CO isotopologue emission
is only marginally affected by the C/O ratio. We therefore conclude that C18O
emission provides estimates of disk masses that are uncertain only to within a
factor of a few, and describe a simplified modeling procedure to obtain gas
disk masses from C18O emission lines.

Protoplanetary disks (1300), Astrochemistry (75), Chemical abundances (224),
CO line emission (264), Star formation (1569)

††journal: ApJ

##  1 Introduction

The question of how much mass is present in protoplanetary disks is
fundamental to understanding disk evolution and planet formation. The rate at
which the total mass declines with time is the primary indicator of the
efficiency of mass and angular momentum transport, and of planet formation and
disk dispersal timescales. The available reservoirs of gas and dust masses,
how they evolve with time and relative to one another, as well as how the disk
chemical composition itself changes are all relevant to determining the
conditions under which rocky and gaseous planets can form.

Dust masses for particles $\lesssim$ 1 cm in size can be determined from
submillimeter continuum emission to within an order of magnitude despite some
uncertainties in opacity, as dust emission from particles with sizes
$\lesssim$ 1 cm is mostly optically thin with the emission proportional to the
dust mass (e.g. Hildebrand, 1983; Ansdell et al., 2016; Pascucci et al., 2016;
Tychoniec et al., 2020).

Gas mass determination, however, has been difficult even to order of magnitude
certainty as there are few optically thin tracers of gas. HD, which is the
isotopologue of the main gas constituent H2, is expected to be optically thin
for typical disks but has been detected only in three disks by Herschel
(Bergin et al., 2013; McClure et al., 2016). Both the lack of facilities
sensitive enough to detect the faint line (e.g. Bergin & Williams, 2017) and
the fact that disk regions containing most of the gas mass are still too cold
to excite the $E_{u}\sim 128$K upper level of HD (e.g. Trapman et al., 2017)
limit our ability to use this molecule as a gas mass tracer. CO is the next
most abundant molecule and its isotopologues are often used as tracer species;
several surveys of disks in star-forming regions have been now conducted with
ALMA measuring both dust continuum emission and C18O line emission (Ansdell et
al., 2016, 2017; Long et al., 2017). However, most recent studies conclude
that C18O may not be a reliable indicator of gas disk mass (e.g. Calahan et
al., 2021b) as ascribed to reasons below.

CO isotopologue emission from disks has been extensively investigated in
recent years. CO emission is bright, and with many accessible transitions and
isotopologues (e.g., 13CO, C18O), can potentially trace a wide range of
temperatures and densities. CO emission can also be spatially mapped using
interferometers such as ALMA and, in principle, a complete survey of CO and
its isotopologues can be used to retrieve and piece together the disk radial
gas surface density distribution (e.g., as has been done for TW Hya or GM Aur,
see Zhang et al., 2017, 2019; Calahan et al., 2021b; Schwarz et al., 2021).
However, CO can freeze out on the surface of cold dust grains, and
characteristic dust temperatures in the outer disk midplanes are such that
these regions are obscured due to a lack of gas phase CO. Even after
accounting for the fact that CO condenses in the midplane (at $T_{d}\sim 20$K,
assuming a CO binding energy of 1150K and no further chemical processing in
the ice), models have typically overproduced C18O fluxes leading to
speculations about a variable CO/H2 ratio in disks. However, a key issue is
that CO emission does not probe total disk mass directly—optically thin CO
isotopologues are photodissociated at the surface and CO freezes out near the
midplane—and models are needed to correlate the observed flux with the
emission; the conclusions made are thus dependent on the model assumptions.
Rarer isotopologues involve multiplications by larger factors to convert to
H2, and even minor uncertainties in their isotopic fractionation chemistry can
get amplified to make such gas mass determinations unreliable.

Recent ALMA suveys have detected line fluxes for CO and its isotopologues and
these are far lower than many model predictions that assume canonical ISM
CO/H2 and dust/gas mass ratios (Ansdell et al., 2016; Miotello et al., 2017;
Long et al., 2017). These low fluxes are usually interpreted using two
scenarios: (1) either that CO is a good measure of gas mass and there is a
significant gas removal leading to low gas masses in these disks, (2) or that
CO is severely under-abundant throughout these disks (i.e. not only in the
disk midplane where CO freezes-out) and that gas phase CO is depleted even at
the surface by some efficient mechanism. For (1), Ansdell et al. (2016), using
the simple model for disk structure and chemistry from Williams & Best (2014)
(developed to interpret SMA data of a survey of 15 disks), conclude that most
observed disks in Lupus have low gas masses, with dust/gas mass ratios $\sim
1$ which implies a factor of 100 reduction from that in the parent cloud
material, indicating relatively rapid disk gas dispersal ; which however is at
odds with the relatively high mass accretion rates inferred in the region
(Manara et al., 2016; Alcalá et al., 2017). These were revised downwards by a
factor of $\sim 10$ using more sophisticated disk models that include
selective photodissociation of C18O (Miotello et al., 2014, 2016, 2017).
However, Herschel observations of HD $(1-0)$ in a few disks find dust/gas mass
ratios that are compatible with ISM values (Bergin et al., 2013; McClure et
al., 2016). These findings have led many authors to favor (2) and therefore
invoke large ($\sim 2$ orders of magnitude) depletion of CO that must act in
addition to processes like CO freeze-out in the cold midplane (e.g. Favre et
al., 2013; Du et al., 2015; Schwarz et al., 2016; Kama et al., 2016; Trapman
et al., 2017). This deficit in observed CO isotopologue emission (compared to
observed HD emission and model CO predictions) has been explained as being due
to a change in C/H and C/O ratios in disks as they evolve (e.g. Du et al.,
2015; Bergin et al., 2016; Kama et al., 2016; Schwarz et al., 2016; Xu et al.,
2017; Krijt et al., 2016, 2018, 2020; Zhang et al., 2019; McClure, 2019;
Calahan et al., 2021b; Zhang et al., 2021; Schwarz et al., 2021) or as the
result of the incorporation of CO into solids on the form of more complex
compounds at the surface of the grains such as CO2, CH3OH and carbon chains
(e.g. Bergin et al., 2014; Furuya & Aikawa, 2014; Aikawa et al., 2015;
Reboussin et al., 2015; Molyarova et al., 2017a; Eistrup et al., 2016; Bosman
et al., 2018; Ruaud & Gorti, 2019). However, the more widely held
interpretation is that pebble formation in disks locks up available carbon and
oxygen and gradually depletes gas phase CO (Krijt et al., 2018, 2020). This is
then believed to change the global C/O and C/H ratios in disks, lowering the
CO gas phase abundance in emitting regions from canonical ISM values. A more
complete discussion of the chemistry involved and the implications for gas
mass is presented in a recent review by Öberg & Bergin (2021).

In this paper, we examine the global CO isotopologue depletion issue in
greater detail and argue against the need for both scenarios (1) and (2). In
order to determine the extent to which CO isotopologue emission can measure
disk mass, one needs to determine both the extent to which C/O changes and to
which C and O abundances need to be altered to produce the observed emission.
We revisit the problem with newer, more complete models of disk chemistry
(Ruaud & Gorti, 2019). Our main motivation for this study is that none of the
previous models have simultaneously considered all the processes critical to
interpreting CO emission (i) a full treatment of grain chemistry and CO
freeze-out, (ii) CO isotopologue selective photodissociation, and (iii) a
self-consistent gas disk density and temperature structure coupled by vertical
pressure equilibrium. Using models that incorporate all these physical and
chemical processes, we conduct a parameter survey where we vary gas disk
masses and the dust/gas mass ratio. We investigate the effect of grain surface
chemistry on CO abundance, the dependence of CO emission on disk parameters,
and also evaluate other mass tracers like HD and attempt to reconcile model
results with available observational data. We show that C18O emission is well
explained with disks having having ISM like elemental abundances and argue
that C18O is a good tracer of disk gas mass over a wide range of disk masses.
We also show that observed CO emission does not by itself indicate any change
in the C/O ratio in disks, and discuss the key model assumptions necessary to
retrieve disk masses from CO isotopologue emission. The structure of the paper
is as follows. After a brief description of the models in §2, the main results
are presented in §3, followed by a comparison with existing models and data in
§4 and prescriptions for simple analyses in §5. Conclusions are presented in
§6.

##  2 Model description

We use the modeling framework described in Ruaud & Gorti (2019) (hereafter
RG19) to compute the disk physical structure and chemistry. A brief summary of
the modeling procedure is provided here and we refer the reader to RG19 for
further details. Stellar properties (mass, spectrum, UV and X-ray spectra) and
the disk gas surface density ($\Sigma_{\mathrm{gas}}(r)$) and dust surface
density ($\Sigma_{\mathrm{dust}}(r)$) are the input parameters. The disk
physical structure is determined by solving for vertical hydrostatic gas
pressure equilibrium simultaneously with a surface density constraint;

$\frac{dP(r,z)}{dz}=-\rho(r,z)\Omega^{2}z;\quad\Sigma_{\mathrm{gas}}(r)=\int\rho(r,z)dz$
(1)

where the pressure $P(r,z)=\rho(r,z)k_{B}T(r,z)/\mu(r,z)$ and
$\rho(r,z),T(r,z),\mu(r,z)$ and $\Omega$ are the gas density, temperature,
mean molecular weight and Keplerian frequency respectively. Gas heating
includes dust collisions, UV, X-rays and other mechanisms, cooling is by line
and continuum emission. At each $(r,z)$ the gas density, temperature and
chemistry and the dust size distribution are all coupled and we find a
consistent solution that matches the local pressure gradient and results in
the specified $\Sigma_{\mathrm{gas}}(r)$ and $\Sigma_{\mathrm{dust}}(r)$. At a
given radius $r$, the dust grain size distribution is determined assuming a
balance between fragmentation and coagulation processes, and dust settling is
treated by calculating the dust scaleheight which is a function of grain size
and local turbulence parameter (for more details see Gorti et al., 2015) and
for which we use a viscous parameter $\alpha=5\times 10^{-3}$. The dust
temperature is computed for each grain (size and composition) at every
$(r,z)$. Since the dust and gas are collisionally coupled with exchange of
momentum and thermal energy, the dust and gas disk structures are determined
together and a solution is iteratively obtained to be consistent with the
constraints set by Equation 1.

The density and temperature structure obtained at every $(r,z)$ of both gas
and dust is used to solve for gas-grain chemistry. In this model, the ice
surface and ice mantle are treated as two separate phases in interaction, and
grain surface chemistry includes diffusion of the chemical species on the
surface, two-body reactions, photoreactions, and thermal and nonthermal
desorption mechanisms (see RG19 and references therein for more details). The
gas-phase chemical network has been extended to account for HD chemistry and
shielding (Bergin et al., 2013) and single and multiple isotope reactions of
13C and 18O. We include reactions listed in Roueff et al. (2015) for 13C and
Loison et al. (2019) for 18O. For all other reaction rates we duplicate the
main isotopologue network and assume similar rate constants and statistical
branching ratios. Isotopic species are restricted to molecules with a maximum
of three carbons and two oxygens. Selective photo-dissociation of CO
isotopologues is considered using the self-shielding parameters from Visser et
al. (2009). Finally, we assume similar binding energies for all ices formed
from C and O isotopes. In total, the network consists of $\sim 800$ species
and $\sim 10^{4}$ reactions.

Emission line fluxes are computed using an escape probability method which is
the same as that used for the cooling calculations in the thermal balance
computations and assuming face-on disks (zero inclination). Non-LTE radiative
transfer includes collisional and radiative processes, and includes radiation
from the dust background (e.g. Hollenbach & McKee, 1979; Tielens & Hollenbach,
1985), which is iteratively computed with the disk physical and chemical
structure. Molecular data are from the Leiden Atomic and Molecular Database
(LAMDA). We use the escape probability method for the analysis of line fluxes
as this allows us to compute the continuum background and treat spatial
variations in dust opacities (relevant for HD emission); moreover, because
this method is also used to compute cooling in the disk structure modeling we
can directly compare the emission flux at a given location $(r,z)$ to the
chemistry, heating and cooling processes.

Using the above models, we conducted a limited parameter survey where the disk
mass is varied in increments of 0.5 dex in the range $3\times
10^{-4}-10^{-1}$M⊙ (i.e. $0.3-100$ MJ) for three different dust/gas mass
ratios $\epsilon=10^{-1},10^{-2}$ and $10^{-3}$. All these models assume an
outer disk radius of 300au. For the canonical value of $\epsilon=10^{-2}$, we
run an additional series where the outer disk radius is set to 100au,
resulting in higher surface densities for a given mass when compared to the
300au series. In total, we run 24 different models. In all cases, the stellar
parameters are kept constant with $M_{*}=1.0M_{\odot},R_{*}=2.0$ R⊙ and UV and
X-ray spectra similar to TW Hya (see RG19). We also assume that there is an
ambient UV field with $G_{0}=1$ in Habing units. Cosmic ray ionization rates
are depth dependent and calculated using the low energy fits derived by
Padovani et al. (2018). Table 1 lists the various assumed model parameters. As
emission lines originate from the surface layers, where photoprocesses drive
the chemistry on short timescales to cause a rapid reset, we assume all
species are initially in atoms and ions at the start of the chemical
simulation for computational efficiency. Chemistry in the disk may possibly
retain some signature of the parent cloud core chemistry (see also RG19), but
this is restricted to the inner midplane regions of the most massive disks and
these regions are not emissive for our targeted lines, further justifying our
assumption of the initial state.

Parameter | Value | Parameter | Value
---|---|---|---
Stellar Mass | 1.0M⊙ | Stellar radius | 2.0 R⊙
Stellar UV |  $10^{31}$ erg s-1 | Stellar X-rays |  $10^{30}$ erg s-1
Dust silicate mass fraction | 0.83 | Dust carbon mass fraction | 0.17
Surface density function | $\Sigma(r)=\Sigma_{0}(r_{c}/r)e^{-r/rc}$ | $r_{c}$ | 100 au
Inner disk radius | 1 au | Outer disk radius | 100 and 300 au
Element | Abundance | Element | Abundance
He | 0.1 | N | 6.2(-5)
H | 1.0 | Si | 1.7(-6)
D | 1.6(-5) | Mg | 1.1(-6)
12C | 1.4(-4) | S | 1.0(-6)
13C | 2.0(-6) | Fe | 1.7(-7)
16O | 3.2(-4) | PAH | 1.00(-9)
18O | 5.74(-7) |  |
Table 1: List of the assumed physical parameters for the disk models and gas-
phase elemental abundances used in the chemical modeling. Gas disk masses and
the dust/gas ratio are varied as described in the text.

##  3 HD and CO isotopologue line emission from models

Panels in the top row of Fig. 1 shows the physical structure of one of the
model disks with a mass $10^{-2}$ M⊙ and a dust/gas ratio $\epsilon=10^{-2}$.
HD and CO isotopologue emission typically originates in regions where the
$A_{V}$ to the star is between $\sim 1-10$. Densities in this region are $\sim
10^{7-9}$ cm-3 and gas and dust temperatures are quite similar for the model
disk shown in Fig. 1. The average dust cross-section per H is reduced by a
factor of nearly 100 in the $1<A_{v}<10$ region compared to the ISM value of
$\sim 5\times 10^{-21}$ cm-2 due to the larger dust grain sizes in disks and
due to settling. The disk structure varies with disk gas and dust masses as
discussed in RG19 in greater detail.

In the disk regions responsible for production of the line emission discussed
here, gas is strongly coupled to the dust via collisions with smaller, hotter
dust grains heating the gas and larger, colder grains cooling it. Dust
collisions are typically a net coolant for the more massive disks with
additional gas heating provided by FUV grain photoelectric heating, and gas is
warmer than the dust. For the lower mass disks, dust collisions are typically
a net heating source with additional cooling provided mainly by CO and
[OI]63$\mu$m emission and in this case gas is slightly cooler than the dust.
For details on the included heating, cooling processes and how they affect
density/temperature distributions, please see Gorti et al. (2011) and RG19.

###  3.1 HD and CO isotopologue chemistry

Figure 1:  Computed physical structure, density maps and line luminosities of
a selection of molecules for a model disk with $\epsilon$ $=10^{-2}$ and
$M_{\text{gas}}=10^{-2}M_{\odot}$. The top panel shows the gas density
($n_{\text{H}}$), gas temperature ($T_{\text{g}}$), the area-weighted mean
dust temperature ($T_{\text{d}}$) and the dust cross-sectional area and the
dust density averaged over the grain size distribution
($\langle\sigma_{\text{d}}n_{\text{d}}\rangle$) as a function of the disk
radius $r$ and the normalized height $z/r$. Solid line contours show the
visual extinction ($A_{V}$) to the star. The middle and lower panels show the
computed density maps of HD, 13CO and C18O and the predicted line emission map
normalised by the total emission of the line (J$=1-0$ for HD, J$=3-2$ for CO
isotopologues).

Figure 1 also shows the computed density and emission maps of HD and the three
CO isotopologues considered here. Although the H/HD transition occurs at lower
$z$ than the H/H2 transition, HD chemistry is similar to that of H2. HD forms
on grains and is photodissociated by UV photons; other reactions involving D
or D+ exchanges are not found to be important in setting the abundance of HD
in the disk. Gas phase HD abundance is also not affected by grain surface
chemistry due to the low binding energy of the HD molecule on grains and HD is
therefore abundant even in the midplane of the disk.

The vertical column densities of gas phase CO and its isotopologues are
determined by the location of the CO photodissociation front at $A_{v}\sim
10^{-2}$ and the height where water ice begins to form at $A_{v}\gtrsim 2-3$.
Gas-phase CO is present in between these two layers (see Fig. 1). In RG19, we
show that photoprocessing of ice, near the surface corresponding to the
vertical water snowline, induces efficient conversion of CO into CO2 ice and
impacts the location of the vertical CO snowline. In this process,
photodissociation of water ice near its snowline leads to the production of OH
which rapidly reacts with CO on grains (before it can desorb) to form the more
bound CO2 ice. The location of the CO snowline is thus determined not by the
CO desorption/condensation balance but by desorption of CO2 that the CO ice
converts to. The vertical CO snowline is therefore shifted higher up from the
midplane because of the higher binding energy of CO2 relative to CO, and the
CO snowline more or less coincides with the location of the vertical water
snowline which is itself determined by photodesorption (see Furuya & Aikawa,
2014; Aikawa et al., 2015; Molyarova et al., 2017a; Bosman et al., 2018;
Trapman et al., 2021, for similar conclusions). This is illustrated in Fig. 2
which shows a comparison of the vertical chemical structure computed when this
effect is taken into account and when it is neglected. The gas density at the
CO snowline location depends on the dust cross section per H nucleus which can
be quite low in disks due to both grain growth and settling. Along the radial
direction, and in the inner disk midplane, UV photons are produced by cosmic
rays and the low UV field results in a CO to CO2 conversion on longer ($\sim$
Myr) timescales. As a result, the CO snowline slowly moves inward until most
of the CO has been converted to CO2 ice in a few Myr111As will be discussed in
§5, the location of the snowline as a function of time can be obtained by
solving the reduced network given in Appendix A (for more details, see RG19).

Figure 2:  Computed vertical gas-phase density of CO and ice density of CO and
CO2 (denoted as sX) as a function of $z/r$ and at $r=50$ au of the disk shown
in Fig. 1. The top panel shows results when CO conversion to CO2 ice was
considered and the bottom panel when this conversion was neglected. Black
dashed lines show the dust temperature. When CO conversion into CO2 ice is
included, the location of the CO snowline is at higher $z$ because of the
higher binding energy of CO2 relative to CO. For this disk and at this radius,
this effect results in all CO being converted to CO2 ice even in the midplane
due to cosmic ray ionization (see Ruaud & Gorti, 2019), and a $\sim 25\%$
decrease of the vertical gas CO column density as compared to the case where
this conversion is neglected.

The vertical extent of gas phase CO and its isotopologues is determined not
only by their condensation front but also by where they are photodissociated
in the surface layers. In the disk atmosphere (i.e. at $z/r\gtrsim 0.4$ for
the model shown in Fig. 1), all CO isotopologues are subject to
photodissociation by stellar and interstellar FUV radiation. Self-shielding
prevents photodissociation of CO and its isotopologues and they begin to rise
in abundance when $z/r\lesssim 0.4$. 13CO chemistry is not particularly
sensitive to selective photodissociation and is closely coupled to 12CO
chemistry through the exchange reaction

${}^{13}\mathrm{C}^{+}+\mathrm{CO}\leftrightarrow\mathrm{C}^{+}+{{}^{13}\mathrm{CO}}+\Delta
E=35K$ (2)

which enhances the abundance of 13CO compared to 12CO (relative to that
expected from selective dissociation effects alone). However, isotope
selective photodissociation has a significant impact on the distribution of
C18O which is not shielded by 12CO or 13CO and is abundant only in a thin
layer located just above the vertical CO snowline (at $z/r\sim 0.15$ for the
model shown in Fig. 1) where it forms by reaction of 18O with CH and CH2.

###  3.2 HD and CO isotopologue line emission

The lower panels of Figure 1 show contributions, as a function of spatial
location, to the corresponding emission of HD $J=1-0$ and the $J=3-2$ line of
the CO isotopologues. For CO isotopologues we focus on the emission of the
$J=3-2$ transition which is usually brighter in disks.

HD traces the main constituent H2 and is abundant throughout the disk,
including the dense midplane regions; however, emission from the $1-0$ line
arises mainly from the warmer upper layers of the disk, with a small region of
emission in the inner disk where the gas is sufficiently warm. In both these
regions, one at the inner disk midplane (at $z/r\lesssim 0.1$ and $r\lesssim
20$ au in Fig. 1) and the other above the disk midplane (at $z/r\gtrsim 0.15$
and $r\gtrsim 20$ au), the temperature of the gas ranges between 30 K and 80
K. As has been noted by previous authors (Bergin et al., 2013; Trapman et al.,
2017; Kama et al., 2020), the relatively higher energy level of the HD $1-0$
transition ($E_{u}\sim 128$K) makes HD emission quite sensitive to the gas
temperature at the cold temperatures typical of disks because the fractional
population in the $J=1$ level, and hence emission of the $J=1-0$ line, is
proportional to $\exp(-128/T_{\text{gas}})$. HD is moreover suppressed by the
strong continuum background in disks at $112\mu$m, which lowers the line flux
that can be detected above the continuum. There is thus a relatively sharp
drop in HD emission at the height where the dust temperature transitions from
the optically thin value (the so-called superheated layer, e.g., Chiang &
Goldreich, 1997). The HD $J=1-0$ line therefore is sensitive to the local
temperature and density structure and originates mainly from the surface
layers of the disk.

Emission of all the CO lines is also confined to the surface, and specifically
to the spatial regions between the CO snowline and the photodissociation
fronts where CO and its isotopologues exist in the gas phase. As is to be
expected, 12CO is optically thick and the $J=3-2$ transition shown in Fig. 1
peaks in a narrow layer near the CO photodissociation front. Depending on the
disk mass and gas column density, 13CO $J=3-2$ could be optically thin but for
the model disk shown here the line is marginally optically thick. As discussed
earlier, 13CO is not significantly affected by selective isotope
photodissociation and the emitting region is therefore more extended. The line
luminosity is dominated by the regions closer to the snowline and with
$r\gtrsim 150$ au for this model disk. Contrary to 13CO, C18O emission mostly
originates from regions with $r\lesssim$ 150 au, where the density is higher,
and is optically thin at most radii even for disks more massive than the
$10^{-2}$M⊙ disk shown here. The emission, like that of HD, is therefore
weighted toward higher densities (lower $z$) and thus higher column densities,
and regions located just above the CO snowline dominate the emergent flux. We
note that the gas temperature, which is within a few K of the dust
temperature, is typically $\gtrsim 30$K above the CO snowline where CO is in
the gas phase. An important consequence therefore is that while both HD and
C18O luminosities are sensitive to the gas density, unlike HD, the lower
energies of the C18O sub-mm transitions make C18O luminosities relatively
insensitive to the disk temperature structure.

###  3.3 Disk mass parameter survey

Figure 3:  HD (1-0),12CO 3-2, 13CO (3-2), C18O(3-2) emission as function of
gas mass for different dust/gas ratios ($\epsilon=10^{-2}$:blue,
$\epsilon=10^{-3}$:green and $\epsilon=10^{-1}$:red). Cyan color data show the
$\epsilon=10^{-2}$ model disks 100au in radial extent and are to be compared
with the blue $\epsilon=10^{-2}$ dataset with radii of 300au. The differences
in emission for the more optically thin lines are mostly due to differences in
columns in the emitting regions. For optically thick lines such as 12CO and
13CO (which is marginally optically thick), differences are mainly due to line
optical depth.

Figure 3 shows the computed luminosity of HD ($1-0$), 12CO, 13CO and
C18O($3-2$) lines for all the models of the parameter survey discussed in §2
as a function of the gas disk mass. We consider three values of the dust/gas
ratio; this is fixed as a function of radius (i.e.,
$\Sigma_{\mathrm{dust}}(r)/\Sigma_{\mathrm{gas}}(r)=\epsilon=10^{-1},10^{-2}$
and $10^{-3}$) but varies with height due to dust settling. We also consider
one set of models with $\epsilon=10^{-2}$ but with a smaller radial extent.
All line luminosities exhibit an increasing trend with gas mass, with a slope
that is shallow for the more optically thick 12CO and 13CO lines and steeper
for HD and C18O. 12CO emission is relatively insensitive to the dust/gas
ratio, mainly because it originates high above the disk (see fig. 1) at unit
line optical depth where the gas temperature is determined by FUV and X-ray
heating mechanisms and line cooling. Because the emission is optically thick,
the CO luminosity depends on the disk spatial extent and decreases
substantially for the series with smaller disk size. 13CO line emission is
marginally optically thick and in general exhibits trends intermediate to very
optically thick 12CO emission and optically thin HD and C18O emission.
Emission from HD and C18O is optically thin at most disk radii; the former
depends on the column density of warm HD and the latter on the column density
of CO above the vertical snowline respectively, and these in turn depend on
the disk physical structure as we describe below.

#### Gas mass variations for a fixed $\epsilon$:

The optically thin emission does not scale linearly with gas disk mass (for a
fixed $\epsilon$) as might be expected because the local gas density and
resulting column densities also depends on the disk structure—denser gas is in
general cooler which affects the emission. This was previously noted by
Trapman et al. (2017) who also found a sublinear change in HD emission with
gas mass. C18O emission, being less sensitive to temperature by virtue of its
lower excitation temperature as previously discussed, exhibits a steeper trend
with gas mass compared to HD, although still sub-linear (i.e. C18O 3-2
luminosity approximately scales with $\propto M_{\text{gas}}^{0.7-0.8}$ while
HD 1-0 luminosity is $\propto M_{\text{gas}}^{0.6-0.7}$). The emitting regions
for HD and C18O approximately lie near the $A_{V}=1$ layer in the disk, at a
few scaleheights. While the overall density in the disk is proportionately
lowered when the disk mass is decreased, a shift in the emitting regions to
denser gas at lower $z$ compensates to produce the sublinear trend with disk
mass. This shift is mainly due to decreased dust/gas coupling in less dense
disks; this causes greater settling of dust to move the $A_{V}=1$ layer and
affects the difference in dust and gas temperatures.

#### Line luminosities with $\epsilon$ variations:

The disk surface ($A_{V}\sim 1$ layer) and flaring are sensitive to the disk
structure, which is again affected by the dust/gas coupling which is
determined by $\epsilon$. The gas temperature is nearly the same as the
(optically thin) dust temperature at the surface in the dense, thermally well
coupled, HD and C18O emitting regions and does not vary much for the different
$\epsilon$ models.For a given disk mass, however, high values of $\epsilon$
move the emitting layer to higher $z$ and lower gas densities affecting the
optically thin emission lines of HD and C18O. For the
$M_{\text{gas}}=10^{-1}$M⊙ disk, both the HD and C18O luminosities increase by
$\sim 30$ when $\epsilon$ changes from $10^{-1}$ to $10^{-3}$. At lower gas
masses, this trend with $\epsilon$ continues for the HD line. However, C18O
line luminosity slope as a function of gas mass is shallower for dustier disks
(red curve in Fig. 3), because in these disks the increased dust (relative to
gas) contributes more toward shielding C18O from FUV photodissociation,
thereby increasing the abundance of C18O and producing more line emission.

Figure 4:  Gas mass fraction traced by HD and C18O as a function of the total
gas mass of each disk, calculated by weighting the mass in a given volume
element by the fractional contribution to the line luminosity from that
element. For both lines, the mass traced is typically less than $\sim 10^{-3}$
or $0.1$% of the total disk mass for most of the cases.

#### Fraction of gas mass traced by HD and C18O line emission

We next compute the fraction of the gas mass that is traced by the HD (1-0)
and C18O (3-2) emission, i.e., we sum the gas mass in each cell at $(r,z)$ by
weighting it by the fractional luminosity contribution from that cell. These
results are shown in Fig. 4. The most striking result apparent from these
figures is that HD and C18O both trace a small fraction ranging between $\sim
10^{-6}-10^{-3}$ of the total gas mass for all the various disk models
considered here. This is mainly due to the fact that the dense midplane is not
emissive; most of the HD and almost all the C18O emission comes from the disk
surface. If the optically thin dust disk surface ($\sim A_{V}=1$) is at some
height $z_{t}=\chi h$, where $h$ is the scaleheight corresponding to the
temperature at the midplane, then the fraction of mass contained above this
height ($z=z_{t}$ to $z=\infty$) is equal to
$\Sigma_{s}(z>z_{t})/\Sigma=1-\mathrm{erf}(\chi/\sqrt{2})$ which is $\sim
5\times 10^{-7}-3\times 10^{-3}$ for values of $\chi=3-5$. This mass fraction
again depends on the disk structure (which determines $\chi$ and $h$) and the
trends observed in Fig. 4 can be traced back to those in Fig. 3, because lower
disk masses in general result in lower values of $\chi$.

Our results therefore indicate that C18O is just as effective in tracing the
disk mass as HD, and in fact could be a better tracer. HD and C18O emission
both depend on the gas density, but C18O is, in addition, relatively
insensitive to the temperature which is $\sim 30-80$K in the emitting regions,
i.e., HD 1-0 emission is proportional to $\exp(-128/T_{\text{gas}})$ while
C18O 3-2 lies at $E_{u}\sim 32$K, and is proportional to
$\exp(-32/T_{\text{gas}})$. Conversion of HD flux to total gas mass, on the
other hand, requires that the temperature structure of the gas be sufficiently
well known. Further, multiple rotational C18O lines are readily accessible by
ALMA, NOEMA, SMA and other sub-mm facilities, breaking some model
degeneracies, whereas the HD (1-0) line at 112$\mu$m is currently
inaccessible. We note here that the results presented above do not invoke any
CO depletion factors and assume standard values (as in the ISM) for the
elemental C and O abundances. We next compare our results with observations,
our discrepancies with several previous works, and possible reasons.

##  4 Comparison with ALMA surveys and previous chemical modeling analyses

###  4.1 Agreement with ALMA CO surveys

|
---|---
Figure 5: 13CO (2-1) and (3-2) emission as function of dust mass for different dust/gas ratios, assuming zero inclination for the disks. Black points show 13CO (3-2) results obtained from the Lupus survey (Ansdell et al., 2016). The point in magenta shows a result obtained from the Cham I survey (Long et al., 2017). The cyan model points are for a smaller 100 AU disk with $\epsilon=10^{-2}$. |
---|---
Figure 6: C18O (2-1) and (3-2) emission as function of dust mass for different
dust/gas ratios. Black points show C18O (3-2) results obtained from the Lupus
survey (Ansdell et al., 2016). The point in magenta shows a result obtained
from the Cham I survey (Long et al., 2017). The cyan model points are for a
smaller 100 AU disk with $\epsilon=10^{-2}$.

In Fig. 5 and 6 we compare our computed CO luminosities with luminosities
derived from surveys in Lupus (Ansdell et al., 2016) and Chamaeleon I (Long et
al., 2017), which are plotted as a function of the inferred mm dust masses. We
focus on disks for which the potential gas mass tracer C18O has been detected,
that is to say, 11 disks for the Lupus survey (Ansdell et al., 2016) and only
one disk for the Chamaeleon I survey (Long et al., 2017). Note that 12CO
luminosities from these surveys are not available in the literature. Of all
the model series, the disk models that reproduce observations the best are
those with a normal interstellar dust/gas ratio and ISM like elemental
abundances. For this set of models, the model predicted line luminosities of
13CO and C18O can explain the ALMA data to within a factor of a few. However,
note that while dust masses for the models are input parameters and hence
known, we have simply used quoted dust disk masses to locate observational
data in the figures and that these include inherent uncertainties in the
continuum flux to dust mass conversion (by up to a factor of $\sim 5$, see
e.g. Ballering & Eisner, 2019; Tychoniec et al., 2020; Ribas et al., 2020;
Macías et al., 2021). Comparison of 13CO and C18O luminosities show that the
models slightly overpredict 13CO or underpredict C18O by a factor of $\sim 2$.
Observed line ratios give 13CO/C18O$\sim 2-4$ while modeled ratios vary
between $3-10$ for all the model series.

Without access to 12CO observations in these disks, which could have helped in
distinguishing between these two possibilities, we first explore the
possibility of 13CO overprediction. Our investigations show that it is
extremely difficult to decrease 13CO emission by a factor of $\sim 2$ by
modifying chemistry without also affecting the emission of 12CO and C18O. In
the models, 13CO chemistry is completely dominated by the exchange reaction
with 12CO (see also e.g. Miotello et al., 2014). Varying input parameters like
stellar UV flux, X-ray luminosity, metal abundances and PAH initial
abundances, or disk surface density distributions, did not affect this result.
Part of the 13CO overprediction compared to ALMA data can be explained by the
fact that our model flux computations assume face-on, zero inclination disks,
while there are a range of inclinations likely present in the observed survey
samples. Whereas inclination will not affect the optically thin C18O line, the
13CO line is marginally optically thick (see Fig. 3) and inclined disks are
expected to have lower observed luminosities. However, disk inclination is
likely to suppress 13CO by less than a factor of 2 on average, and even if
this is included, our models will still overestimate 13CO by a small factor.

Another possibility for decreasing 13CO emission while keeping C18O emission
roughly constant would be to truncate the disk. As shown seen in Fig. 1, 13CO
emission originates from the outer disk (i.e. $r\gtrsim 150$ au) while C18O
emission originates in regions with $r\lesssim 150$ au. However, as seen in
seen in Fig. 5 and 6 for disks with $\epsilon=10^{-2}$, decreasing the disk
radii from $R=300$ au and $R=100$ au decreases 13CO emission by a factor $\sim
3-5$ but it also decreases C18O emission by a factor of $\sim 2-3$. Thus a
more compact disk can only partly explain the factor of $\sim 2$
overprediction in the 13CO emission. Other possibilities include currently
missing chemical pathways able to compete with the exchange reaction with 12CO
and thus lower 13CO abundance in the disk (see also Section 4.3). We note that
recent work by Trapman et al. (2021) also find that 13CO is over-produced by a
factor of $\sim 2$ by their models.

Investigations on the possibility of C18O underprediction show that increasing
C18O emission by a factor $\sim 2$ could be achieved by small variations in
the chemistry. Only a small fraction of the available 18O is in the form of
C18O in the models and photochemistry dominates its destruction. As a result,
a small increase in the formation rate of C18O can help in increasing C18O
column density and resulting emission. As discussed in §3, C18O is
predominantly formed by reaction of 18O with CH and CH2 and this rate has not
explicitly been measured. In our model we assume that the rates of reaction
with 18O are similar to that with 16O, for which rate constants are available
only at room temperature or above (Loison et al., 2014). A test for the
$\epsilon=10^{-2}$ model series shows that increasing these rates by a factor
of 2 is enough to enhance C18O emission by a factor of $\sim 1.5$. The
conversion rate of CO into CO2 is also uncertain; in this case, decreasing the
conversion of CO into CO2 ice could increase the abundance of gas phase CO
close to the vertical CO snowline and result in an increase in the C18O column
density and associated emission. A test on the $\epsilon=10^{-2}$ models shows
that completely neglecting the conversion of CO ice into CO2 ice enhances C18O
emission by a factor of $\sim 1.5-2$. Finally, because photochemistry
dominates the destruction of C18O, uncertainties in the C18O shielding factors
could also affect the computed C18O emission.

To summarize, unlike the conclusions from associated modeling in previous
studies (e.g. Ansdell et al., 2016; Miotello et al., 2017), our results can
recover the observed 13CO and C18O line fluxes without invoking increased
dust/gas ratios or CO depletion factors and are only discrepant by a small
factor of $\sim 2$ from the observations. As figure 6 shows, C18O emission can
be explained with disks having ISM like elemental abundances (see Table 1) for
C and O, and with ISM like dust/gas ratios, $\epsilon=10^{-2}$.

###  4.2 Discrepancies with earlier models

Our results substantially differ from those obtained from other similar
thermochemical studies, and here we investigate possible reasons for these
differences. Some of the earliest modeling (e.g., Favre et al., 2013)
calculated C18O fluxes higher than observed and reported CO being
underabundant by more than two orders of magnitude. Later studies showed that
isotope selective photodissociation affects the CO isotopologue fluxes, and
that CO underabundances (or alternately, gas depletion) could be only about an
order of magnitude (Miotello et al., 2016, 2017). However, our results do not
indicate any order-of-magnitude discrepancies.

While it is difficult to directly assess the key differences without formal
bechmarking due to the complexity of the codes, we attempt to identify some
differences below which apply to many currently used thermochemical codes
including DALI (Bruderer et al., 2012), ProDiMO (Woitke et al., 2016) and
RAC2D (Du et al., 2015). Our analysis revealed the following main differences:

  * -

Most models (including DALI, ProDiMO and RAC2D) use a fixed isothermal
Gaussian profile for the density structure of the disk in computing
hydrostatic pressure equilibrium and do not update this to be consistent with
the chemistry and their computed gas temperatures.

  * -

The dust disk structure is not coupled to the vertical gas density or pressure
profile, instead key quantities like the scaleheight and flaring angle are
parameterized independently.

  * -

In many cases, a simplified grain surface chemistry is used where only
accretion, desorption and hydrogenation reactions are taken into account. CO
is assumed to freeze on grains when dust is at $\sim 20$K, with no additional
chemistry. Some models further do not treat selective photodissociation of CO
isotopologues.

We next investigate how the above assumptions affect the calculated C18O model
fluxes. C18O emission is mostly optically thin, and therefore, the luminosity
in a line is determined by the population of the higher level of the
transition. For the $J=3-2$ line, the upper level population follows
$N_{J=3}=N_{\mathrm{C^{18}O}}\exp(-32\text{K}/T_{\mathrm{gas}})/Z(T_{\mathrm{gas}})\propto
N_{\mathrm{C^{18}O}}\exp(-32/T_{\mathrm{gas}})/T_{\mathrm{gas}}$; where
$N_{\mathrm{C^{18}O}}$ is the total number of C18O molecules and $Z(T)$ is the
partition function. For typical emitting regions $T_{\mathrm{gas}}\sim
20-100$K, and the level population and hence line luminosity, peaks near the
excitation temperature ($T\sim 32$K) and weakly declines with gas temperature.
The emission is, therefore, relatively insensitive to gas temperature. The
column density of C18O molecules is consequently the only likely factor
responsible for the differences in C18O emission in our models compared other
studies.

As discussed in §3.2, the C18O column density is set by the CO snowline (at
$z_{\mathrm{ice}}$) and the C18O photodissociation layer (at
$z_{\mathrm{PDR}}$) and is equal to
$\int_{z_{\mathrm{ice}}}^{z_{\mathrm{PDR}}}n_{\mathrm{C^{18}O}}\ dz$. Miotello
et al. (2014) studied the effects of isotope selective photodissociation and
determined that this effect reduces the column density of C18O (by lowering
$z_{\mathrm{PDR}}$) and resulting emission by almost an order of magnitude.
The dominant contribution to the C18O columns however arises from the lower
heights due to the rapid decline in density with $z$. For models that include
isotope selective photodissociation, the location of the CO snowline and the
density here are therefore the main discrepant factors.

The significance of the CO snowline is evident; the assumption that CO freezes
at dust temperature $\sim$ 20K vs. 35K due to conversion to CO2 ice moves the
snowline deeper into the disk midplane where densities are high and will
result in a clear over-estimate of the C18O column density at a given radius
(see Fig. 2) and over-estimate the C18O luminosity. We find that this chemical
conversion effect is one key difference between our models and earlier
results. Recent DALI models by Trapman et al. (2021) that include CO
conversion also find a factor of $\sim 3-5$ difference in the C18O and 13CO
fluxes from earlier DALI models by Miotello et al. (2017).

More fundamentally, the location of the snowline is in itself a function of
the disk structure. Figure 7 shows the comparisons at $r=50$ au between our
vertical hydrostatic equilibrium model (VSHE, top left panel) and models where
the density profile is a fixed gaussian with two choices of parametrization.
The density profile is assumed to be determined by a power law form for the
scale height $h$ and two values typically adopted in the literature ($h/r\sim
0.1$ and $h/r\sim 0.2$) are used. The density distribution for the two
gaussian profiles (top right, $h/r=0.1$; bottom left, $h/r=0.2$) are evidently
different from the pressure equilibrium profile (top left panel of Fig. 7).
Further, the dust temperature transitions at different heights in the disk
depending on the adopted value of the scaleheight. These differences result in
marked changes in $z_{\mathrm{PDR}}$ (due to density) and $z_{\mathrm{ice}}$
(due to dust temperature profile) and therefore very different C18O densities
when compared with our self-consistent models. The bottom right panel of Fig.
7 shows the corresponding range in column density as a function of height
which is seen to be quite sensitive to the adopted parametrization of the
density profile. We thus find that calculating the disk density structure with
reasonable accuracy, in tandem with the dust radiative transfer and in a self-
consistent manner is important for calculating the C18O line luminosities.

Figure 7:  Computed vertical structure as a function of $z/r$ and at $r=50$au
of a disk with $\epsilon=10^{-2}$ and
$\Sigma_{\text{gas}}(r=50\mathrm{au})\sim 2.0$ g.cm-2 obtained with our full
thermo-chemical disk model (a), and a parametric expression which assumes that
the vertical density structure follows a Gaussian profile with $h/r=0.1$ (b)
and $h/r=0.2$ (c). The black solid lines are for the gas density, the blue
dashed lines for the dust temperature and the red dotted lines for the
computed C18O density. The last panel (d) shows the range of vertical C18O
column density as a function of height computed for these different models.
The effect of the conversion of CO into CO2 ice at the surface of the grains
is taken into account for all these models.

The HD emission lines in our model are also systematically lower by a factor
$\sim 2$ as compared to the emission calculated by Trapman et al. (2017) using
the DALI code, and this again may be directly related to differences in disk
physical structure. HD emission also arises from above the $A_{V}=1$ layer
where our densities are lower as described above, but is in addition sensitive
to the column density of warm gas (unlike C18O emission which is only
sensitive to column density). Since the bulk of the disk gas is cold,
differences in the temperature (e.g. mean HD emitting temperature of 45 K vs.
35K) could account for the factor of two difference in emission between the
models. HD emission is further sensitive to other factors like the dust/gas
ratio, which affect the disk physical structure and the detectability of the
line above the continuum. The disk structure in our model is strongly coupled
to disk properties (surface density of gas and dust are the only varied
parameters), but many quantities including the dust scale height are treated
as independent variable parameters in Trapman et al. (2017) making direct
comparisons with our work difficult. We note however that in a test model,
Trapman et al. (2017) also find a factor of 2 decrease in the computed HD
(1-0) emission when the density and temperature structure are iterated a few
times by solving the equation of hydrostatic equilibrium (see their Appendix
C, note they do not run these models to convergence). We therefore conclude
that the main difference can be attributed to their isothermal density
structure and our density/temperature coupling through the hydrostatic
equilibrium condition.

###  4.3 The TW Hya disk and HD vs. CO discrepancies

Most arguments for CO depletion and underabundances of both C and O come from
the copious observational data available for the disk around the nearby star,
TW Hya (e.g., Bergin et al., 2016). This is the only disk for which HD has
been detected with sufficient S/N, and for which there is therefore a
chemistry-independent assessment of gas mass (Bergin et al., 2013). In this
section we compare our results with observations for the TW Hya disk using the
structure determined in Gorti et al. (2011) using the thermo-chemical model of
Gorti & Hollenbach (2004, 2008) which solves for vertical hydrostatic
equilibrium. This model was based on the dust model by Calvet et al. (2002)
and reproduces the spectral energy distribution (SED) and can match the
integrated emission of several surface emission lines including H2 S(1) and
S(2), HCO+ (4-3), CO rotational lines, CO rovibrational lines, OH line
emission, and [Ne II], [Ne III], [O I] fine structure emission and optical
forbidden lines, as well as the HD flux observed with Herschel (Bergin et al.,
2013; McClure et al., 2016) to within a factor of $\sim 2-3$. In Gorti et al.
(2011), no depletion factors for C or O were assumed, and normal ISM-like
elemental abundances along with the standard value of $\epsilon=10^{-2}$ were
used. Using this disk physical structure (densities and temperatures), here we
compute gas grain chemistry as in RG19 and line luminosities for 12CO, 13CO,
C18O and a few other species (Fig. 8) that were not previously considered.

With the notable exception of hydrocarbons, to be discussed later, the the new
chemical model reproduces observed CO isotopologue lines to within a factor of
a few despite using standard values for $\epsilon$ and the elemental C and O
abundances. 12CO and C18O fluxes are very well matched (see Fig. 8), but as in
the Lupus survey comparison there is systematic overestimation of the 13CO
(3-2) and (6-5) emission by a factor of $\sim 2-3$ (note TW Hya is nearly
face-on, so inclination does not explain this discrepancy, see discussion in
§3.1).

HD $1-0$ emission is under-estimated by a factor of $\sim$ 2, as in earlier
work (Bergin et al., 2013) that used the same Gorti et al. (2011) disk
structure, and there is thus an inconsistency of a factor of $\sim 2$ with
respect to the C18O emission. Previous authors (e.g. Favre et al., 2013; Du et
al., 2015; Kama et al., 2016; Trapman et al., 2017; Bergin & Williams, 2017)
have, however, invoked large depletion factors ranging from $10-100$ of carbon
and oxygen to reconcile HD and CO observations in TW Hya. For the various
reasons discussed in §4.2, abandoning parametrization of the disk structure
and adopting vertical hydrostatic equilibrium derived density/temperature
distributions may alleviate the need for such severe depletion of C and O. Our
results nevertheless indicate that CO isotopologue emission and HD emission
could be discrepant, but by a small factor of $\sim 2-3$ at most. To some
extent, the lower HD could be mitigated assuming a slightly enhanced elemental
abundance of D (we use $X(\mathrm{D})=1.6\times 10^{-5}$ compared to
$X(\mathrm{D})\sim 2\times 10^{-5}$ used in e.g Trapman et al., 2017; Calahan
et al., 2021b). Considering dust settling (which was not treated in Gorti et
al., 2011) may also slightly decrease the dust background and increase the
line flux over the continuum. We note that we did not recompute the TW Hya
disk model with settling as this would involve reconciling the spectral energy
distribution with the gas disk structure again and involves considerable
effort. This will be the subject of a forthcoming study. Another possibility
is that the disk mass is higher. However, the system is an old ($\sim$ 8 Myr)
transition disk with evidence for planet formation and possibly embedded
planets (e.g. Teague et al., 2019). Disk masses higher than the $0.06$M⊙
indicated here would moreover imply this evolved disk is gravitationally
unstable with a mass comparable to the central star. We therefore consider
significantly higher gas disk masses to be unlikely, but note that a mass
$\sim 2$ times higher than in our model has been inferred by Powell et al.
(2019) with a different approach using resolved multiwavelength continuum
data.

While our models with ISM elemental abundances of C and O (and C/O $\sim
0.45$) can successfully explain most line emission data (also see Gorti et
al., 2011), including HD, calculated emission from C2H and C3H2 is far lower
than observed. A change in the C/O ratio to a value $>1$ is a straightforward
explanation as this leaves free carbon (not bound in CO) available for
hydrocarbon production as expected from chemical considerations; C/O$>1$ has
thus been invoked in order to explain abundances of hydrocarbons inferred from
observations (Bergin et al., 2016; Kama et al., 2016). As seen in Fig. 8,
increasing the C/O ratio via oxygen depletion (i.e. C/H=O/H=$1.4\times
10^{-4}$) or carbon enhancement (i.e. C/H=O/H=$3.2\times 10^{-4}$) both give
nearly similar results for all molecules and only marginally affect the CO
isotopologue line emission. Assuming C/H=O/H=$1.4\times 10^{-4}$ slightly
increases C18O emission because lowering O decreases free oxygen available to
form water ice and affects the conversion of CO to CO2 ice on grain surfaces
as it also lowers the abundance of the OH radical. This shifts the vertical CO
snowline to lower $z$ and results in an higher column density of CO and thus
higher C18O emission (see §3.1). A further decrease in oxygen to a C/O ratio
of 2, decreases the abundance of C18O throughout the disk due to limited
availability of oxygen and slightly decreases C18O emission. Enhancing carbon
to affect the C/O ratio results in slightly higher C18O luminosities due to an
overall increase in CO isotopologue abundance. 12CO and 13CO lines are
optically thick in this disk and are not affected by changes in the C/O ratio.

Hydrocarbon emission, on the other hand, is drastically affected by changes in
the C/O ratio, with a factor of $\gtrsim 10$ increase in line luminosities due
to the increased availability of carbon. This is due to the fact that the
hydrocarbon abundances needed to produce the observed emission are only a
small fraction of the elemental carbon available (e.g.,
$X(\text{C}_{2}\text{H})\sim 10^{-8}\sim 10^{-4}X(\text{C})$). While C3H2
emission is well explained with C/O=1 (which is independent of the assumption
made on either enhancing carbon or depleting oxygen), C2H emission is still
under-produced by both the C/O $\sim$ 1 models. Increasing C/O to 2 results in
an improvement in matching C2H emission, but significantly over-produces C3H2
emission. We can thus constrain any possible change in the C/O ratio in the TW
Hya disk to $\sim 1-2$, which causes only in a small change in the CO
isotopologue emission lines.

The main conclusion that can be drawn from the models is that lowering the
elemental abundances of C and O by large factors is not required to explain
line emission from the TW Hya disk, and hydrocarbon emission favors a possibly
small change in the C/O ratio; this could moreover be restricted to the
surface of the disk where emission lines originate. Furthermore, unlike that
of carbon chains, 13CO emission is only marginally affected by a change in the
C/O ratio.

We note, however, that spatially resolved line observations are available for
the TW Hya disk, and many of the models invoking elemental depletion Zhang et
al. (2017, 2019); Calahan et al. (2021b) have attempted to recover the radial
profile of CO isotopologue emission with parametric fitting. Data obtained
from the recent ALMA MAPS program (Oberg et al., 2021), and in particular
constraints obtained in favorably inclined disks also allow inferences on the
vertical origin of the emission (e.g. Zhang et al., 2021; Calahan et al.,
2021a; Schwarz et al., 2021). For a more definitive conclusion on the extent
to which the C/O ratio changes in disks we need to self-consistently model the
structure of individual disks in greater detail to simultaneously reproduce
all the infrared line emission and all spatially resolved sub-millimeter line
emission; this will be subject of future work.

Figure 8:  Computed line luminosities for TW Hya and comparison with
observations of HD (Bergin et al., 2013), CO (Schwarz et al., 2016; Calahan et
al., 2021b), C2H (Bergin et al., 2016), C3H2 (Cleeves et al., 2021). The
reference model (C/O=0.45) assumes ISM-like elemental abundances, i.e. C/H =
$1.4\times 10^{-4}$ and O/H = $3.2\times 10^{-4}$.

##  5 Simplified models for estimating gas disk masses using C18O line fluxes

We have thus far argued that C18O lines are reliable tracers of the disk gas
mass; the analysis presented in §4 demonstrated that C18O emission traces a
similar total disk mass fraction as HD, is less sensitive to uncertainties in
disk structure and temperature, and that consistent physico-chemical modeling
indicates discrepancies are limited to a factor $\sim 2$ when compared to HD
for the TW Hya disk. However, our models are inordinately complex and
estimating disk masses from observations by detailed modeling of disks (as for
TW Hya) is impractical. From our model results, we have identified two main
modeling components that need to be considered in order to infer spatial
density distributions and estimate gas masses to within a factor of a few from
C18O emission lines (and 13CO lines):

  * -

A physical model that provides a reasonable estimate of the vertical density
structure at the $A_{V}\sim 1-3$ layer near the water condensation front.

  * -

A minimal chemical network including the chemistry of CO isotopologues,
selective photodissociation, freeze-out and conversion of CO into CO2 ice at
the surface of grains.

###  5.1 Simplified vertical density structure of the disk

Determining the disk physical structure requires solutions to the equations of
vertical pressure balance coupled to the equations of radiative transfer.
Pressure is determined by gas temperature and density, and the gas temperature
is further coupled to the chemistry and dust properties, making consistent
disk solutions a formidable computational effort. Therefore, the problem
becomes much more tractable by decoupling the density structure calculation
from the chemistry. Models often simplify the procedure by assuming a
parametric Gaussian (gas assumed isothermal) vertical density profile that is
kept fixed, and the gas temperature and chemistry are iterated to convergence.
Very often, when modeling individual disks, the procedure adopted is justified
by finding a dust temperature that can reproduce observed continuum spectral
energy distributions (SEDs). However, the dust temperature is not isothermal
with height and transitions from a cooler midplane to a warmer surface near
the $A_{V}=1$ layer (see several older works, for e.g., Calvet et al., 1991;
Chiang & Goldreich, 1997; D’Alessio et al., 1998; Dullemond et al., 2001). As
discussed earlier (also see RG19), most emission lines originate above
$A_{V}=1$ with excitation determined by the temperature at the surface, while
the midplane temperature characterizes the regions where most of the gas mass
resides and sets the gas vertical density structure. Therefore, characterizing
the density profile with a single gaussian characteristic of the surface
emission will lead to an underestimate of the density in the cooler midplane
regions (see Fig. 7). We therefore advocate against parametrizations based on
the notion of a single scaleheight.

For modeling C18O isotopologue emission, in particular, the density structure
is more critical than the temperature structure (see §4.1) and we find that a
dust disk model can reasonably approximate the full solution (see e.g.,
McClure et al., 2016; Molyarova et al., 2017b, for a similar approach) and
provide a useful alternative approach. We advocate a density structure based
on the dust temperature profile that is consistent with pressure equilibrium.
We use the publicly available Monte Carlo radiative transfer code RADMC-3D
(Dullemond et al., 2012) and follow a very simple procedure to illustrate this
method. For an input surface density distribution, we set an initial guess
density distribution with height and calculate the temperature with RADMC-3D.
In our simple setup, we use a size averaged dust opacity which yields a single
effective (radiative) dust temperature and equate this to the gas temperature.
The density distribution is re-calculated according to the vertical pressure
gradient and re-normalized to the surface density (see Eq. 1). The new density
distribution is used again as input to RADMC-3D to get a new temperature
distribution and this process is repeated to convergence. We note that dust
settling is not taken into account in this procedure although we find that the
resulting change to the dust structure does impact the C18O by a small factor.

###  5.2 Simplified CO chemistry

The disk physical structure (i.e. density and temperature as function of
$(r,z)$) obtained as described above can then be used to compute the CO
chemistry. The minimal chemical network for estimating disk gas masses from CO
isotopologue observations needs to include H2 and CO isotopologue gas phase
chemistry, self-shielding, selective photodissociation, CO freeze-out and
subsequent conversion of CO into CO2 ice at the surface of grains. For H2 and
CO chemistry, we find that the reduced network proposed by Gong et al. (2017)
for hydrogen and carbon chemistry gives good results when compared with our
full network. It includes 18 chemical species coupled through $\sim 50$
reactions, including photoreactions as well as self- and mutual-shielding for
H2 (Draine & Bertoldi, 1996) and CO (Visser et al., 2009). This network then
needs to be extended to take into account the chemistry of 13C and 18O and in
particular reactions listed in Roueff et al. (2015) for 13C and Loison et al.
(2019) for 18O. For all other reactions, the chemical network can simply be
duplicated assuming similar rate constants and statistical branching ratios.
Isotope selective photodissociation can reduce C18O emission by up to a factor
of $\sim 10$ (Miotello et al., 2014) and therefore needs to be included;
shielding functions are provided in Visser et al. (2009). Finally, we find
that $\sim 10$ reactions are sufficient for CO grain surface chemistry. The
set of reactions is given in Appendix A and takes into account CO freeze-out
and its subsequent conversion into CO2 ice as well as water ice formation and
ice photochemistry, which are both important to determine the location of the
vertical CO snowline.

###  5.3 C18O emission from simplified models

In Fig. 9 we show a comparison of the C18O emission computed using the full
thermo-chemical model and the simple model described above. For most of the
disks explored here, the use of a simple dust model leads to $\lesssim 2-3$
difference in the computed C18O (2-1) and (3-2) emission compared to our full
thermo-chemical model, and provides a computationally efficient and reasonably
accurate procedure to retrieve gas disk masses from observational data. The
reason that this procedure works reasonably well is because C18O emission
mainly arises from regions where dust and gas are collisionally coupled and
where gas temperature deviates only slightly from that of the dust. The
RADMC-3D-derived model is found to over-estimate the C18O emission slightly,
due to differences in density caused by differences in the true gas
temperature from the dust mean radiative value and consequently the column
density of C18O in the disk.

We emphasize that while the above procedure may be a good approximation to
compute emission lines from layers close to the water snowline, it fails at
greater heights where gas and dust decouple and their temperatures deviate,
and is not appropriate for modeling 12CO, 13CO, as well as near and mid-
infrared emission lines. The disk structure also depends on the dust
properties (grain size distribution, opacity, settling) and uncertainties in
these will consequently affect derived gas disk masses as well. Nevertheless,
we expect the above method to be a good approximation for gas masses.

|
---|---
Figure 9: Comparison between C18O (2-1) and (3-2) emission computed with the
full thermo-chemical model (solid lines) and the simple RADMC-3D-based dust
disk model (dashed lines), for different dust/gas ratios. The filled region
denotes a factor of 2 deviation from the emission obtained with the detailed
model. For most disks, results obtained from the simple model lie close or
within this region.

##  6 Conclusions

Using self-consistent thermo-chemical gas and dust disk models coupled with
detailed chemical models including dust settling and gas-grain chemistry, we
reassess the potential of using CO isotopologue emission for tracing gas
masses of protoplanetary disks. A parameter survey was conducted by varying
gas disk masses, dust/gas mass ratio and disk outer radius to study the
effects on the computed CO isotopologue and HD emission. Our main conclusion
is that, using appropriate disk models, C18O line emission can be used to
infer gas disk masses to accuracies comparable to what is achieved by disk
dust mass estimations. C18O lines are readily accessible with ALMA, can
further be spatially resolved to retrieve the surface density distribution,
and multiple rotational transitions can be used to constrain the gas
temperature. In addition to its accessibility, C18O is moreover not sensitive
to the typically unknown gas temperature structure of the disk, and traces a
similar disk mass fraction compared to HD, making it a better indicator of
disk mass. Our results can be summarized as follows:

  1. 1.

The disk physical structure has a significant impact on the computed line
emission of CO isotopologues and HD. In particular, in order to determine gas
masses using C18O emission lines it is necessary to make a reasonable estimate
of the vertical density structure just above the CO freeze-out layer where
most of the emission originates. We argue that parametrized density profiles
are inadequate and pressure equilibrium needs to be considered in disk
structure models.

  2. 2.

We find that HD and C18O both trace only a very small fraction of the total
gas mass ranging from $\sim 10^{-6}-10^{-3}$ for all the various disk models
considered here. This is because both HD and C18O emission originate from
regions above the disk midplane where most of the mass resides.

  3. 3.

Based on the above finding, we conclude that C18O is a overall a good tracer
of the gas disk mass for typical disk conditions. HD and C18O emission both
depend on the gas density, but contrary to HD emission which shows an
exponential dependence to the temperature of the gas, C18O emission is
relatively insensitive to the temperature in the emitting layer.

  4. 4.

An important finding of our study is that we can explain CO isotopologue line
emission from surveys using ISM-like elemental abundances, i.e. C/H =
$1.4\times 10^{-4}$ and O/H = $3.2\times 10^{-4}$, and a normal interstellar
dust/gas ratio, i.e. $\epsilon=10^{-2}$. We reach the same conclusion for the
TW Hya disk, where modeled HD and C18O emission are only discrepant by a
factor of $\sim 2$. This result substantially differs from those obtained by
other studies where strong depletion factors of C and O are usually invoked.
We find that the assumption of an isothermal profile for the density can lead
to large deviations in the density at C18O emitting heights and therefore also
the computed column densities and hence line emission.

  5. 5.

In agreement with previous studies, we also find that explaining the
observation of species like small hydrocarbons in TW Hya seems to require a
factor of $\sim 2$ enhancement in elemental carbon or a factor of $\sim 2$
reduction in elemental oxygen as compared to the ISM, with their ratios
ranging between C/O$\sim 1-2$. We also note that hydrocarbon emission does not
provide any information on a potential global depletion of carbon and oxygen
but only on their relative abundance (i.e. C/O) ratio.

  6. 6.

Given that detailed models of disks as those presented here are complex, we
provide a simplified approach to computing C18O line fluxes and determine gas
disk masses using a dust-derived disk structure.

We are thankful to the anonymous referee for a careful and thorough report
which improved the quality of this manuscript. MR and UG acknowledge support
by the National Aeronautics and Space Administration through the NASA
Astrobiology Institute under Cooperative Agreement Notice NNH13ZDA017C issued
through the Science Mission Directorate. UG also acknowledges support from an
award through NASA XRP 80NSSC20K0273.

##  Appendix A Reduced network for CO ice chemistry and snowline location

The location of the CO snowline can be obtained by solving the set of
reactions:

$\displaystyle{\rm H/H_{2}}$ $\displaystyle\rightarrow$ $\displaystyle{\rm
sH/sH_{2}}~{}(k_{\text{acc}})$ $\displaystyle{\rm sH/sH_{2}}$
$\displaystyle\rightarrow$ $\displaystyle{\rm H/H_{2}}~{}(k_{\text{des}})$
$\displaystyle{\rm O}$ $\displaystyle\rightarrow$ $\displaystyle{\rm
sO}~{}(k_{\text{acc}})$ $\displaystyle{\rm sO+sH}$ $\displaystyle\rightarrow$
$\displaystyle{\rm sOH}~{}(k_{\text{reac}})$ $\displaystyle{\rm
sOH+sH/sH_{2}}$ $\displaystyle\rightarrow$ $\displaystyle{\rm
sH_{2}O}~{}(k_{\text{reac}})$ $\displaystyle{\rm sO}$
$\displaystyle\rightarrow$ $\displaystyle{\rm O}~{}(k_{\text{des}})$
$\displaystyle{\rm sH_{2}O}+\gamma$ $\displaystyle\rightarrow$
$\displaystyle{\rm sOH+sH}~{}(k_{\text{diss}})$ $\displaystyle{\rm CO}$
$\displaystyle\rightarrow$ $\displaystyle{\rm sCO}~{}(k_{\text{acc}})$
$\displaystyle{\rm sCO}$ $\displaystyle\rightarrow$ $\displaystyle{\rm
CO}~{}(k_{\text{des}})$ $\displaystyle{\rm sCO+sOH}$
$\displaystyle\rightarrow$ $\displaystyle{\rm
sCO_{2}+sH}~{}(k_{\text{reac}};\Delta E=150\text{K})$ $\displaystyle{\rm
sCO_{2}}$ $\displaystyle\rightarrow$ $\displaystyle{\rm CO}+{\rm
O}~{}(k_{\text{des}})$

In this reduced network, we neglect the desorption of sOH and sH2O from the
grains and assume that CO2 is instantly photodissociated back to CO upon
desorption (last reaction). In practice, this simplified network does not
allow an analytical solution and must be solved numerically. In this set of
reactions, $k_{\text{acc}}$ is the accretion rate and is given by

$k_{\text{acc}}(i)=v_{i}\langle\sigma_{d}n_{d}\rangle$ (A1)

where $v_{i}$ is the thermal velocity of gas species $i$ and
$\langle\sigma_{d}n_{d}\rangle$ the product of the dust cross-sectional area
and the dust density averaged over the grain size distribution. Note that this
expression assumes that the sticking coefficient is equal to 1.
$k_{\text{reac}}$ is the reaction rate at the surface and is computed by

$k_{\text{reac}}(ij)=\kappa_{ij}\frac{k_{\text{hop}}(i)+k_{\text{hop}}(j)}{4\langle\sigma_{d}n_{d}\rangle
N_{\text{sites}}}$ (A2)

where $\kappa_{ij}$ is the probability that the reaction occurs,
$k_{\text{hop}}(i)$ the thermal hopping rate of species $i$ and
$N_{\text{sites}}\sim 10^{15}$ cm-2 the number of adsorption sites per cm2.
The thermal hopping rate is given by

$k_{\text{hop}}(i)=\nu_{i}\exp\big{(}-E_{\text{diff}}(i)/T_{\text{d}}\big{)}$
(A3)

where $T_{d}$ is the dust temperature,
$\nu_{i}=\sqrt{2N_{\text{sites}}kE_{\text{bind}}(i)/\pi^{2}m_{i}}$ is the
vibrational frequency of the species $i$, $E_{\text{bind}}$ its binding
energy, $m_{i}$ its mass and $E_{\text{diff}}(i)$ the barrier against
diffusion. $\kappa_{ij}$ is equal to 1 for all reactions that don’t have a
barrier. For reactions with barrier $\kappa_{ij}$ is computed as the result of
the competition among reaction and diffusion (Garrod & Pauly, 2011; Ruaud et
al., 2016)

$\kappa_{ij}=\frac{\max(\nu_{i},\nu_{j})\kappa_{ij,0}}{\max(\nu_{i},\nu_{j})\kappa_{ij,0}+k_{\text{hop},i}+k_{\text{hop},j}}$
(A4)

where $\kappa_{ij,0}$ is expressed as the $\exp(-\Delta E_{ij}/T_{\text{d}})$
with $\Delta E_{ij}$ the activation barrier of the reaction. For reactions
involving H and H2 we use the quantum mechanical probability for tunneling
through a rectangular barrier of thickness $a=1$Å;
$\kappa_{ij,0}=\exp[-2(a/\hbar)(2\mu\Delta E_{ij})^{1/2}]$, with $\mu$ the
reduced mass of the system. Finally, $k_{\text{des}}$ is the desorption rate
by both thermal and non-thermal processes (i.e. photodesorption in this case)

$k_{\text{des}}(i)=\Big{[}\nu_{i}\exp\big{(}-E_{\text{bind}}(i)/T_{\text{d}}\big{)}+\frac{Y_{\text{pd}}(i)F_{\text{FUV}}}{4N_{\text{sites}}N_{\text{act}}}\Big{]}\eta_{\text{des}}$
(A5)

where $E_{\text{bind}}$ is the binding energy of species $i$, $Y_{\text{pd}}$
its photodesorption yield, $F_{\text{FUV}}$ is the incident FUV photon flux on
the grain, $N_{\text{act}}$ the number of active surface layers (taken to be
equal to 2 in our model) and $\eta_{\text{des}}$ the fraction of molecular
species in the top desorbable surface layers and given by

$\eta^{\text{des}}=\begin{cases}1,&\text{if we have less than $N_{\text{lay}}$
monolayers of ice}\\\ \frac{4\langle\sigma_{d}n_{d}\rangle
N_{\text{sites}}N_{\text{act}}}{n_{\text{ice}}},&\text{otherwise.}\end{cases}$
(A6)

Finally, $k_{\text{diss}}$ is the photodissociation rate by both direct FUV
photons and CR generated photons and is given by

$k_{\text{diss}}(i)=G_{0}\alpha_{i}\exp\big{(}-\beta_{i}A_{V}\big{)}+\gamma_{i}\zeta_{\text{H}_{2}}f_{\text{H}_{2}}$
(A7)

where $G_{0}$ is the integrated FUV flux with respect to that of the ISRF in
units of the Draine field in the 6-13.6 eV range, $\alpha$ is the unattenuated
photorate, $\gamma$ the cosmic ray induced photodissociation efficiency,
$\zeta_{\text{H}_{2}}$ the cosmic ray ionization rate of H2 and
$f_{\text{H}_{2}}=n(\text{H}_{2})/[n(\text{H})+2n(\text{H}_{2})]$ the
molecular fraction. For water ice, we use $\alpha=10^{-9}$ s-1, $\beta=1.8$
and $\gamma=10^{3}$. Binding energies used for each of the species included in
the reduced network are given in Table 2. For diffusion energy barriers we
assume a fixed fraction of the binding energy of $0.4\times E_{\text{bind}}$.
For more details on the different parameters used here see Ruaud et al. (2016)
and Ruaud & Gorti (2019).

Table 2: Binding energies and diffusion energy barriers used for each of the species included in the reduced network. Species | Binding energy (K) | Reference
---|---|---
sH | 650 | 1
sH2 | 440 | 2
sO | 800 | 3
sCO | 1150 | 3
sCO2 | 2575 | 3
sOH | 2850 | 3
sH2O | 5700 | 3

## References

  * Aikawa et al. (2015) Aikawa, Y., Furuya, K., Nomura, H., & Qi, C. 2015, ApJ, 807, 120
  * Al-Halabi & van Dishoeck (2007) Al-Halabi, A., & van Dishoeck, E. F. 2007, MNRAS, 382, 1648
  * Alcalá et al. (2017) Alcalá, J. M., Manara, C. F., Natta, A., et al. 2017, A&A, 600, A20
  * Amiaud et al. (2007) Amiaud, L., Dulieu, F., Fillion, J. H., Momeni, A., & Lemaire, J. L. 2007, J. Chem. Phys., 127, 144709
  * Ansdell et al. (2017) Ansdell, M., Williams, J. P., Manara, C. F., et al. 2017, AJ, 153, 240
  * Ansdell et al. (2016) Ansdell, M., Williams, J. P., van der Marel, N., et al. 2016, ApJ, 828, 46
  * Ballering & Eisner (2019) Ballering, N. P., & Eisner, J. A. 2019, AJ, 157, 144
  * Bergin et al. (2014) Bergin, E. A., Cleeves, L. I., Crockett, N., & Blake, G. A. 2014, Faraday Discussions, 168, 61
  * Bergin et al. (2016) Bergin, E. A., Du, F., Cleeves, L. I., et al. 2016, ApJ, 831, 101
  * Bergin & Williams (2017) Bergin, E. A., & Williams, J. P. 2017, The Determination of Protoplanetary Disk Masses, ed. M. Pessah & O. Gressel, Vol. 445, 1
  * Bergin et al. (2013) Bergin, E. A., Cleeves, L. I., Gorti, U., et al. 2013, Nature, 493, 644
  * Bosman et al. (2018) Bosman, A. D., Walsh, C., & van Dishoeck, E. F. 2018, A&A, 618, A182
  * Bruderer et al. (2012) Bruderer, S., van Dishoeck, E. F., Doty, S. D., & Herczeg, G. J. 2012, A&A, 541, A91
  * Calahan et al. (2021a) Calahan, J. K., Bergin, E. A., Zhang, K., et al. 2021a, arXiv e-prints, arXiv:2109.06202
  * Calahan et al. (2021b) Calahan, J. K., Bergin, E., Zhang, K., et al. 2021b, ApJ, 908, 8
  * Calvet et al. (2002) Calvet, N., D’Alessio, P., Hartmann, L., et al. 2002, ApJ, 568, 1008
  * Calvet et al. (1991) Calvet, N., Patino, A., Magris, G. C., & D’Alessio, P. 1991, ApJ, 380, 617
  * Chiang & Goldreich (1997) Chiang, E. I., & Goldreich, P. 1997, ApJ, 490, 368
  * Cleeves et al. (2021) Cleeves, L. I., Loomis, R. A., Teague, R., et al. 2021, arXiv e-prints, arXiv:2102.09577
  * D’Alessio et al. (1998) D’Alessio, P., Cantö, J., Calvet, N., & Lizano, S. 1998, ApJ, 500, 411
  * Draine & Bertoldi (1996) Draine, B. T., & Bertoldi, F. 1996, ApJ, 468, 269
  * Du et al. (2015) Du, F., Bergin, E. A., & Hogerheijde, M. R. 2015, ApJ, 807, L32
  * Dullemond et al. (2001) Dullemond, C. P., Dominik, C., & Natta, A. 2001, ApJ, 560, 957
  * Dullemond et al. (2012) Dullemond, C. P., Juhasz, A., Pohl, A., et al. 2012, RADMC-3D: A multi-purpose radiative transfer tool, , , ascl:1202.015
  * Eistrup et al. (2016) Eistrup, C., Walsh, C., & van Dishoeck, E. F. 2016, A&A, 595, A83
  * Favre et al. (2013) Favre, C., Cleeves, L. I., Bergin, E. A., Qi, C., & Blake, G. A. 2013, ApJ, 776, L38
  * Furuya & Aikawa (2014) Furuya, K., & Aikawa, Y. 2014, ApJ, 790, 97
  * Garrod & Herbst (2006) Garrod, R. T., & Herbst, E. 2006, A&A, 457, 927
  * Garrod & Pauly (2011) Garrod, R. T., & Pauly, T. 2011, ApJ, 735, 15
  * Gong et al. (2017) Gong, M., Ostriker, E. C., & Wolfire, M. G. 2017, ApJ, 843, 38
  * Gorti & Hollenbach (2004) Gorti, U., & Hollenbach, D. 2004, ApJ, 613, 424
  * Gorti & Hollenbach (2008) —. 2008, ApJ, 683, 287
  * Gorti et al. (2015) Gorti, U., Hollenbach, D., & Dullemond, C. P. 2015, ApJ, 804, 29
  * Gorti et al. (2011) Gorti, U., Hollenbach, D., Najita, J., & Pascucci, I. 2011, ApJ, 735, 90
  * Hildebrand (1983) Hildebrand, R. H. 1983, QJRAS, 24, 267
  * Hollenbach & McKee (1979) Hollenbach, D., & McKee, C. F. 1979, ApJS, 41, 555
  * Kama et al. (2016) Kama, M., Bruderer, S., Carney, M., et al. 2016, A&A, 588, A108
  * Kama et al. (2020) Kama, M., Trapman, L., Fedele, D., et al. 2020, A&A, 634, A88
  * Krijt et al. (2020) Krijt, S., Bosman, A. D., Zhang, K., et al. 2020, ApJ, 899, 134
  * Krijt et al. (2016) Krijt, S., Ciesla, F. J., & Bergin, E. A. 2016, ApJ, 833, 285
  * Krijt et al. (2018) Krijt, S., Schwarz, K. R., Bergin, E. A., & Ciesla, F. J. 2018, ApJ, 864, 78
  * Loison et al. (2014) Loison, J.-C., Wakelam, V., Hickson, K. M., Bergeat, A., & Mereau, R. 2014, MNRAS, 437, 930
  * Loison et al. (2019) Loison, J.-C., Wakelam, V., Gratier, P., et al. 2019, MNRAS, 485, 5777
  * Long et al. (2017) Long, F., Herczeg, G. J., Pascucci, I., et al. 2017, ApJ, 844, 99
  * Macías et al. (2021) Macías, E., Guerra-Alvarado, O., Carrasco-González, C., et al. 2021, A&A, 648, A33
  * Manara et al. (2016) Manara, C. F., Rosotti, G., Testi, L., et al. 2016, A&A, 591, L3
  * McClure (2019) McClure, M. K. 2019, A&A, 632, A32
  * McClure et al. (2016) McClure, M. K., Bergin, E. A., Cleeves, L. I., et al. 2016, ApJ, 831, 167
  * Miotello et al. (2014) Miotello, A., Bruderer, S., & van Dishoeck, E. F. 2014, A&A, 572, A96
  * Miotello et al. (2016) Miotello, A., van Dishoeck, E. F., Kama, M., & Bruderer, S. 2016, A&A, 594, A85
  * Miotello et al. (2017) Miotello, A., van Dishoeck, E. F., Williams, J. P., et al. 2017, A&A, 599, A113
  * Molyarova et al. (2017a) Molyarova, T., Akimkin, V., Semenov, D., et al. 2017a, ApJ, 849, 130
  * Molyarova et al. (2017b) —. 2017b, ApJ, 849, 130
  * Öberg & Bergin (2021) Öberg, K. I., & Bergin, E. A. 2021, Phys. Rep., 893, 1
  * Oberg et al. (2021) Oberg, K. I., Guzman, V. V., Walsh, C., et al. 2021, arXiv e-prints, arXiv:2109.06268
  * Padovani et al. (2018) Padovani, M., Ivlev, A. V., Galli, D., & Caselli, P. 2018, A&A, 614, A111
  * Pascucci et al. (2016) Pascucci, I., Testi, L., Herczeg, G. J., et al. 2016, ApJ, 831, 125
  * Powell et al. (2019) Powell, D., Murray-Clay, R., Pérez, L. M., Schlichting, H. E., & Rosenthal, M. 2019, ApJ, 878, 116
  * Reboussin et al. (2015) Reboussin, L., Wakelam, V., Guilloteau, S., Hersant, F., & Dutrey, A. 2015, A&A, 579, A82
  * Ribas et al. (2020) Ribas, Á., Espaillat, C. C., Macías, E., & Sarro, L. M. 2020, A&A, 642, A171
  * Roueff et al. (2015) Roueff, E., Loison, J. C., & Hickson, K. M. 2015, A&A, 576, A99
  * Ruaud & Gorti (2019) Ruaud, M., & Gorti, U. 2019, ApJ, 885, 146
  * Ruaud et al. (2016) Ruaud, M., Wakelam, V., & Hersant, F. 2016, MNRAS, 459, 3756
  * Schwarz et al. (2016) Schwarz, K. R., Bergin, E. A., Cleeves, L. I., et al. 2016, ApJ, 823, 91
  * Schwarz et al. (2021) Schwarz, K. R., Calahan, J. K., Zhang, K., et al. 2021, arXiv e-prints, arXiv:2109.06228
  * Teague et al. (2019) Teague, R., Bae, J., Huang, J., & Bergin, E. A. 2019, ApJ, 884, L56
  * Tielens & Hollenbach (1985) Tielens, A. G. G. M., & Hollenbach, D. 1985, ApJ, 291, 722
  * Trapman et al. (2021) Trapman, L., Bosman, A. D., Rosotti, G., Hogerheijde, M. R., & van Dishoeck, E. F. 2021, arXiv e-prints, arXiv:2103.05654
  * Trapman et al. (2017) Trapman, L., Miotello, A., Kama, M., van Dishoeck, E. F., & Bruderer, S. 2017, A&A, 605, A69
  * Tychoniec et al. (2020) Tychoniec, Ł., Manara, C. F., Rosotti, G. P., et al. 2020, A&A, 640, A19
  * Visser et al. (2009) Visser, R., van Dishoeck, E. F., & Black, J. H. 2009, A&A, 503, 323
  * Williams & Best (2014) Williams, J. P., & Best, W. M. J. 2014, ApJ, 788, 59
  * Woitke et al. (2016) Woitke, P., Min, M., Pinte, C., et al. 2016, A&A, 586, A103
  * Xu et al. (2017) Xu, R., Bai, X.-N., & Öberg, K. 2017, ApJ, 835, 162
  * Zhang et al. (2017) Zhang, K., Bergin, E. A., Blake, G. A., Cleeves, L. I., & Schwarz, K. R. 2017, Nature Astronomy, 1, 0130
  * Zhang et al. (2019) Zhang, K., Bergin, E. A., Schwarz, K., Krijt, S., & Ciesla, F. 2019, ApJ, 883, 98
  * Zhang et al. (2021) Zhang, K., Booth, A. S., Law, C. J., et al. 2021, arXiv e-prints, arXiv:2109.06233