#  DiskMINT: A Tool to Estimate Disk Masses with CO Isotopologues

Dingshan Deng Lunar and Planetary Laboratory, the University of Arizona,
Tucson, AZ 85721, USA Dingshan Deng<EMAIL_ADDRESS>Maxime Ruaud NASA
Ames Research Center, Moffett Field, CA 94035, USA  Carl Sagan Center, SETI
Institute, Mountain View, CA 94043, USA  Uma Gorti NASA Ames Research Center,
Moffett Field, CA 94035, USA  Carl Sagan Center, SETI Institute, Mountain
View, CA 94043, USA  Ilaria Pascucci Lunar and Planetary Laboratory, the
University of Arizona, Tucson, AZ 85721, USA

###### Abstract

CO is one of the most abundant molecules in protoplanetary disks, and
optically thin emission from its isotopologues has been detected in many of
them. However, several past works have argued that reproducing the relatively
low emission of CO isotopologues requires a very low disk mass or significant
CO depletion. Here, we present a Python code, DiskMINT, which includes gas
density and temperature structures that are both consistent with the thermal
pressure gradient, isotope-selective chemistry, and conversion of CO into
$\mathrm{CO_{2}}$ ice on grain-surfaces. The code generates a self-consistent
disk structure, where the gas disk distribution is obtained from a Spectral
Energy Distribution (SED)-derived dust disk structure with multiple grain
sizes. We use DiskMINT to study the disk of RU Lup, a high-accreting star
whose disk was previously inferred to have a gas mass of only $\sim 1.5\times
10^{-3}\,M_{\odot}$ and gas-to-dust mass ratio of $\sim 4$. Our best-fit model
to the long-wavelength continuum emission can explain the total
$\mathrm{C^{18}O}$ luminosity as well as the $\mathrm{C^{18}O}$ velocity and
radial intensity profiles, and obtains a gas mass of $\sim 1.2\times
10^{-2}\,M_{\odot}$, an order of magnitude higher than previous results. A
disk model with parametric Gaussian vertical distribution that better matches
the IR-SED can also explain the observables above with a similarly high gas
mass $\sim 2.1\times 10^{-2}\,M_{\odot}$. We confirm the conclusions of Ruaud
et al. (2022) that optically thin $\mathrm{C^{18}O}$ rotational lines provide
reasonable estimates of the disk mass and can therefore be used as gas disk
tracers.

Protoplanetary disks(1300); Astrochemistry(75); Chemical abundances(224); CO
line emission(262); Planet formation(1241)

##  1 Introduction

Disks of gas and dust around young stars (hereafter, protoplanetary disks) are
the sites of planet formation, and their mass is fundamental to understanding
when and how planets and small bodies form. While the gas content sets limits
on the potential masses of forming giant planets, the dust mass constrains the
masses and formation times for the cores of gaseous planets and terrestrial
planets. The gas-to-dust mass ratio ($\Delta_{\mathbf{gd}}$), moreover, is an
indicator of the relative rates of planet formation and gas disk dispersal and
indicates the stage of disk evolution and planet formation (e.g., Miotello et
al. 2022 for a recent review). Ideally, independent and reliable dust and gas
mass estimations are needed to infer the disk physics and evolution, but
measuring both masses is complicated and challenging.

Dust masses ($M_{\mathrm{dust}}$) are estimated by the dust thermal emission
at (sub)millimeter wavelengths, which is sensitive to particles with sizes
$\lesssim 1\,\mathrm{cm}$ and is mostly optically thin (e.g., Ansdell et al.
2016; Pascucci et al. 2016). However, $M_{\mathrm{dust}}$ estimates rely on
the dust opacity $\kappa_{\nu}$ which depends on the composition of dust
grains and their size distribution. Therefore, estimates of
$M_{\mathrm{dust}}$ from a single flux measurement strongly depend on the
assumptions made on the dust properties (e.g., Miotello et al., 2022).
Improved estimates of $M_{\mathrm{dust}}$ can be made by fitting the spectral
energy distribution (SED) at long wavelengths ($\gtrsim 100\,\micron$) where
the emission is typically optically thin (e.g., Woitke et al., 2019).

Gas masses ($M_{\mathrm{gas}}$) are more difficult to estimate since there are
very few optically thin gas emission lines that may trace the disk mass
reservoir. $\mathrm{H_{2}}$ is the most abundant molecule in the gas phase in
the disk, but its emission is faint. This is because $\mathrm{H_{2}}$ is a
light, homonuclear molecule with no permanent dipole moment and hence has only
transitions at high energy levels ($E_{u}\sim$ few 100$-$1000K), while the
majority of the gas in the disk around T-Tauri stars is far colder ($\sim
30\,\mathrm{K}$). The less abundant isotopologue $\mathrm{HD}$ is favored to
measure $M_{\mathrm{gas}}$, although it also traces relatively warm gas
(needed to excite the first rotational level of $\mathrm{HD}$ at $E_{u}\sim
128\,\mathrm{K}$), and therefore has some limitations on its suitability as a
mass tracer (Trapman et al., 2017; Ruaud et al., 2022).

Carbon monoxide ($\mathrm{CO}$) is the most abundant molecule after
$\mathrm{H_{2}}$ and is co-spatially distributed with $\mathrm{H_{2}}$ at the
disk surface. In the disk mid-plane, $\mathrm{CO}$ freezes out on the dust
grain surface (when $T_{\mathrm{dust}}\lesssim 20\,\mathrm{K}$) where it can
be processed into more refractory ices. With its high detectability at
(sub)millimeter wavelengths in disks, $\mathrm{CO}$ and its isotopologues have
long been considered among the best tracers of gas disk mass. However, recent
Atacama Large Millimeter/submillimeter Array (ALMA) observations of Class-II
disks have cast doubts about its ability as a mass tracer because model-
predicted line emissions of $\mathrm{CO}$ and its isotopologues are higher
than observed even after accounting for the fact that $\mathrm{CO}$ freezes-
out in the mid-plane (e.g., Ansdell et al., 2016; Miotello et al., 2017; Long
et al., 2017). This raises questions as to whether CO chemical abundances in
disks differ from that in the interstellar medium (ISM) or whether the disk
gas masses are low. Furthermore, the $\mathrm{CO}$-based $M_{\mathrm{gas}}$
were smaller by $\sim 1-2$ orders of magnitude compared with the
$\mathrm{HD}$-based values for the few disks where $\mathrm{HD}$ has also been
detected (e.g., Bergin et al., 2013; McClure et al., 2016; Trapman et al.,
2017). Thus, some works have argued for higher gas masses but large-scale
depletion of $\mathrm{CO}$ due to dynamical processes that sequester
$\mathrm{CO}$ into forming planetesimals and proto-planets (e.g., Bergin &
Williams, 2017; Bosman & Banzatti, 2019; Sturm et al., 2022).

A different solution was proposed recently by Ruaud et al. (2022) (hereafter
RGH22), who argued that by including (a) the density distribution given by
self-consistent vertical hydrostatic pressure equilibrium, (b) isotopologue-
selective chemistry, and (c) grain-surface chemistry where $\mathrm{CO}$ to
$\mathrm{CO_{2}}$ conversion is a key reaction, the apparent discrepancy
between the $\mathrm{HD}$ and $\mathrm{C^{18}O}$ derived masses can be
resolved. They concluded that CO chemistry in disks is in fact similar to that
in the ISM and that the optically thin lines from $\mathrm{C^{18}O}$ can be
used as a gas mass tracer. Although they could retrieve typical
$\mathrm{C^{18}O}$ fluxes observed for the Lupus sample, they did not consider
individual disks in detail or compare the profile and radial distribution of
the line emission.

In this work, we develop a tool to estimate the disk mass: DiskMINT (Disk
Model for INdividual Targets). It uses the dust temperature-based approach
suggested in RGH22: generating a self-consistent gas disk structure on top of
a SED-derived dust disk. It also uses a reduced chemical network that properly
captures the conversion of $\mathrm{CO}$ into $\mathrm{CO_{2}}$ ice. The tool
is tested in considerable detail for the Class II source RU Lup. We select RU
Lup because this disk has been previously inferred (Miotello et al., 2017) to
have a low gas mass of $M_{\mathrm{gas}}\sim 1.5\times 10^{-3}\,M_{\odot}$
with $\Delta_{\mathrm{gd}}\sim 4$ which is at odds with the large mass
accretion rate onto the star (Alcalá et al., 2017) and the large disk size
(Huang et al., 2018, 2020).

The paper is organized as follows. First, we describe the modeling procedure
in Section 2. Then, we summarize the stellar parameters, observational data,
and model setup for RU Lup in Section 3, followed by the results and
discussion in Section 4. We present our summary and outlook in Section 5.

##  2 Model Description

DiskMINT is a dust temperature-based disk model, and uses the recommendations
made by RGH22. From their analysis using a full thermo-chemical model that
includes isotope-selective photodissociation and 3-phase grain-surface
chemistry, RGH22 identified two main components that can be used to construct
a simplified model to accurately simulate $\mathrm{C^{18}O}$ emission. The two
components are: (a) a self-consistent disk physical structure, based on the
dust temperature $T_{\mathrm{d}}$ and imposing vertical hydrostatic pressure
equilibrium (hereafter VHSE) to calculate densities consistent with this
vertical temperature; and (b) a reduced chemical network that includes
isotope-selective photodissociation and grain-surface chemistry that accounts
for conversion of $\mathrm{CO}$ into $\mathrm{CO_{2}}$ ice (see Appendix A of
RGH22). Since $\mathrm{C^{18}O}$ traces the vertical layer where gas
temperature $T_{\mathrm{g}}$ is still very similar to $T_{\mathrm{d}}$, RGH22
found that a simplified dust disk structure model (which does not consider the
self-consistent gas temperature $T_{\mathrm{g}}$ computed from full thermal
equilibrium) can be used to estimate the $\mathrm{C^{18}O}$ emission. As such,
we build DiskMINT based on this simplified model.

The overall method adopted in our analysis is summarized in the flow chart
shown in Figure 1. Two main steps are involved in obtaining a self-consistent
disk model that fits the $\mathrm{C^{18}O}$ line and continuum data. The goal
of Step 1 is to find a density structure — based on the dust temperature
profile — that is self-consistent with pressure equilibrium and fits the SED.
This is achieved by iteration: starting from an arbitrary initial density,
computing the dust temperature using RADMC-3D (Version 2.0, Dullemond et al.,
2012), determining the resulting gas temperature, solving for vertical
hydrostatic pressure equilibrium, and subsequently updating the density and
temperatures in iterations until convergence. Step 2 computes the
$\mathrm{C^{18}O}$ abundance distribution via the reduced chemical network
that includes isotopologue-selective dissociation and CO to CO2 ice conversion
on grains. It then computes the $\mathrm{C^{18}O}$ line emission using the
radiative transfer tool LIME (Line Modelling Engine Version 1.9.5, Brinch &
Hogerheijde 2010) and compares it to the observed line emission profiles. If
the agreement is poor, then the initial parameters (e.g., surface density
distribution $\Sigma$, gas-to-dust mass ratio $\Delta_{\mathrm{gd}}$) are
modified to repeat the entire modeling procedure from Step 1 until a
satisfactory match with both SED and line emission is obtained. Details about
the two steps are provided in the following subsections.

Figure 1:  Flow chart summarizing our modeling approach. Free parameters are
listed at the top of the chart, and stellar parameters are fixed. These free
parameters include the disk’s inner radius ($r_{\mathrm{in}}$), tapering-off
radius ($r_{\mathrm{tap}}$), power-law index of the surface density
distribution ($p$), dust opacity ($\kappa_{\mathrm{dust}}(a)$), dust disk mass
($M_{\mathrm{dust}}$), and the gas-to-dust mass ratio
($\Delta_{\mathrm{gd}}$). In this model, we assume that the dust and gas are
well-coupled. Therefore, the gas temperature ($T_{\mathrm{gas}}$) represents a
cross-section weighted mean value of the dust grain temperature
($T_{\mathrm{dust}}$), adding the contribution from viscous heating
($T_{\mathrm{visc}}$). Similarly, the gas mass density ($\rho_{\mathrm{gas}}$)
and number density ($n_{\mathrm{gas}}$) are derived from the dust density
($\rho_{\mathrm{dust}}$) by multiplying it with the factors of
$\Delta_{\mathrm{gd}}$ and $\Delta_{\mathrm{gd}}/\mu$, respectively, where
$\mu$ represents the mean molecular mass.

###  2.1 Model Step 1: Finding a Self-consistent Disk Structure that Fits the
SED

The main input parameters for this step (apart from the stellar parameters)
are the surface density distribution, dust size distribution and opacity, and
the disk gas-to-dust ratio.

#### Surface Density

The surface density distribution is assumed to be that of a viscously evolving
disk (e.g., Hartmann et al., 1998) and is specified as:

$\Sigma(r)=\Sigma_{1}(\frac{r}{1\,\mathrm{AU}})^{-p}\exp{\left[-(\frac{r}{r_{\mathrm{tap}}})^{2-\gamma}\right]};(r_{\mathrm{in}}<r<r_{\mathrm{out}})$
(1)

where $r$ is the radial distance from the star, $\Sigma_{1}$ is the surface
density at $1\,\mathrm{AU}$ that is scaled according to the chosen disk mass,
$\gamma$ is the tapering-off exponent and $p$ is the power-law index. We
further assume that $\gamma=p$, which represents the self-similar viscous
solution. The inner radius cut-off, $r_{\mathrm{in}}$, is assumed to be the
dust sublimation radius while the outer radius $r_{\mathrm{out}}$ is chosen to
be much larger than $r_{\mathrm{tap}}$ to ensure that all the mass is
included.

#### Dust Properties

As discussed later, the gas temperature is computed assuming an equilibrium
value from collisional heating/cooling by dust grains. To determine the gas
temperature accurately, we use multiple dust sizes and calculate the size-
dependent dust temperature. The dust species are divided into multiple grain-
size bins equally distributed in log-space. The dust number density follows a
power-law distribution: $n(a)\propto a^{-q}\ \mathrm{with}\
a\in[a_{\mathrm{min}},a_{\mathrm{max}}]$ where $a$ is the dust grain size and
$q$ is the exponent describing the size distribution. We adopt a dust
composition consisting of 64% astronomical silicates and 36% graphite by
volume fraction, which is representative of ISM dust with the ratio of visual
extinction to reddening $R_{V}=5.5$ (Weingartner & Draine, 2001, hereafter
WD01). A similar composition has been adopted in many previous disk models
(e.g., Ansdell et al., 2016; Miotello et al., 2016; Woitke et al., 2019). The
dsharp_opac package from Birnstiel et al. (2018) is used to compute the
wavelength dependence of dust opacity, and the optical constants are those of
astrosilicate from WD01 and graphite from Draine (2003).

#### Gas-to-dust ratio

In order to determine the vertical hydrostatic pressure equilibrium solution,
the gas pressure gradient and hence the gas density are needed. In DiskMINT,
the surface density distribution of gas and dust can in principle be specified
separately as a function of radius, and this determines the local gas-to-dust
ratio $\Delta_{\mathrm{gd}}(r)$. However, for our modeling of RU Lup, we
assumed a constant value throughout the disk for simplicity, which, as we show
later, can already match the $\mathrm{C^{18}O}$ data.

The vertical dust density distribution, $\rho_{{\rm d}}(r,z)$, is initially
set as an arbitrary Gaussian profile. This is then distributed according to
the mass fraction in each grain size bin to obtain $\rho_{{\rm d}}(r,z,a)$.
RADMC-3D (Dullemond et al., 2012) is used to compute the dust temperature
$T_{\mathrm{d}}(r,z,a)$ for each grain size bin. We first determine the gas
temperature $T_{\mathrm{g}}(r,z)$ balancing collisional energy exchange with
dust grains; this contribution is denoted as $T_{\mathrm{g,d}}$. Since RU Lup
is a high accretor, the near and mid-infrared SED can be affected by viscous
heating (Boss & Yorke, 1996). RADMC-3D does not currently include this viscous
heating term, and we hence add this as a separate contribution to the gas
($T_{\mathrm{g,v}}$), and the dust as described later below.

$T_{\mathrm{g,d}}$ is estimated from the following equation balancing dust
heating and cooling

$\begin{array}[]{cc}\displaystyle\sum_{T_{\mathrm{d}}(a)>T_{\mathrm{g,d}}}&A_{\mathrm{H}}n_{\mathrm{d}}(a)\pi
a^{2}n_{\mathrm{H}}\bar{v}_{\mathrm{H}}2k_{B}[T_{\mathrm{d}}(a)-T_{\mathrm{g,d}}]\\\
=\displaystyle\sum_{T_{\mathrm{d}}(a)<T_{\mathrm{g,d}}}&A_{\mathrm{H}}n_{\mathrm{d}}(a)\pi
a^{2}n_{\mathrm{H}}\bar{v}_{\mathrm{H}}2k_{B}[T_{\mathrm{g,d}}-T_{\mathrm{d}}(a)],\end{array}$
(2)

where $T_{\mathrm{d}}(a)$ is the dust temperature at the grain size $a$,
$A_{\mathrm{H}}$ is the mean accommodation coefficient, $n_{\mathrm{d}}(a)$ is
the dust number density distribution, $n_{\mathrm{H}}$ is the gas number
density, $\bar{v}_{\mathrm{H}}$ is the gas thermal velocity, and $k_{B}$ is
the Boltzmann constant. This thermal balance equation simplifies to

$\begin{array}[]{cc}\displaystyle\sum_{T_{\mathrm{d}}(a)>T_{\mathrm{g,d}}}n_{\mathrm{d}}(a)a^{2}\
[T_{\mathrm{d}}(a)-T_{\mathrm{g,d}}]\\\
=\displaystyle\sum_{T_{\mathrm{d}}(a)<T_{\mathrm{g,d}}}n_{\mathrm{d}}(a)a^{2}\
[T_{\mathrm{g,d}}-T_{\mathrm{d}}(a)],\end{array}$ (3)

and only the terms related to dust size remain.

The gas temperature contributed by dust grain collisions ($T_{\mathrm{g,d}}$)
is thus a cross-section weighted mean value between the hot (small) and cold
(large) dust grain temperatures, and therefore the number of grain size bins
($N$) used could potentially affect the accuracy of the gas temperature
evaluation. We adopt $N=20$ as we find that this results in gas temperature
deviations (caused by $N$) to be less than 5%.

We next estimate the temperature due to a balance between accretion heating
and radiative cooling (e.g., Armitage, 2022). Viscous heating is given by
$(9/4)\nu\Sigma\Omega_{k}^{2}$ where $\nu$ is the kinematic viscosity, and
$\Omega_{K}$ is the Keplerian angular frequency. Cooling is given by
$2\sigma_{\mathrm{SB}}T^{4}$, and for a disk accreting in steady state the
accretion rate ($\dot{M}_{\mathrm{acc}}\sim 3\pi\nu\Sigma$) we have

$T_{\mathrm{g,v}}=\left[\frac{3GM_{\star}\dot{M}_{\mathrm{acc}}}{8\pi\sigma_{\rm
SB}r^{3}}\times(1-\sqrt{\frac{r_{\star}}{r}})\right]^{\frac{1}{4}},$ (4)

where $G$ is the gravitational constant and $\sigma_{\mathrm{SB}}$ is the
Stefan–Boltzmann constant. The resulting gas temperature is determined by
adding the two temperatures in quadrature and is therefore given by

$T_{\mathrm{g}}(r,z)=\left[T_{\mathrm{g,d}}^{4}(r,z)+T_{\mathrm{g,v}}^{4}(r,z)\right]^{\frac{1}{4}}.$
(5)

Viscous heating dominates only at the mid-plane in the inner disk ($\lesssim
10\,\mathrm{AU}$) for typical disk densities (also see, e.g., D’Alessio et
al., 1998).

Once the gas temperature is computed, the new density structure is calculated
from the pressure gradient by solving

$\frac{dP(r,z)}{dz}=-\rho_{\mathrm{gas}}(r,z)\Omega^{2}z,$ (6)

where $P(r,z)$, $\rho_{\mathrm{gas}}$ and $\Omega$ are the gas pressure, gas
density (assumed to be the total dust density times a constant
$\Delta_{\mathrm{gd}}$) and Keplerian frequency, respectively. For the next
iteration, the dust density profile with $z$ is rescaled with this vertical
gas density profile, and re-normalized to the surface density at this radius.
The dust temperatures are re-calculated with the new dust density distribution
using RADMC-3D. The steps above are recomputed until convergence is achieved
at the iteration $m$:
$|\left[\rho_{\mathrm{gas},m}(r,z)-\rho_{\mathrm{gas},m-1}(r,z)\right]/\rho_{\mathrm{gas},m-1}(r,z)|<5\%$
for regions with $\rho_{\mathrm{gas}}>10^{-20}\,\mathrm{g\,cm^{-3}}$
(corresponding to $n_{\rm H}\gtrsim 10^{3}\,\mathrm{cm^{-3}}$). The error
tolerance was chosen as a reasonable compromise between accuracy and speed of
computation ($\sim 4$ hours to achieve convergence when running with 24
threads with 2.10 GHz CPUs). Lower tolerances did not significantly change the
results.

The above procedure results in a dust and gas density and temperature
distribution which are all self-consistent with the local vertical pressure
gradient. We described viscous heating for gas above, but this term is also
relevant for heating dust grains. Since this is difficult to incorporate into
the RADMC-3D code, we include this effect by adding it to the dust grains
before computing the SED. This is done by considering the gas as a thermal
reservoir that equilibrates the dust temperature in regions where dust and gas
are highly coupled. In practice, we estimate the extent of this mid-plane
region as the region where the temperature differences between the
hottest/smallest grain and coldest/largest grain are small enough as
$|\left(T_{\mathrm{d}}(a_{\mathrm{min}})-T_{\mathrm{d}}(a_{\mathrm{max}})\right)/T_{\mathrm{d}}(a_{\mathrm{max}})|<10\%$.
$T_{\mathrm{d}}(a)=T_{\mathrm{g}}$ is set for all grain sizes in this coupled
region. We then run RADMC-3D to compute the SED and compare it with the
observed SED.

We vary the disk dust parameters until a satisfactory match to the SED is
obtained. The dust opacity $\kappa_{\nu}$ and the dust mass
$M_{\mathrm{dust}}$ are two main parameters affecting the synthetic SED:
Changing $\kappa_{\nu}$ alters the slope of the long-wavelength portion of the
SED and $M_{\mathrm{dust}}$ moves the flux density up and down. In practice,
we find the best fit $\kappa_{\nu}$ by comparing the slope of the dust opacity
$\beta_{\mathrm{abs}}=-\frac{d\log(\kappa_{\mathrm{abs}})}{d\log{\lambda}}$
with the slope of the SED at long wavelength
$\alpha_{\mathrm{SED}}=-\frac{d\log{F_{\nu}}}{d\log{\lambda}}$ ($\lambda\gg
100\,\micron$) based on the relation between the two slopes
$\beta_{\mathrm{abs}}=\alpha_{\mathrm{SED}}-2$. When the dust composition is
fixed, we first vary the maximum particle size $a_{\mathrm{max}}$ and keep the
slope of the number density distribution with size fixed to $q=3.5$, which is
the value expected in collisional equilibrium (Birnstiel et al., 2011). If the
upper limit of $a_{\mathrm{max}}=1\,\mathrm{cm}$ is reached while varying
$a_{\mathrm{max}}$, then $q$ is varied to find the best match of the slope.
After the best-fit $\kappa_{\nu}$ is found, the $M_{\mathrm{dust}}$ is derived
by matching the absolute value of the flux density at long wavelengths.

###  2.2 Model Step 2: Computing the $\mathrm{C^{18}O}$ Line Emission and
Profile

The next step in our modeling approach is to run the reduced chemical network
described in RGH22 to obtain the $\mathrm{C^{18}O}$ abundance with $(r,z)$.
The photodissociation rates (for our application target RU Lup) are computed
from the UV _HST_ /COS median-resolution spectrum obtained by France et al.
(2014) (see also Figure 2 for average photometric values from this spectrum).
We assume all gas is molecular in the disk structure calculation but
explicitly solve for the chemistry by specifying the corresponding H nuclei
density ($n_{\mathrm{H}}$) for the chemical network. This means that all
molecular abundances in the chemical network are defined by their density
ratio compared to the density of H nuclei. Finally, the gaseous abundances of
$\mathrm{C^{18}O}$ and the disk structure are inputs to LIME (Brinch &
Hogerheijde, 2010) to compute the non-LTE(local thermal equilibrium) synthetic
$\mathrm{C^{18}O}$ (2-1) and (3-2) emission.

The model parameters are varied until the synthesized SED and
$\mathrm{C^{18}O}$ line emission match the observations. We fix $\kappa_{\nu}$
and $M_{\mathrm{dust}}$ to the values determined in the SED fitting, and
explore a range of gas-to-dust ratios $\Delta_{\mathrm{gd}}=5,10,50,100$
(which covers the low disk $\Delta_{\mathrm{gd}}$ reported in the literature
up to the ISM value) to generate a grid of $M_{\mathrm{gas}}$. Since the self-
consistent VHSE solution depends on the gas mass (which varies with
$\Delta_{\mathrm{gd}}$ in the gas mass grid), the vertical density structure
of each of these models slightly differ. However, the SED at long wavelengths
traces the optically thin thermal emission from the large grains and remains
the same as it is not sensitive to the vertical dust density distribution. The
derived $M_{\mathrm{dust}}$ therefore remains unaltered even as
$\Delta_{\mathrm{gd}}$ is varied. We start from the beginning for each grid
point and find that we do not need to re-fit the SED, hence we calculate the
dust thermal structure with RADMC-3D, solve the VHSE through iterations, and
then derive the $\mathrm{C^{18}O}$ abundance by the reduced chemical network.
Next, the $\mathrm{C^{18}O}$ line luminosity ($L_{\mathrm{C^{18}O}}$) is
computed to compile a $L_{\mathrm{C^{18}O}}$ vs. $M_{\mathrm{gas}}$ relation.
The best-fit $M_{\mathrm{gas}}$ is then determined as the value where the
modeling relation ($L_{\mathrm{C^{18}O}}$ vs. $M_{\mathrm{gas}}$) intersects
the luminosity inferred from the observations. Finally, we run the model with
best-fit $M_{\mathrm{gas}}$ again also from the beginning to verify the
estimate found above.

In this work, we not only compare total line luminosities as in RGH22 but also
match the velocity profile and radial distribution of the $\mathrm{C^{18}O}$
(2-1) line. These are generated by the Python package GoFish (Teague, 2019)
from the simulated LIME image and follow the same procedure used on
observational data. The slope of the surface density distribution and the gas-
to-dust ratio as a function of radius are parameters that can be changed to
improve the fit on the line profile, if necessary.

##  3 Application to RU Lup

###  3.1 The Highly Accreting RU Lup Star and Its Dust and Gas disk

RU Lup (Sz 83, 2MASS J15564230-3749154) is a K7-type star located at a
distance of $158.9\,\mathrm{pc}$ (Gaia Collaboration et al., 2018) and a
member of the Lupus II star-forming region (Comerón, 2008). RU Lup has the
highest mass accretion rate ($\sim 10^{-7}\,M_{\odot}$/yr, Alcalá et al. 2017)
and is one of the most active stars in the region with large irregular
variations in both spectroscopy and photometry from ultraviolet (UV) to
infrared (IR) wavelengths (e.g., Hughes et al., 1994; Herczeg et al., 2005;
Gahm et al., 2013). The stellar mass estimates range from 0.2 to 1.2
$M_{\odot}$ (e.g., Alcalá et al., 2017; Andrews et al., 2018; Yen et al.,
2018). Here, we adopt the value of $0.7\,M_{\odot}$ from more recent
evolutionary models (Alcalá et al., 2017) over the dynamical mass of 0.2
$M_{\odot}$. This is because the disk of RU Lup is close to face-on which
introduces a large uncertainty in the dynamical mass (Yen et al., 2018). As
one of the most extensively observed Class II objects in Lupus, photometry and
spectra are available from the UV to radio wavelengths resulting in the multi-
wavelength spectral energy distribution (SED) shown in Figure 2, where average
photometry is reported for multi-epoch observations.

Figure 2:  Observational data that used in this work. Top panel: spectral
energy distribution (SED). Photometry in colored markers and Spitzer/IRS
spectrum in magenta. A combined BT-settl stellar photosphere spectrum (Allard
et al., 2003, 2011) with a 10,500 K blackbody radiation fitting the UV data-
points is shown as a black line. Lower panels: ALMA $\mathrm{C^{18}O}$ (2-1)
velocity profile (lower left) and radial profile (lower right), data from
Huang et al. (2020) with uncertainties calculated by GoFish. Linear
interpolation (grey line) is adopted to recover the channel contaminated by a
dark cloud (dashed black line) in the line spectrum. SED references:
photometry from the _HST_ spectrum in France et al. (2014); optical photometry
from Gras-Velázquez & Ray (2005), APASS (Henden et al., 2016), _Gaia_ (Gaia
Collaboration et al., 2018); infrared photometry from 2MASS (Skrutskie et al.,
2006), _WISE_ (Wright et al., 2010), _AKARI_ (Ishihara et al., 2010); and
(sub)millimeter data from Ansdell et al. (2016), DSHARP (Andrews et al.,
2018), and Lommen et al. (2007, 2009).

A large-scale, complex proto-planetary disk has also been recently revealed by
ALMA. The millimeter dust disk appears symmetric with multiple annular gaps
and rings and extends out to a radius of $\sim 63\,\mathrm{AU}$ (Huang et al.,
2018, 2020). In contrast, CO emission has a more asymmetric morphology. Huang
et al. (2020) identified a Keplerian disk with a radius of $\sim
120\,\mathrm{AU}$, similar in size to that inferred via scattered light
(Avenhaus et al., 2018), surrounded by an envelope extending out to $\sim
260\,\mathrm{AU}$ with spiral arms and clumps.

However, the $\mathrm{C^{18}O}$ emission, which we focus on and aim to model
in this work, is less complex. The $\mathrm{C^{18}O}$ emission is symmetric,
only traces the Keplerian disk, and has a radius of $\lesssim
100\,\mathrm{AU}$. The lower panels of Figure 2 show the $\mathrm{C^{18}O}$
(2-1) line profile and radial intensity cut from publicly available datacubes
(Huang et al., 2020) generated using GoFish. We choose the same aperture and
wavelength range used in Huang et al. (2020), $1.5\arcsec\sim
240\,\mathrm{AU}$ and $1.75-7.25\,\mathrm{km/s}$, as the maximum extent to
include all the emitting areas and channels when computing the line profile
and radial profile. In the line profile, there is clear dark cloud
contamination at $5\,\mathrm{km/s}$ (dashed line). Linear interpolation (grey
point and line) is utilized to recover the disk emission in this channel,
which brings the integrated total flux of $\mathrm{C^{18}O}$ (2-1) from
$0.34\pm 0.03\,\mathrm{Jy\,km/s}$ to $0.37\pm 0.03\,\mathrm{Jy\,km/s}$. For
the radial profile, deprojection is applied using the disk position angle
$\mathrm{PA}=121\arcdeg$ and inclination $i=18.8\arcdeg$ (Huang et al., 2018).

The dust and gas mass of RU Lup, hence the gas-to-dust mass ratio
$\Delta_{\mathrm{gd}}$, have been previously estimated using continuum
millimeter emission and CO isotopologue emission. Ansdell et al. (2016)
measured the $\mathrm{{}^{13}CO}$(3-2) and $\mathrm{C^{18}O}$(3-2) line fluxes
and compared them to a grid of simple disk models by Williams & Best (2014):
They inferred a gas disk mass of $\sim 2.7_{-1.7}^{+7.3}\times
10^{-3}\,M_{\odot}$ and a gas-to-dust mass ratio $\Delta_{\mathrm{gd}}\sim
8.9_{-5.7}^{+24.4}$. Miotello et al. (2014, 2016, 2017) included isotope-
selective dissociation in the thermo-chemical physical code DALI (Bruderer et
al., 2012; Bruderer, 2013) and used the same line luminosities to infer an
even lower gas disk mass ($\sim 1.5_{-1.0}^{+2.5}\times 10^{-3}\,M_{\odot}$)
and gas-to-dust ratio ($\Delta_{\mathrm{gd}}\sim 3.8_{-2.6}^{+6.2}$). Clearly,
the low inferred disk mass and gas-to-dust mass ratio are hard to reconcile
with the large dust and gas disk of RU Lup and the high accretion rate onto
the star; we therefore re-examine the dust and gas mass constraints using the
DiskMINT modeling approach.

###  3.2 Specific Models

Two models are considered in this work with different vertical density
distributions: (a) the VHSE model uses a self-consistent vertical hydrostatic
pressure equilibrium solution; (b) the Gaussian model uses a parameterized
Gaussian vertical structure. Both models share the same surface density
distribution and use the same dust grains (same $\kappa_{\nu}$ and
$M_{\mathrm{dust}}$) determined by fitting the long wavelength portion of the
SED ($\lambda\gtrsim 100\,\micron$). The Gaussian model additionally fits the
IR wavelengths ($10\,\micron\lesssim\lambda\lesssim 100\,\micron$) by assuming
the pressure gradient to be a free parameter and thus varying pressure scale
height in the Gaussian structure:
$H_{p}=H_{p,\mathrm{100}}(r/100\,\mathrm{AU})^{\alpha}$ with free
characteristic height $H_{p,100}$ and flaring index $\alpha$. This is the
approach taken in a few recent studies to estimate disk masses and
$\Delta_{\mathrm{gd}}$ (e.g., Woitke et al., 2019; Zhang et al., 2021).

The model input parameters are presented in Table 1. Stellar mass $M_{\star}$,
radius $r_{\star}$ as well as mass accretion rate $\dot{M}_{\mathrm{acc}}$ are
fixed and taken from the literature (see Section 3.1). The inner radius
$r_{\mathrm{in}}$ is fixed at the dust sublimation radius, and the tapering-
off radius is set as the dust outer radius $r_{\mathrm{tap}}\sim
r_{\mathrm{dust}}$ given in Huang et al. (2018). The dust opacity
$\kappa_{\nu}$ is computed by dsharp_opac: It uses a dust composition
described in Section 2.1, and has fixed volume fraction from $a_{\rm min}\sim
1.0\times 10^{-6}\,\mathrm{cm}$ through to $a_{\rm max}$, in which $a_{\rm
max}$ and the power law index $q$ are free parameters. The other two free
parameters are the dust disk mass $M_{\mathrm{dust}}$, and the gas-to-dust
mass ratio $\Delta_{\mathrm{gd}}$.

The synthetic imaging setup for the models is obtained from observations
(summarized in Section 3.1). The output synthetic image is created with a
pixel size of 0.04″, and with source distance, $i$ and $\mathrm{PA}$. The
image has $151\times 151$ pixels to include all disk emission within 3.0″. The
dust continuum emission is also included in the synthetic image, and then the
continuum is subtracted in the final line imaging datacube. Then, the LIME
output image is convolved with a beam of $0.32\times 0.32\,\arcsec$ to get the
final synthesized image.

Table 1: Model Parameters

Parameter | Symbol | Value
---|---|---
Dust Properties |  |  |  |
Volume fraction |  |  | 64% Silicate | 36% Graphite
minimum size | $a_{\mathrm{min}}$ |  $1\times 10^{-6}$ cm
maximum size | $a_{\mathrm{max}}$ | free parameter
exponential slope | $q$ | free parameter
Radial Structure |  |  |  |
inner radius of the disk | $r_{\mathrm{in}}$ | 0.035 AU
tapering-off radius | $r_{\mathrm{tap}}$ | 63 AU
surface density slope | $p$ | 1
Vertical Structure |  |  | VHSE | Gaussian
Characteristic Scale Height | $H_{p,\mathrm{100}}$ | solved | free parameter
Flaring Index | $\alpha$ | solved | free parameter

Note. — In principle, all parameters in this table could be varied to fit the
observations. However, only $a_{\mathrm{max}}$ and $q$ are changed here for
the VHSE model as the default settings for other parameters could already give
a good fit. $H_{p,\mathrm{100}}$ and $\alpha$ are also set free for the
Gaussian model while the vertical structure for the VHSE model is solved self-
consistently from pressure equilibrium. The best-fit free parameters are
summarized in Table 2.

##  4 Results and Discussion

The inferred dust parameters, dust and gas masses are summarized in Table 2.
One of the main results of this work is that our model can explain RU Lup’s
long-wavelength ($\lambda\gtrsim 100\micron$) SED, the $\mathrm{C^{18}O}$
(2-1) and (3-2) line luminosities, and the velocity and radial profiles, with
a higher $M_{\mathrm{gas}}$ and thus higher $\Delta_{\mathrm{gd}}$ than
previously inferred. We present details on these models in Section 4.1.
Effects of $\mathrm{CO}\leftrightarrow\mathrm{CO_{2}}$ conversion on grain-
surface and differences between the Gaussian and VHSE models are discussed in
Section 4.2. Our VHSE model under-estimates the strong IR excess of RU Lup by
a factor of $\lesssim 3$, and we discuss possible reconciliations in Section
4.3.

Table 2: Model Main Results  Model | $H_{p,\mathrm{100}}$ | $\alpha$ | $a_{\mathrm{max}}$ | $q$ | $M_{\mathrm{dust}}$ | $M_{\mathrm{gas}}$ | $\Delta_{\mathrm{gd}}$
---|---|---|---|---|---|---|---
| (AU) |  | (cm) |  | $(M_{\odot})$ | $(M_{\odot})$ |
VHSE | - | - | 0.3 | 3.5 | $4.0\times 10^{-4}$ | $1.2\times 10^{-2}$ | 30
Gaussian | 30 | 1.1 | 0.3 | 3.5 | $4.0\times 10^{-4}$ | $2.1\times 10^{-2}$ | 52

###  4.1 $\mathrm{C^{18}O}$ Emission Indicates a Relatively High Gas Disk Mass
for RU Lup

In DiskMINT, the disk density structure is based on the dust temperature
profile, and the dust disk is constructed by fitting the SED (See Section
2.1). The SED fits of the two models (VHSE and Gaussian) introduced in Section
3.2 are shown in the top panel of Figure 3. Both models share the same dust
grain properties described in Table 1 and the best-fit free parameters are
reported in Table 2. The best-fit maximum grain size $a_{\mathrm{max}}$, $q$
parameters and dust disk mass $M_{\mathrm{dust}}$ are the same for both
models: $a_{\mathrm{max}}=0.3\,\mathrm{cm}$, $q=3.5$ and
$M_{\mathrm{dust}}\sim 4.0\times 10^{-4}\,M_{\odot}$. Since the pressure scale
height is determined using free parameters to match the SED in the Gaussian
model, this model provides a better fit to the IR SED. To find the best
parameters, we start from the best-fit pressure scale height with
$H_{p,100}=20.9324\,\mathrm{AU}$ and $\alpha=1.1301$ reported in Woitke et al.
(2019) for the disk of RU Lup, and generate a grid of
$H_{p,100}=15,20,25,30\,\mathrm{AU}$ and $\alpha=1.05,1.10,1.15,1.20$.
Although we use a different dust composition and updated parameters for the
central star, we find a relatively close pressure scale height with
$H_{p,100}=30\,\mathrm{AU}$ and $\alpha=1.10$ and a very similar synthetic SED
for the Gaussian model. For the VHSE models, the procedure of iteration to
determine the vertical density structure to be consistent with the temperature
profile sets the pressure scale height; there is no simple power law to
describe the scale height thus the parameters $H_{p,100}$ and $\alpha$ are not
valid.

Figure 3:  Top panel: observed photometry compared to modeled SEDs. The
uncertainties of the observation are shown in the gray shade. The VHSE and
Gaussian models are shown with brown and blue lines, respectively. Lower
panels:$\mathrm{C^{18}O}$ (2-1) velocity profile (left) and radial cut (right)
in black compared to modeled profiles. Color-coding for the models as in the
top panel.

We then run the reduced chemical network and LIME to obtain the synthetic
$\mathrm{C^{18}O}$ luminosity, which is compared with the observation to
obtain the best-fit $\Delta_{\mathrm{gd}}$ and hence the $M_{\mathrm{gas}}$
(see Section 2.2). The synthetic $\mathrm{C^{18}O}$ (2-1) and (3-2)
luminosities vs. $M_{\mathrm{gas}}$ using different $\Delta_{\mathrm{gd}}$ are
presented in Figure 4. There are four data points on each modeling line
representing $\Delta_{\mathrm{gd}}=5,10,50,100$ (points in Figure 4), and one
additional best-fit model (‘$\times$’ in Figure 4), which is obtained at the
cross point with the observation of $\mathrm{C^{18}O}$ (2-1) luminosity at
left panels. The best-fit gas masses for both models are within a factor of
two: $M_{\mathrm{gas}}\sim 1.2\times 10^{-2}\,M_{\odot}$ and $\sim 2.1\times
10^{-2}\,M_{\odot}$ for the VHSE and Gaussian models, respectively.

In addition to matching the luminosity, the line spectrum and radial
distribution from the synthetic $\mathrm{C^{18}O}$ line images are also
compared with the observations. The lower panels of Figure 3 present the
$\mathrm{C^{18}O}$ (2-1) spectra and radial distribution for different models
generated from simulated LIME datacubes with GoFish and the same setup (see
Section 3.2) as for the observational data. These panels demonstrate that the
VHSE model with default input parameters (Table 1) also matches the
$\mathrm{C^{18}O}$ (2-1) line profile and radial cut. The Gaussian model can
reproduce the line luminosity and matches the $\mathrm{C^{18}O}$ (2-1) line
velocity profile relatively well, but its emission is more compact than the
VHSE model with the intensity peaking closer to the host star. We note that
the models also fit the $\mathrm{C^{18}O}$ (3-2) luminosity (Ansdell et al.,
2016), as shown in the right panel of Figure 4. The $\mathrm{C^{18}O}$ (3-2)
line emission has a similar velocity profile and radial cut, but it is a
factor of $\sim 5$ more luminous than the (2-1) line. For both models, even
better fits may be achieved by changing the surface density distribution and
by including a radial-dependent gas-to-dust ratio, but we did not consider
these modifications necessary for RU Lup.

The Gaussian model has an emission profile that is less radially extended
compared with the observations. This is because it has a very puffed-up
density distribution which appears necessary to fit the IR SED:
$H_{p,100}=30\,\mathrm{AU}$ by this work (also $H_{p,100}\sim 21\,\mathrm{AU}$
reported in Woitke et al. (2019) which is a better match to the SED). Since
the scale height is parameterized as a power-law, this implies that the
flaring index in the outer $\mathrm{C^{18}O}$ emitting regions of the model
disk is also higher. The increased flaring moves the $\mathrm{C^{18}O}$
emitting layer closer to the star and higher. It is nearly a factor of $\sim
3-4$ higher than the VHSE disk at $r\sim 50\,\mathrm{AU}$ where most of the
$\mathrm{C^{18}O}$ emission comes from (Figure 5). Although it is hard to
obtain the height of the $\mathrm{C^{18}O}$ emitting layer for the RU Lup disk
due to its small inclination angle, this unrealistically puffed-up disk
scenario – with the emitting layer as high as $z/r\sim 0.8$ at $r\sim
50\,\mathrm{AU}$ – is at odds with recent observations which instead find the
$\mathrm{C^{18}O}$ emitting layer of Class II disks to be at $z/r\sim 0.1$ for
$r<100\,\mathrm{AU}$ (Paneque-Carreño et al., 2023). We note that, in
principle, if the height of the emitting layer could be measured as it has
been in some disks, then this information could be used to fit the
$\mathrm{C^{18}O}$ radial and velocity profiles for the Gaussian model. We
also find that if we assume the scale height obtained from the VHSE model and
repeat the Gaussian modeling for RU Lup, it results in a combination of
($M_{\mathrm{gas}}$, $L_{\mathrm{C^{18}O}}$) similar to the best-fit VHSE
model, although the synthetic IR SED is no longer an improved match to the
data.

While it may be possible to fit all of the observational data using a Gaussian
disk model, determining the emission scale height requires very high spatial
resolution observations and only works for disks with favorable inclination
angles. In their absence, the disk structure parameterization can deviate
substantially from reality as we show for RU Lup. On the other hand, the VHSE
model is physically motivated, determines the scale height at each radius via
coupling of the disk density and temperature structure, and can simultaneously
fit the radial and velocity distribution of flux. Hence, we believe it is a
more reliable indicator of conditions in the disk.

In summary, our VHSE model fits the SED, total $\mathrm{C^{18}O}$ line
emission, velocity, and radial profiles from recent observations (e.g.,
Ansdell et al., 2016; Huang et al., 2020) with relatively high
$M_{\mathrm{gas}}$ and $\Delta_{\mathrm{gd}}$ ($\sim 30$) in comparison with
the previously inferred $\Delta_{\mathrm{gd}}$ of $\sim 4$. Using a Gaussian
vertical distribution, our model also derives a similarly high
$M_{\mathrm{gas}}$ within a factor of $\sim 2$ of the one obtained from the
VHSE model. Thus, we conclude that the RU Lup disk is not significantly low in
its gas mass and nor has it undergone any substantial change in CO chemistry
due to changes in C/H and O/H caused by planet formation processes. We confirm
the conclusions of RGH22, and find that optically thin $\mathrm{C^{18}O}$
lines provide reasonable estimates of the disk mass. We also note that the
RGH22 models compare favorably not only with the $\mathrm{C^{18}O}$ fluxes,
but also with the $\mathrm{{}^{13}CO}$, $\mathrm{CO}$, and atomic carbon
forbidden line $\mathrm{[CI]}$ fluxes for a sample of large disks (R$\geq$200
AU) (Pascucci et al., 2023), and cold water emission as well (Ruaud & Gorti,
submitted).

Figure 4:  Synthetic $\mathrm{C^{18}O}$ (2-1) and (3-2) luminosity vs.
$M_{\mathrm{gas}}$. In both panels, the black horizontal line marks the
luminosity inferred from the data and grey shades show their uncertainties:
(2-1) line from Huang et al. (2020) and (3-2) from Ansdell et al. (2016). The
VHSE models introduced in this work are shown in points and line in brown. The
Gaussian models are presented in squares and line in blue. In each modeling
scenario in this work, four models with $\Delta_{\mathrm{gd}}=5,10,50,100$
were simulated, and one additional best-fit model was simulated at the cross-
match ‘$\times$’ point with observations.  Figure 5: Gas density distribution,
gas temperature, $\mathrm{C^{18}O}$ density distribution, $\mathrm{C^{18}O}$
(2-1) and (3-2) line emitting layers (normalized to the total luminosity) are
presented from left to right, VHSE model top and Gaussian model bottom row. In
the first two panels, dashed black lines mark the boundaries of the C18O
location where $n(\mathrm{C^{18}O})\geq 0.1$ and emitting region, and upper &
lower white lines represent the $A_{\rm V}\sim 1\,\&\,10$, respectively. While
the $\mathrm{C^{18}O}$ (3-2) flux is higher than the (2-1), the emitting
layers are very similar.  Figure 6:  Synthetic $\mathrm{C^{18}O}$ (2-1) and
(3-2) luminosity vs. $M_{\mathrm{gas}}$. The observation results are presented
in black. Two best-fit models are simulated at the cross-match point with
observations represented by ‘$\times$’ in brown and blue following the color
in Figure 4. All possible results from full VHSE models in Ruaud et al. (2022)
(RGH22) are shown in magenta shades. The Gaussian models without grain-surface
chemistry from Miotello et al. (2017) are shown in green shades. The median
value of all possible estimated luminosity for each $M_{\mathrm{gas}}$ is
presented in dashed magenta and green lines, respectively. Our model predicts
a relatively large gas disk mass ($\sim 1\times 10^{-2}\,M_{\odot}$) lines
between the model grid of RGH22 and Miotello et al. (2017).

###  4.2 Comparisons with Literature Values

Our work is the first to focus on specifically modeling CO isotopologue
emission from RU Lup, and matches the SED, $\mathrm{C^{18}O}$ line luminosity,
spectrum and radial profile. In this section, we compare our source-specific
model with the grids generated in previous works and discuss possible
explanations for the different results in gas masses and gas-to-dust ratios.

First, we comment on the differences between the VHSE model presented here and
those in RGH22. Here our dust-temperature based model gives a factor of $\sim
2$ larger $L_{\mathrm{C^{18}O}}$ compared with the full VHSE thermo-chemical
model by RGH22 (all possible results are shown in Figure 6 magenta regions).
This is approximately consistent with the differences found by RGH22 for the
dust and gas temperature based modeling, and a similar result of a factor of
$\sim 2$ difference was found at these $M_{\mathrm{gas}}$ values.

We also note a few additional differences. We use similar grain-surface
chemistry (as the reduced network was adopted from tests conducted in RGH22),
but use the dust temperature to set our gas temperature whereas RGH22 computed
the gas temperature. We also do not include settling, while in RGH22 most of
the $a\gtrsim 100\,\micron$ settles and plays a negligible role in the thermal
balance, because the balance is dominated by the small grains that has higher
density (Equation 3). Another important difference is the dust composition
used in our models vs. RGH22. In this work, we adopt a combination of
astrosilicate and graphite based on WD01 – that is similar to the dust
composition used in Miotello et al. (2016, 2017) – while RGH22 used a mix of
olivine (76% by volume) and amorphous carbon (24% by volume); more
importantly, we construct the dust disk by fitting the SED of RU Lup. How
different dust compositions affect the disk structure, temperature and grain-
surface chemistry, and how they could be better constrained are out of the
scope of this paper and will be the subject of future work.

We find similar differences in the models for RU Lup from Miotello et al.
(2017), although the dust composition used in our model is similar to theirs.
Their $M_{\mathrm{dust}}$ estimation was derived from the flux at mm-
wavelength and not by fitting the SED, but the dust mass estimation of the two
models converge to the same $M_{\mathrm{dust}}\sim 4\times
10^{-4}\,M_{\odot}$. However, their best-fit value of $M_{\mathrm{gas}}$ is a
factor of $\sim 8$ smaller than the VHSE result and $\sim 14$ smaller than our
Gaussian gas disk model estimate. This can be partially attributed to the fact
that the grid of Gaussian disk models used in Miotello et al. (2017) are not
tailored to RU Lup. For example, as noted earlier, the scale height parameters
adopted impact the inferred line luminosity and therefore the mass estimate.
For the range of scale height parameters (together with other free parameters)
used in Miotello et al. (2017), the mass estimates in fact range from $4\times
10^{-3}\,M_{\odot}$ to $4.8\times 10^{-4}\,M_{\odot}$. Another contributor is
the $\mathrm{CO}\leftrightarrow\mathrm{CO_{2}}$ grain-surface chemistry
conversion which is not accounted for in Miotello et al. (2017); this could
bring a discrepancy of a factor $\sim 2-3$ as noted by Trapman et al. (2021)
and RGH22.

We would like to note that there could be other processes at work that may
deplete gas-phase CO at the surface, e.g., vertical diffusion of gas into the
icy midplane where it may freeze out, although the extent to which this occurs
will also depend on the ability of small grains to form and transport ices
back into the surface layers (e.g., Krijt et al., 2020; Powell et al., 2022).
However, to correctly consider those processes require a full 2D transport
model including the particle dynamics. Such simulations are not suitable for
detailed modeling of observational data on individual targets, as they require
knowledge of the disk’s history; in fact, modeling presented here may help
decipher disk conditions at different evolutionary stages from observations
and inform the development of theoretical transport models.

In summary, the derived $M_{\mathrm{gas}}$ for RU Lup in this work lies
between the model grids from RGH22 and Miotello et al. (2017), see Figure 6.
Our model is the first one that is specifically built for RU Lup. We also fit
the SED, $\mathrm{C^{18}O}$ line spectrum, and radial distribution, while both
previous models only matched the luminosity using a grid of models which
resulted in larger uncertainties on the derived parameters. We thus
demonstrate that DiskMINT is a promising tool for modeling individual disks
and deriving more robust disk mass estimates.

###  4.3 The Missing IR Emission in VHSE Models

As discussed so far, the VHSE model successfully reproduces the
$\mathrm{C^{18}O}$ observations including the line velocity profile and radial
distribution. While the VHSE model presented in this work is capable of
fitting the entire SED of the average $\sim 1-3$ Myr-old disk (e.g., the
median Taurus SED from Furlan et al., 2006), and can also match all available
continuum photometry of RU Lup beyond $\gtrsim 100\,\micron$, it
underestimates the infrared emission from the disk of RU Lup by a factor of
$\sim 2$ between $\sim 2-60\,\micron$ (Figure 7 upper panel).

We first check and confirm that this IR continuum underestimation does not
affect the gas mass determination from the $\mathrm{C^{18}O}$ (2-1) line.
RADMC-3D simulations show that the IR continuum emission comes from within a
radial distance of $10\,\mathrm{AU}$ (see the cumulative dust emission in
Figure 7 lower panel), but the $\mathrm{C^{18}O}$ line emission mostly arises
from the outer disk radius (Figure 3 lower right panel). There is therefore a
deficit of dust emission from within $\lesssim 10$ AU, indicating a possible
missing physical process in our simple disk models.

This lack of strong IR emission in VHSE models has also been noted previously,
e.g., Woitke et al. (2016) for T-Tauri stars and Davies et al. (2018) for
Herbig Ae/Be stars. Moreover, RU Lup has one of the strongest IR excesses, a
factor of $\sim 2$ higher than the upper boundary of the Taurus median SED
(Figure 7), a region of similar age to Lupus (Comerón, 2008; Kenyon et al.,
2008).

Figure 7:  Top panel: SED fitting results for the VHSE model. The same
notations used in Figure 3 are adopted. The median SED from Taurus (Furlan et
al., 2006) is added here with a blue shade. Lower panel: normalized cumulative
dust continuum fluxes for the best-fit VHSE models as a function of the
distance from the star at different wavelengths. Note that even at 90 µm 80%
of the flux arises within $10\,\mathrm{AU}$.

One obvious shortcoming of our VHSE models is that we ignore gas thermal
processes that are important, especially at the surface of the disk at small
radii. Here other heating processes – notable stellar high energy X-ray and UV
photons – will heat the gas to higher temperatures. When densities are high,
gas and dust are better coupled which leads to more small dust at higher
elevations, increasing the IR excess.

Another intriguing possibility is that small dust grains ($\sim\micron$ size)
are uplifted by a wind in the inner part of the disk (Pascucci et al. 2022 for
a recent review on disk winds). This would lead to hotter dust at a higher
scale height and thus increase the IR emission (Bans & Königl, 2012). A
parametric wind and disk model has previously been used to fit the strong IR
excess from an Herbig disk (Fernandes et al., 2018). Interestingly, RU Lup has
a well-known inner wind detected via optical forbidden lines (e.g., Fang et
al., 2018; Banzatti et al., 2019; Whelan et al., 2021). It is quite likely
that the wind (if dense, i.e., $n\gtrsim 10^{6-7}$ cm-3) can loft small
amounts of dust to greater heights and can explain the factor of $\sim 2$
deficit in the IR excess we find with the VHSE models. The hypothesis of a
wind lifting small dust and increasing the IR excess warrants further
exploration.

##  5 Summary and Outlook

We developed a dust temperature-based self-consistent vertical hydrostatic
pressure equilibrium disk model, DiskMINT, to compute gas disk masses.
DiskMINT is a Python code built on RADMC-3D and LIME for the continuum and gas
line radiative transfer, respectively; and it includes a reduced chemical
network suggested in RGH22 to determine the $\mathrm{C^{18}O}$ distribution.
With DiskMINT, we introduce a target-based approach to estimate the disk mass
in considerable details, where we fit the SED and also the $\mathrm{C^{18}O}$
line emission.

We further test it on RU Lup, whose disk was previously inferred to have just
over a Jupiter-mass gas disk ($M_{\mathrm{gas}}\sim 1.5\times
10^{-3}\,M_{\odot}$) and a gas-to-dust mass ratio $\Delta_{\mathrm{gd}}$ of
only $\sim 4$. We show that our model can match the long wavelength portion of
the SED, the total $\mathrm{C^{18}O}$ (2-1) and (3-2) line luminosity as well
as the $\mathrm{C^{18}O}$ (2-1) velocity and radial profiles with an order of
magnitude higher mass ($M_{\mathrm{gas}}\sim 1.2\times 10^{-2}\,M_{\odot}$)
and gas-to-dust ratio ($\Delta_{\mathrm{gd}}\sim 30$). We also test a Gaussian
vertical density distribution that fits the SED better from IR- to millimeter-
wavelengths and considers $\mathrm{CO}\leftrightarrow\mathrm{CO_{2}}$
conversion. We find this Gaussian model that can match the line luminosity
with even higher gas mass ($M_{\mathrm{gas}}\sim 2.1\times
10^{-2}\,M_{\odot}$) and gas-to-dust ratio ($\Delta_{\mathrm{gd}}\sim 52$) but
consider its large vertical height unrealistic. We also find that the VHSE
model underestimates the IR SED ($\lambda\sim 2-60\,\micron$) by a factor of
$\sim 2$, which may indicate the need for considering more detailed gas
thermal balance in the inner disk and/or an inner dusty wind from RU Lup.

With our target-based approach, the RU Lup’s estimated disk mass is better-
constrained, and it is larger than the Minimum Mass Solar Nebula (Hayashi,
1981). The larger mass is more in agreement with the young-age, high accretion
rate, large disk size, and a lack of strong radial substructures in the disk
of RU Lup. Our derived $\Delta_{\mathrm{gd}}$ is just a factor of a few lower
than the ISM value of $\sim 100$. This may indicate the disk has lost some of
its gas within $\lesssim 1\,\mathrm{Myr}$, or alternately, that CO is depleted
by a factor of few. If the CO is depleted however by a factor of $\sim 10$ for
RU Lup as suggested by Zhang et al. (2020) for disks in Lupus star forming
region – based on the data by Ansdell et al. 2016 and attributed by them to
the coupling of physical and chemical processes – then the true
$M_{\mathrm{gas}}$ would be as high as $\sim 0.1-0.2\,M_{\odot}$. Given that
the stellar mass of RU Lup is $\sim 0.2-1.2\,M_{\odot}$, such a massive disk
would be gravitationally unstable, and we, therefore, consider large depletion
factors unlikely.

In summary, a better understanding of disk physics and evolution could be
achieved by modeling target-by-target and obtaining better-constrained disk
masses for more disks of different ages. The procedure of fitting the long-
wavelength portion of the SED in combination with the $\mathrm{C^{18}O}$ line
emission demonstrated in this work could be easily implemented on other
targets with sufficient photometric data. The DiskMINT code is also released
(Deng, 2023) and available in the public
repository111https://github.com/DingshanDeng/DiskMINT, so that the community
can extend this approach to other disks.

#### Acknowledgments

The authors thank C.P. Dullemond for helpful discussions and assistance on
building our wrapper based on RADMC-3D, thank J. Barnes, A. Youdin and the
anonymous referee for helpful suggestions and comments. DD, IP and UG
acknowledge support from the NASA/XRP research grant 80NSSC20K0273 which made
this work possible. Support for MR’s research was provided by NASA’s Planetary
Science Division Research Program, through ISFM work package ‘The Production
of Astrobiologically Important Organics during Early Planetary System
Formation and Evolution’ at NASA Ames Research Center.

## References

  * Alcalá et al. (2017) Alcalá, J. M., Manara, C. F., Natta, A., et al. 2017, A&A, 600, A20, doi: 10.1051/0004-6361/201629929
  * Allard et al. (2003) Allard, F., Guillot, T., Ludwig, H.-G., et al. 2003, in Brown Dwarfs, ed. E. Martín, Vol. 211, 325
  * Allard et al. (2011) Allard, F., Homeier, D., & Freytag, B. 2011, in Astronomical Society of the Pacific Conference Series, Vol. 448, 16th Cambridge Workshop on Cool Stars, Stellar Systems, and the Sun, ed. C. Johns-Krull, M. K. Browning, & A. A. West, 91.  https://arxiv.org/abs/1011.5405
  * Andrews et al. (2018) Andrews, S. M., Terrell, M., Tripathi, A., et al. 2018, ApJ, 865, 157, doi: 10.3847/1538-4357/aadd9f
  * Andrews et al. (2018) Andrews, S. M., Huang, J., Pérez, L. M., et al. 2018, The Astrophysical Journal Letters, 869, L41, doi: 10.3847/2041-8213/aaf741
  * Ansdell et al. (2016) Ansdell, M., Williams, J. P., van der Marel, N., et al. 2016, The Astrophysical Journal, 828, 46, doi: 10.3847/0004-637X/828/1/46
  * Armitage (2022) Armitage, P. J. 2022, arXiv e-prints, arXiv:2201.07262.  https://arxiv.org/abs/2201.07262
  * Avenhaus et al. (2018) Avenhaus, H., Quanz, S. P., Garufi, A., et al. 2018, The Astrophysical Journal, 863, 44, doi: 10.3847/1538-4357/aab846
  * Bans & Königl (2012) Bans, A., & Königl, A. 2012, ApJ, 758, 100, doi: 10.1088/0004-637X/758/2/100
  * Banzatti et al. (2019) Banzatti, A., Pascucci, I., Edwards, S., et al. 2019, ApJ, 870, 76, doi: 10.3847/1538-4357/aaf1aa
  * Bergin & Williams (2017) Bergin, E. A., & Williams, J. P. 2017, The Determination of Protoplanetary Disk Masses (Cham: Springer International Publishing), doi: 10.1007/978-3-319-60609-5_1
  * Bergin et al. (2013) Bergin, E. A., Cleeves, L. I., Gorti, U., et al. 2013, Nature, 493, 644, doi: 10.1038/nature11805
  * Birnstiel et al. (2011) Birnstiel, T., Ormel, C. W., & Dullemond, C. P. 2011, A&A, 525, A11, doi: 10.1051/0004-6361/201015228
  * Birnstiel et al. (2018) Birnstiel, T., Dullemond, C. P., Zhu, Z., et al. 2018, The Astrophysical Journal, 869, L45, doi: 10.3847/2041-8213/aaf743
  * Bosman & Banzatti (2019) Bosman, A. D., & Banzatti, A. 2019, A&A, 632, L10, doi: 10.1051/0004-6361/201936638
  * Boss & Yorke (1996) Boss, A. P., & Yorke, H. W. 1996, ApJ, 469, 366, doi: 10.1086/177786
  * Brinch & Hogerheijde (2010) Brinch, C., & Hogerheijde, M. R. 2010, Astronomy & Astrophysics, 523, A25, doi: 10.1051/0004-6361/201015333
  * Bruderer (2013) Bruderer, S. 2013, Astronomy and Astrophysics, 559, doi: 10.1051/0004-6361/201321171
  * Bruderer et al. (2012) Bruderer, S., Dishoeck, E. F. V., Doty, S. D., & Herczeg, G. J. 2012, Astronomy and Astrophysics, 541, doi: 10.1051/0004-6361/201118218
  * Comerón (2008) Comerón, F. 2008, in Handbook of Star Forming Regions, Volume II, ed. B. Reipurth, Vol. 5 (ASP Monograph Publications), 295
  * D’Alessio et al. (1998) D’Alessio, P., Cantö, J., Calvet, N., & Lizano, S. 1998, ApJ, 500, 411, doi: 10.1086/305702
  * Davies et al. (2018) Davies, C. L., Kraus, S., Harries, T. J., et al. 2018, The Astrophysical Journal, 866, 23, doi: 10.3847/1538-4357/aade51
  * Deng (2023) Deng, D. 2023, DiskMINT: a tool to estimate disk masses with CO isotopologues, v1.0.0.beta, Zenodo, doi: 10.5281/zenodo.8117966
  * Draine (2003) Draine, B. T. 2003, ApJ, 598, 1026, doi: 10.1086/379123
  * Dullemond et al. (2012) Dullemond, C. P., Juhasz, A. ., Pohl, A. ., et al. 2012, {RADMC}-{3D}: {A} multi-purpose radiative transfer tool - {NASA}/{ADS}.  https://ui.adsabs.harvard.edu/abs/2012ascl.soft02015D/abstract
  * Fang et al. (2018) Fang, M., Pascucci, I., Edwards, S., et al. 2018, ApJ, 868, 28, doi: 10.3847/1538-4357/aae780
  * Fernandes et al. (2018) Fernandes, R. B., Long, Z. C., Pikhartova, M., et al. 2018, The Astronomical Journal, 856, 103, doi: 10.3847/1538-4357/aaaae7
  * France et al. (2014) France, K., Schindhelm, R., Bergin, E., Roueff, E., & Abgrall, H. 2014, ApJ, 784, 127, doi: 10.1088/0004-637X/784/2/127
  * Furlan et al. (2006) Furlan, E., Hartmann, L., Calvet, N., et al. 2006, ApJS, 165, 568, doi: 10.1086/505468
  * Gahm et al. (2013) Gahm, G. F., Stempels, H. C., Walter, F. M., Petrov, P. P., & Herczeg, G. J. 2013, A&A, 560, A57, doi: 10.1051/0004-6361/201322750
  * Gaia Collaboration et al. (2018) Gaia Collaboration, A., Brown, A. G. A., Vallenari, A., et al. 2018, Astronomy and Astrophysics, 616, A1, doi: 10.1051/0004-6361/201833051
  * Gras-Velázquez & Ray (2005) Gras-Velázquez, À., & Ray, T. P. 2005, A&A, 443, 541, doi: 10.1051/0004-6361:20042397
  * Hartmann et al. (1998) Hartmann, L., Calvet, N., Gullbring, E., & D’Alessio, P. 1998, The Astrophysical Journal, 495, 385, doi: 10.1086/305277
  * Hayashi (1981) Hayashi, C. 1981, Progress of Theoretical Physics Supplement, 70, 35, doi: 10.1143/PTPS.70.35
  * Henden et al. (2016) Henden, A. A., Templeton, M., Terrell, D., et al. 2016, VizieR Online Data Catalog, II/336
  * Herczeg et al. (2005) Herczeg, G. J., Walter, F. M., Linsky, J. L., et al. 2005, AJ, 129, 2777, doi: 10.1086/430075
  * Huang et al. (2018) Huang, J., Andrews, S. M., Dullemond, C. P., et al. 2018, ApJ, 869, L42, doi: 10.3847/2041-8213/aaf740
  * Huang et al. (2020) Huang, J., Andrews, S. M., Öberg, K. I., et al. 2020, ApJ, 898, 140, doi: 10.3847/1538-4357/aba1e1
  * Hughes et al. (1994) Hughes, J., Hartigan, P., Krautter, J., & Kelemen, J. 1994, The Astronomical Journal, 108, 1071, doi: 10.1086/117135
  * Ishihara et al. (2010) Ishihara, D., Onaka, T., Kataza, H., et al. 2010, A&A, 514, A1, doi: 10.1051/0004-6361/200913811
  * Kenyon et al. (2008) Kenyon, S. J., Gómez, M., & Whitney, B. A. 2008, in Handbook of Star Forming Regions, Volume I, ed. B. Reipurth, Vol. 4 (ASP Monograph Publications), 405
  * Krijt et al. (2020) Krijt, S., Bosman, A. D., Zhang, K., et al. 2020, ApJ, 899, 134, doi: 10.3847/1538-4357/aba75d
  * Lommen et al. (2009) Lommen, D., Maddison, S. T., Wright, C. M., et al. 2009, A&A, 495, 869, doi: 10.1051/0004-6361:200810999
  * Lommen et al. (2007) Lommen, D., Wright, C. M., Maddison, S. T., et al. 2007, A&A, 462, 211, doi: 10.1051/0004-6361:20066255
  * Long et al. (2017) Long, F., Herczeg, G. J., Pascucci, I., et al. 2017, ApJ, 844, 99, doi: 10.3847/1538-4357/aa78fc
  * McClure et al. (2016) McClure, M. K., Bergin, E. A., Cleeves, L. I., et al. 2016, ApJ, 831, 167, doi: 10.3847/0004-637X/831/2/167
  * Miotello et al. (2014) Miotello, A., Bruderer, S., & Dishoeck, E. F. V. 2014, Astronomy and Astrophysics, 572, doi: 10.1051/0004-6361/201424712
  * Miotello et al. (2022) Miotello, A., Kamp, I., Birnstiel, T., Cleeves, L. I., & Kataoka, A. 2022, Setting the Stage for Planet Formation: Measurements and Implications of the Fundamental Disk Properties, arXiv, doi: 10.48550/ARXIV.2203.09818
  * Miotello et al. (2016) Miotello, A., van Dishoeck, E. F., Kama, M., & Bruderer, S. 2016, Astronomy & Astrophysics, 594, A85, doi: 10.1051/0004-6361/201628159
  * Miotello et al. (2017) Miotello, A., van Dishoeck, E. F., Williams, J. P., et al. 2017, Astronomy & Astrophysics, 599, A113, doi: 10.1051/0004-6361/201629556
  * Paneque-Carreño et al. (2023) Paneque-Carreño, T., Miotello, A., van Dishoeck, E. F., et al. 2023, A&A, 669, A126, doi: 10.1051/0004-6361/202244428
  * Pascucci et al. (2022) Pascucci, I., Cabrit, S., Edwards, S., et al. 2022, arXiv e-prints, arXiv:2203.10068.  https://arxiv.org/abs/2203.10068
  * Pascucci et al. (2016) Pascucci, I., Testi, L., Herczeg, G. J., et al. 2016, ApJ, 831, 125, doi: 10.3847/0004-637X/831/2/125
  * Pascucci et al. (2023) Pascucci, I., Skinner, B. N., Deng, D., et al. 2023, arXiv e-prints, arXiv:2307.02704, doi: 10.48550/arXiv.2307.02704
  * Powell et al. (2022) Powell, D., Gao, P., Murray-Clay, R., & Zhang, X. 2022, Nature Astronomy, 6, 1147, doi: 10.1038/s41550-022-01741-9
  * Ruaud et al. (2022) Ruaud, M., Gorti, U., & Hollenbach, D. J. 2022, The Astrophysical Journal, 925, 49, doi: 10.3847/1538-4357/ac3826
  * Skrutskie et al. (2006) Skrutskie, M. F., Cutri, R. M., Stiening, R., et al. 2006, AJ, 131, 1163, doi: 10.1086/498708
  * Sturm et al. (2022) Sturm, J. A., McClure, M. K., Harsono, D., et al. 2022, A&A, 660, A126, doi: 10.1051/0004-6361/202141860
  * Teague (2019) Teague, R. 2019, The Journal of Open Source Software, 4, 1632, doi: 10.21105/joss.01632
  * Trapman et al. (2021) Trapman, L., Bosman, A. D., Rosotti, G., Hogerheijde, M. R., & van Dishoeck, E. F. 2021, Astronomy & Astrophysics, 649, A95, doi: 10.1051/0004-6361/202039200
  * Trapman et al. (2017) Trapman, L., Miotello, A., Kama, M., Dishoeck, E. F. V., & Bruderer, S. 2017, Astronomy and Astrophysics, 605, A69, doi: 10.1051/0004-6361/201630308
  * Weingartner & Draine (2001) Weingartner, J., & Draine, B. 2001, The Astrophysical Journal, 548, 296, doi: 10.1086/318651
  * Whelan et al. (2021) Whelan, E. T., Pascucci, I., Gorti, U., et al. 2021, ApJ, 913, 43, doi: 10.3847/1538-4357/abf55e
  * Williams & Best (2014) Williams, J. P., & Best, W. M. 2014, Astrophysical Journal, 788, 59, doi: 10.1088/0004-637X/788/1/59
  * Woitke et al. (2016) Woitke, P., Min, M., Pinte, C., et al. 2016, Astronomy and Astrophysics, 586, doi: 10.1051/0004-6361/201526538
  * Woitke et al. (2019) Woitke, P., Kamp, I., Antonellini, S., et al. 2019, Publications of the Astronomical Society of the Pacific, 131, 64301, doi: 10.1088/1538-3873/aaf4e5
  * Wright et al. (2010) Wright, E. L., Eisenhardt, P. R. M., Mainzer, A. K., et al. 2010, AJ, 140, 1868, doi: 10.1088/0004-6256/140/6/1868
  * Yen et al. (2018) Yen, H.-W., Koch, P. M., Manara, C. F., Miotello, A., & Testi, L. 2018, A&A, 616, A100, doi: 10.1051/0004-6361/201732196
  * Zhang et al. (2020) Zhang, K., Schwarz, K. R., & Bergin, E. A. 2020, ApJ, 891, L17, doi: 10.3847/2041-8213/ab7823
  * Zhang et al. (2021) Zhang, K., Booth, A. S., Law, C. J., et al. 2021, The Astrophysical Journal Supplement Series, 257, 5, doi: 10.3847/1538-4365/ac1580