exoALMA. XV. Interpreting the Height of CO Emission Layer

Giovanni P. Rosotti, Cristiano Longarini, Teresa Paneque-Carreño, Gianni Cataldi, Maria Galloway-Sprietsma, Sean M. Andrews, Jaehan Bae, Marcelo Barraza-Alfaro, Myriam Benisty, Pietro Curone, et al.

► To cite this version:

Giovanni P. Rosotti, Cristiano Longarini, Teresa Paneque-Carreño, Gianni Cataldi, Maria Galloway-Sprietsma, et al. . exoALMA. XV. Interpreting the Height of CO Emission Layer. The Astrophysical Journal Letters, 2025, 984, 10.3847/2041-8213/adc42e . insu-05099872

HAL Id: insu-05099872

https://insu.hal.science/insu-05099872v1

Submitted on 6 Jun 2025

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
OPEN ACCESS

<figure><img src="image_2.png" /></figure>

CrossMark

exoALMA. XV. Interpreting the Height of CO Emission Layer

Giovanni P. Rosotti<fnref n="1" />, Cristiano Longarini<fnref n="2" />,<fnref n="1" />, Teresa Paneque-Carreño<fnref n="3" />,<fnref n="4" />,<fnref n="5" />, Gianni Cataldi<fnref n="6" />, Maria Galloway-Sprietsma<fnref n="7" />, Sean M. Andrews<fnref n="8" />, Jaehan Bae<fnref n="7" />, Marcelo Barraza-Alfaro<fnref n="9" />, Myriam Benisty<fnref n="10" />,<fnref n="11" />, Pietro Curone<fnref n="1" />,<fnref n="12" />, Ian Czekala<fnref n="13" />, Stefano Facchini<fnref n="1" />, Daniele Fasano<fnref n="10" />, Mario Flock<fnref n="11" />, Misato Fukagawa<fnref n="6" />, Himanshi Garg<fnref n="14" />, Cassandra Hall<fnref n="15" />,<fnref n="16" />,<fnref n="17" />, Jane Huang<fnref n="18" />, John D. Ilee<fnref n="19" />, Andrés F. Izquierdo<fnref n="7" />,<fnref n="4" />,<fnref n="3" />,<fnref n="27" />, Kazuhiro Kanagawa<fnref n="20" />, Geoffroy Lesur<fnref n="21" />, Giuseppe Lodato<fnref n="1" />, Ryan A. Loomis<fnref n="22" />, Ryuta Orihara<fnref n="20" />, Christophe Pinte<fnref n="21" />,<fnref n="14" />, Daniel J. Price<fnref n="14" />, Jochen Stadler<fnref n="10" />, Richard Teague<fnref n="9" />, Gaylor Wafflard-Fernandez<fnref n="21" />, Andrew J. Winter<fnref n="10" />,<fnref n="11" />, Lisa Wölfer<fnref n="9" />, Hsi-Wei Yen<fnref n="23" />, Tomohiro C. Yoshida<fnref n="6" />,<fnref n="24" />, and Brianna Zawadzki<fnref n="25" />,<fnref n="26" />

Received 2024 November 28; revised 2025 February 24; accepted 2025 March 13; published 2025 April 28

Abstract

The availability of exquisite data and the development of new analysis techniques have enabled the study of emitting heights in protoplanetary disks. In this paper, we introduce a simple model linking the emitting height of CO to the disk surface density and temperature structure. We then apply the model to measurements of the emitting height and disk temperature conducted as part of exoALMA, integrated with additional legacy measurements from the MAPS Large Programme, to derive CO column densities and surface density profiles (assuming a CO abundance) for a total of 14 disks. A unique feature of the method we introduce to measure surface densities is that it can be applied to optically thick observations, rather than optically thin as conventionally done. While we use our method on a sample of well-studied disks where temperature structures have been derived using two emission lines, we show that reasonably accurate estimates can be obtained also when only one molecular transition is available. With our method, we obtain independent constraints from ¹²CO and ¹³CO, and we find they are in general good agreement using the standard ¹²C/¹³C isotopic ratio. The masses derived from our method are systematically lower compared with the values derived dynamically from the rotation curve if using an interstellar matter (ISM) CO abundance, implying that CO is depleted by a median factor ~20 with respect to the ISM value, in line with other works that find that CO is depleted in protoplanetary disks.

Unified Astronomy Thesaurus concepts: Planet formation (1241); Protoplanetary disks (1300)

1. Introduction

In the last few years, the telescope Atacama Large Millimeter/
submillimeter Array (ALMA) has rapidly transformed the
observational field of protoplanetary disks. Thanks to the
combination of high spatial and spectral resolution and the high
sensitivity, ALMA has allowed us to study in detail disk
kinematics (see C. Pinte et al. 2023 for a review). In parallel, new
analysis tools have been developed (S. Casassus & S. Pérez 2019;

<figure><img src="image_3.png" /></figure>

<footnote n="1">¹ Dipartimento di Fisica, Università degli Studi di Milano, Via Celoria 16, 20133 Milano, Italy; giovanni.rosotti@unimi.it</footnote>

<footnote n="2">²Institute of Astronomy, University of Cambridge, Madingley Rd., CB30HA, Cambridge, UK</footnote>

<footnote n="3">³European Southern Observatory, Karl-Schwarzschild-Str. 2, D-85748 Garching bei München, Germany</footnote>

<footnote n="4">⁴Leiden Observatory, Leiden University, P.O. Box 9513, NL-2300 RA Leiden, The Netherlands</footnote>

<footnote n="5">⁵ Department of Astronomy, University of Michigan, 1085 South University Avenue, Ann Arbor, MI 48109, USA</footnote>

<footnote n="6">⁶National Astronomical Observatory of Japan, 2-21-1 Osawa, Mitaka, Tokyo 181-8588, Japan</footnote>

<footnote n="7">⁷Department of Astronomy, University of Florida, Gainesville, FL 32611, USA</footnote>

<footnote n="8">⁸Center for Astrophysics | Harvard & Smithsonian, Cambridge, MA 02138, USA</footnote>

<footnote n="9">⁹Department of Earth, Atmospheric, and Planetary Sciences, Massachusetts Institute of Technology, Cambridge, MA 02139, USA</footnote>

<footnote n="10">¹⁰Université Côte d'Azur, Observatoire de la Côte d'Azur, CNRS, Laboratoire Lagrange, 06304 Nice, France</footnote>

<footnote n="11">¹¹Max-Planck-Institute for Astronomy (MPIA), Königstuhl 17, 69117 Heidelberg, Germany</footnote>

<footnote n="12">¹²Departamento de Astronomía, Universidad de Chile, Camino El Observatorio 1515, Las Condes, Santiago, Chile</footnote>

<footnote n="13">¹³School of Physics & Astronomy, University of St. Andrews, North Haugh, St. Andrews KY16 9SS, UK</footnote>

<footnote n="14">¹⁴School of Physics and Astronomy, Monash University, Clayton VIC 3800, Australia</footnote>

<footnote n="15">¹⁵Department of Physics and Astronomy, The University of Georgia, Athens, GA 30602, USA</footnote>

<footnote n="16">¹⁶Center for Simulational Physics, The University of Georgia, Athens, GA 30602, USA</footnote>

<footnote n="17">¹⁷Institute for Artificial Intelligence, The University of Georgia, Athens, GA, 30602, USA</footnote>

<footnote n="18">¹⁸Department of Astronomy, Columbia University, 538 W. 120th St., Pupin Hall, New York, NY, USA</footnote>

<footnote n="19">¹⁹School of Physics and Astronomy, University of Leeds, Leeds, LS2 9JT, UK</footnote>

<footnote n="20">²⁰College of Science, Ibaraki University, 2-1-1 Bunkyo, Mito, Ibaraki 310-8512, Japan</footnote>

<footnote n="21">²¹Univ. Grenoble Alpes, CNRS, IPAG, 38000 Grenoble, France</footnote>

<footnote n="22">²²National Radio Astronomy Observatory, 520 Edgemont Rd., Charlottesville, VA 22903, USA</footnote>

<footnote n="23">²³Academia Sinica Institute of Astronomy & Astrophysics, 11F of Astronomy-Mathematics Building, AS/NTU, No.1, Sec. 4, Roosevelt Rd., Taipei 10617, Taiwan</footnote>

<footnote n="24">²⁴Department of Astronomical Science, The Graduate University for Advanced Studies, SOKENDAI, 2-21-1 Osawa, Mitaka, Tokyo 181-8588, Japan</footnote>

<footnote n="25">²⁵Department of Astronomy, Van Vleck Observatory, Wesleyan University, 96 Foss Hill Drive, Middletown, CT 06459, USA</footnote>

<footnote n="26">²⁶Department of Astronomy & Astrophysics, 525 Davey Laboratory, The Pennsylvania State University, University Park, PA 16802, USA</footnote>

<footnote n="27">²⁷ NASA Hubble Fellowship Program Sagan Fellow.</footnote>

<footnote>Original content from this work may be used under the terms
of the [Creative Commons Attribution 4.0 licence](https://creativecommons.org/licenses/by/4.0/). Any further
distribution of this work must maintain attribution to the author(s) and the title
of the work, journal citation and DOI.</footnote>
R. Teague 2019; A. F. Izquierdo et al. 2021) to make the most use of the data.

One application of spatially and spectrally resolved observations is to measure the height of the emission surface for molecular lines (K. A. Rosenfeld et al. 2013; C. Pinte et al. 2018), and the temperature at the emitting surface since the emission is optically thick. These kinds of measurements are now available for a growing sample of disks (C. J. Law et al. 2021, 2022, 2023) and also for multiple molecules (e.g., T. Paneque-Carreño et al. 2023; F. Urbina et al. 2024), while initial efforts focused mostly on the CO molecule on account of its brightness.

This paper is part of the exoALMA Large Programme (R. Teague et al. 2025), which consists in observations designed to perform a kinematic analysis targeting 15 protoplanetary disks. As part of a dedicated paper series, A. Izquierdo et al. (2025) and M. Galloway-Sprietsma et al. (2025) derive emitting heights and temperatures for the disks of the sample with two different techniques. In this work, we adopt the emitting heights derived parametrically by A. Izquierdo et al. (2025) and the temperatures derived by M. Galloway-Sprietsma et al. (2025).

While these measurements are becoming routine, comparatively little effort has been invested to investigate what sets the emitting height and what the measurements tell us about the disk structure and composition. Addressing this issue is the focus of this paper. In particular, we aim to develop a model to quantitatively link a given density and temperature structure with the emitting height. Because the observations put constraints on both the latter and the temperature structure, our aim is then to invert the problem and use the observational constraints to measure the disk density.

The paper is structured as follows. In Section 2, assuming hydrostatic equilibrium in the vertical direction, we construct an analytical model to link the emitting height with the disk surface density, with assumptions on the disk vertical temperature profile. The model depends on a free parameter describing photodissociation, which we calibrate against thermochemical models in Section 3. We then apply the model to the exoALMA sample (integrated with additional disks from MAPS) in Section 4, deriving surface densities for our targets. We discuss the implications of our results in Section 5, in particular with respect to the issues of CO abundance and isotopic ratio, and we finally draw our conclusions in Section 6.

## 2. Model for Emission Height

In this work, we use the Eddington–Barbier approximation (D. Barbier 1943; but see F. Paletou 2018 for an interesting historical perspective) to make the assumption that the CO emission height traces the location where $\tau \simeq 2/3$, where $\tau$ is measured from the observer to the emitting layer. We call $N_\tau$ the column density of CO for which this happens. Unless the disk is viewed nearly edge-on, and using the fact that the disk is vertically thin, the line of sight can be approximated as almost vertical in the disk frame, and therefore probes a single radius. We thus need to estimate the optical depth of the emission line along vertical rays to interpret the observed heights.

In addition, we also need to take into account that in the upper layers of the disk molecules are photodissociated (e.g., Y. Aikawa et al. 2002) by the UV radiation (either emitted from the star or from the interstellar field). This implies that the uppermost layers of the disk do not provide any contribution to the line optical depth, since they contain virtually no CO. Thermochemical models find that the CO abundance rises

quickly over a vertically narrow region in which CO self-shields effectively from UV radiation; therefore, it is a good approximation, and often employed (e.g., C. Qi et al. 2011; J. P. Williams & W. M. J. Best 2014; K. M. Flaherty et al. 2015; C. Toci et al. 2021), to assume that there is no contribution to the line optical depth up to some height $z_{ph}$, the photodissociation layer, while for $z < z_{ph}$ CO has a constant abundance $x_{CO}$. While precise constraints on photodissociation would require detailed thermochemistry to be followed exactly, in this work, we make the simplifying assumption that a fixed vertical column of the molecule $N_{ph}$ is required to self-shield from UV radiation, consistent with what is found in more complex calculations (e.g., E. F. van Dishoeck & J. H. Black 1988; R. Visser et al. 2009). This column is typically found to be somewhere around $10^{16} \text{cm}^{-2}$, but we will compare with thermochemical models in Section 3 to find the best $N_{ph}$ value. Note that, under the assumption that CO self-shields, in this paper we define this quantity as the equivalent column of CO required to self-shield, and not as the total gas column—this is because the latter would depend on CO abundance, whereas $N_{ph}$ we employ does not. Combining the two requirements, we expect the emission layer to be at a height where the column density has a critical value $N_{crit} = N_{ph} + N_\tau$. We now focus in Section 2.1 on the calculation of $N_\tau$ and then discuss in the rest of Section 2 our model for the emitting height and how to turn measurements of the emitting height into measurements of the surface density, taking into account that disks have a vertical density and temperature stratification. We will leave $N_{ph}$ as a free parameter in this section and discuss the procedure we followed to set its value in Section 3.

### 2.1. Line Opacity and Critical Column Probed by the Emitting Layer

We recap here standard radiative transfer formulae (see, e.g., F. F. S. van der Tak et al. 2007; J. J. Condon & S. M. Ransom 2016). Under the local thermodynamic equilibrium (LTE) condition, which thermochemical models show to hold true for low J transitions of the CO molecule, the opacity coefficient$^{28}$ at a frequency $\nu$ of a line with emission frequency $\nu_{ul}$, corresponding to the transition from an upper level with the quantum degeneracy $g_u$ to a lower level with the quantum degeneracy $g_l$, is given by

$$ \kappa(\nu) = \frac{c^2}{8\pi\nu_{ul}^2} A_{ul} \frac{g_u}{g_l} x_l n_{mol} \left[ 1 - \exp\left(-\frac{h\nu_{ul}}{k_b T}\right) \right] \phi(\nu), \quad (1) $$

where $A_{ul}$ is the Einstein coefficient of the transition, $x_l$ is the fraction of molecules in the lower level of the transition, $n_{mol}$ is the volume density of the molecule, $k_b$ is the Boltzmann constant, and $\phi(\nu)$ is the line profile. Because we are in LTE, the fraction of molecules in a given energy level is

$$ x_i = \frac{g_i}{Z(T)} \exp\left(-\frac{E_i}{k_b T}\right), \quad (2) $$

where $E_i$ is the energy of the level, and $Z(T)$ is the partition function. We estimate the partition function using the expansion $Z = 1/3 + k_b T/\hbar B_0$, where $B_0$ is the rotational constant of the molecule. We verified that there is no

<footnote>$^{28}$ We follow here the notation of J. J. Condon & S. M. Ransom (2016) in which $\tau = \int \kappa ds$. Note that in this notation the symbol $\kappa$ is not a mass opacity.</footnote>
measurable difference in our results from explicitly calculating the partition function by summing over all the rotational energy levels in the Leiden molecular database (F. L. Schöier et al. 2005) for the range of temperatures we are interested in ($T > 10$ K).

The last factor to discuss is the line profile. Lines in protoplanetary disks are significantly Doppler shifted by Keplerian rotation, but here, we are interested in the optical depth at a specific point of the disk. We can therefore neglect rotation, as long as the velocity does not change significantly along the line of sight. While disks do have vertical velocity gradients (T. Takeuchi & D. N. C. Lin 2002; G. P. Rosotti et al. 2020a; P. Martire et al. 2024), these are typically small (e.g., much smaller than radial velocity gradients), justifying our assumption. We also note that in this work we are interested in estimating the height of the emitting layer, and not the total optical depth of the line (where the Doppler shift from the upper layers to the midplane, or to the back surface, can be significant). In practice, since disks have a strong vertical density gradient as a result of hydrostatic equilibrium, most of the optical depth of the line is generated in the proximity of the emitting layer. A corollary of this fact is that the integrations we will write in the rest of this section, although formally over the whole line of sight up to the emitting layer, are dominated by the region close to the emission layer that we probe observationally. Assuming no Doppler shift across the spatial length-scale of the layer is therefore a good assumption.

We also neglect turbulent broadening, since even in cases where this has been measured (K. Flaherty et al. 2020, 2024; T. Paneque-Carreño et al. 2024) it is at most comparable to the thermal value, and smaller in all other cases (see G. P. Rosotti 2023 for a summary). Finally, although pressure broadening has been detected in TW Hya (T. C. Yoshida et al. 2022) and C. T. Yoshida et al. (2025) report another detection as part of this paper series, pressure broadening is a discernible contribution only in the wings of the line profile coming from the disk midplane, while here we are concerned with the line center originating in a location where the density is significantly lower. Therefore, in what follows, we only consider thermal broadening; this translates to a Gaussian line profile, which at the line center evaluates to

$$ \phi(\nu)|_{\nu_{ul}} = \frac{c}{\nu_{ul}} \frac{1}{\sqrt{\pi} v_D}, \qquad (3) $$

so that the line profile is normalized to 1. The thermal line width $v_D$ is

$$ v_D = \sqrt{2k_b T / m_{mol}}, \qquad (4) $$

where $m_{mol}$ is the mass of the molecule. To see more in detail why we can neglect turbulent line broadening, consider that the turbulent line width would add in quadrature to the thermal line width. Since the opacity (Equation (1)) is inversely proportional to the line width, only a turbulence width much larger than the thermal width would significantly change the results we will now discuss.

Although these are textbook relations, to illustrate what they imply in the context of protoplanetary disks, we show<sup>29</sup> in Figure 1 $N_τ$, the column at which $τ = 2/3$, assuming an

<figure><img src="image_4.png" /><figcaption>**Figure 1.** $N_τ$ the CO column has $τ = 2/3$ for the $J = 3-2$ and $J = 2-1$ transition of $^{12}$CO and $^{13}$CO, as a function of temperature. The photodissociation column value we show is only a reference, and it will be estimated more precisely in Section 3.</figcaption></figure>

isothermal slab of material. To compute this, we use that the optical depth is $τ = \int \kappa(T)ds$, where $κ(T)$ is given by Equation (1). In the isothermal slab, the expression can be rewritten as $σ(T) \int n_{mol} ds = σ(T) N_{mol}$, where we gathered all the terms except the density of the molecule in a factor $σ(T)$— dimensionally, this factor is an area and can be regarded as the cross section to radiation. This highlights that for an isothermal slab the optical depth depends only on the column of material. Because of the transitions we will later use in this paper, we show the results for both the $J = 3-2$ and $J = 2-1$ transitions for $^{12}$CO and $^{13}$CO. There is very little difference between the two isotopologues, while the difference is more substantial for the two transitions, but we still expect them to trace roughly the same region. We note that $N_τ$ is a strong function of temperature, and accurate estimates of the temperature are therefore critical in order to estimate it correctly. This dependence is driven by several factors in Equation (1) that depend on temperature, namely, thermal broadening, the partition function, and the Boltzmann factor appearing in the fractional occupancy of the levels. In general, for these low $J$ transitions, all these factors tend to decrease the line opacity at higher temperatures, since the molecules are spread over a wide velocity range and more quantum levels. In addition, stimulated emission (the negative exponential factor) further reduces the net opacity, since its relevance is proportional to the difference in the occupation fraction between upper and lower levels, but this becomes progressively smaller at high temperatures. This drives the general trend that the critical column increases with temperature. Because disks have a radial (as well as vertical) temperature gradient, this implies that the emitting layer does not trace a fixed column of material, but it traces different columns of material at different radii. Finally, $N_τ$ is of the same order as $N_{ph}$ (at least for low temperatures), implying that we need to accurately estimate $N_{ph}$ for the model to be accurate.

## 2.2. The “Direct” Problem: Computing the Emitting Height Given a Surface Density Profile

Armed with Equation (1), we are now equipped to solve what we call the “direct” problem, which consists of determining the emitting height given a surface density profile

<footnote n="29"><sup>29</sup> We benchmarked our calculation against RADEX (F. F. S. van der Tak et al. 2007) and confirm it yields the same results, except for a factor ~1.06. This is known and descends from the fact that RADEX assumes a rectangular, rather than Gaussian, line profile.</footnote>
$\Sigma(R)$. This is normally done through radiative transfer models, using either simplified assumptions for photodissociation (e.g., J. P. Williams & W. M. J. Best 2014; C. Pinte et al. 2018; C. Toci et al. 2021) or consistent thermochemistry (e.g., S. Bruderer 2013; P. Woitke et al. 2016). To do this, we first need an expression for the vertical density profile $\rho(z)$. In the commonly used vertically isothermal case, this is simply

$$ \rho = \rho_0 \exp\left(-\frac{z^2}{2H^2}\right), \qquad (5) $$

where $H = c_s/\Omega$ is the disk scale-height, and $\rho_0 = \Sigma/\sqrt{2\pi H}$ is the density at the disk midplane. In the more general vertically stratified case, the vertical density profile is given by (e.g., G. P. Rosotti et al. 2020b; P. Martire et al. 2024)

$$ \rho(z) = \rho_0 \frac{c_{s, \text{midplane}}^2}{c_s^2(z)} \exp\left(-\int_0^z \frac{\Omega_k^2(z')z'dz'}{c_s^2(z')}\right), \qquad (6) $$

where $c_s(z)$ is the vertical profile of the sound speed, and $\Omega_k^2(z) = GM_*/(r^2 + z^2)^{3/2}$ is the Keplerian frequency at a height $z$. Because most of the mass is at the disk midplane, we find (see also M. Flock et al. 2013) that for reasonable choices of the vertical temperature the profile $\rho_0 = \Sigma/\sqrt{2\pi H}$ still holds very well (with a typical accuracy of a few percent), as long as the scale-height is estimated using the temperature value at the disk midplane. In practice, however, we integrate Equation (6) to link $\rho_0$ with $\Sigma$.

Once the vertical density profile is known, we can find the photodissociation height $z_{ph}$, which is given implicitly by

$$ -\int_{+\infty}^{z_{ph}} n_{\text{CO}}(z) dz = N_{\text{ph}}, \qquad (7) $$

where $n_{CO}$ is the number density of CO molecules. To link it with the mass density $\rho$, we introduce a factor $x_{CO}$ to parameterize the CO molecular abundance in the disk; with this substitution:

$$ n_{\text{CO}} = \frac{\rho x_{\text{CO}}}{\mu m_h}, \qquad (8) $$

where $\mu = 2.35$ is the mean molecular weight. We will discuss the value of $x_{CO}$ in Section 2.4. At this point, we can finally integrate the equation of radiative transfer to find the height $z_{em}$ where the line reaches an optical depth of two-thirds:

$$ -\int_{z_{\text{ph}}}^{z_{\text{em}}} \kappa(n_{\text{CO}}, T(z)) \frac{dz}{\cos i} = 2/3, \qquad (9) $$

where $\kappa(T(z))$ is the line opacity we derived in Section 2.1. This is a function of vertical height since it depends on the density of the CO molecule, and a further dependence is introduced by the dependence of the temperature with the vertical coordinate. The factor $\cos i$ (with $i$ the disk inclination) accounts for the inclined line of sight to the observer. To highlight more the dependence on the vertical coordinate, we can also rewrite this as

$$ -\int_{z_{\text{ph}}}^{z_{\text{em}}} \sigma(T(z)) n_{\text{CO}}(z) \frac{dz}{\cos i} = 2/3, \qquad (10) $$

where $\sigma(T)$ is the cross section we introduced in Section 2.1.

For the vertically stratified case, these integrals, and therefore the estimation of the emitting layer, can be computed numerically. For the isothermal case, we can do a further step

and derive an analytical solution; Equation (9) evaluates to

$$ \begin{aligned} N_{\text{crit},12} &= \int_{z_{\text{em},12}}^{\infty} n_{\text{CO}}(z) dz \\ &= \int_{z_{\text{em},12}}^{\infty} \frac{x_{\text{CO}}}{\mu m_h} \rho_0 \exp\left(-\frac{z^2}{2H^2}\right) dz \\ &= \frac{x_{\text{CO}} \rho_0}{\mu m_h} \sqrt{\frac{\pi}{2}} H \text{erfc}\left(\frac{z_{\text{em},12}}{\sqrt{2} H}\right), \end{aligned} \qquad (11) $$

where $N_{crit,12} = N_r \cos i + N_{ph}$ is the critical column for $^{12}$CO, and erfc is the complementary error function, defined as

$$ \text{erfc}(x) = \sqrt{\frac{2}{\pi}} \int_{x}^{\infty} e^{-t^2} dt. \qquad (12) $$

Equation (11) is an implicit equation for the emitting height; the solution can be written formally as

$$ z_{\text{em},12} = \sqrt{2} H \text{erfc}^{-1} \left( \frac{\mu m_H}{x_{\text{CO}} \rho_0} \sqrt{\frac{2}{\pi}} \frac{N_{\text{crit},12}}{H} \right). \qquad (13) $$

Since $\text{erfc}^{-1}$ is a decreasing function, the equation shows that, as we may naively expect, the emitting height increases both with the midplane density and with the scale-height $H$. Finally, we note that $\rho H \propto \Sigma$, where $\Sigma$ is the disk surface density, and we thus expect the emitting height to increase also with the surface density, and correspondingly the disk mass.

Regardless of the vertical density profile, Equation (9) has the self-shielding column $N_{ph}$ and the CO abundance $x_{CO}$ as free parameters. We will calibrate the former against thermochemical models in Section 3 and discuss the latter in Section 2.4, as well as constrain it from observations in Section 5.3.

## 2.3. The “Inverse” Problem: Retrieving the Surface Density Profile Given an Observed Emitting Height

The “direct” problem we exposed above illustrates the usefulness of our model, since it allows determining the emitting height with simple calculations rather than setting up a full thermochemical calculation, provided that suitable parameterizations for the temperature and CO abundance be available. However, when interpreting observations, one is normally in the opposite situation: the observations provide an emitting height, and we are interested in knowing what constraints it gives on the surface density profile of the disk. We call this the “inverse” problem. We can regard this problem as an extrapolation from the conditions at the emitting height down to the disk midplane. We detail in what follows how we solve the problem, for the isothermal and vertically stratified case.

### 2.3.1. The Isothermal Case

We first note that the inversion is relatively trivial for the isothermal case (Equation (11)). This yields

$$ \rho_0 = \frac{\mu m_H}{x_{\text{CO}}} \sqrt{\frac{2}{\pi}} \frac{N_{\text{crit},12}}{H} \left( \text{erfc} \left( \frac{z_{\text{em},12}}{\sqrt{2} H} \right) \right)^{-1}. \qquad (14) $$

We can write an analogous expression for $^{13}$CO:

$$ \rho_0 = \frac{\mu m_H}{x_{\text{CO}}} \sqrt{\frac{2}{\pi}} \frac{N_{\text{crit},13}}{H} \left( \text{erfc} \left( \frac{z_{\text{em},13}}{\sqrt{2} H} \right) \right)^{-1}, \qquad (15) $$
where we have assumed that ¹³CO is a factor ¹²C/¹³C ≈ 77 (T. L. Wilson & R. Rood 1994) less abundant, and we have allowed for the general possibility that $N_{crit,13}$ is different from $N_{crit,12}$. Physically speaking, this arises since the two molecules have slightly different parameters for the transition. We will also use this difference to account for the fact that observationally the two molecules have different brightness temperatures,<fnref n="30" /> which leads to different values of $N_τ$ for the two isotopologues. We note that recently T. Paneque-Carreño et al. (2023) presented a similar analysis. In that case, however, the authors assumed a surface density profile and used the constraint provided by the observed emitting height of ¹²CO to derive the disk scale-height, while in this work we use the emitting height to set constraints on the surface density profile.

The two expressions can be evaluated provided that we know the scale-height of the disk H. Indeed, the scale-height is needed since, as we explained when introducing the “inverse” problem, we need to extrapolate from the emitting height to the midplane. We thus need more information in the form of independent constraints on the disk temperature. For the isothermal model, we discuss two different approaches to retrieve this information. The first consists in using the line brightness temperature under the assumption that the emission is optically thick, which is often used in disks to measure the disk temperature.

To introduce the second approach, we note that Equations (14) and (15) constitute a system of two equations and two unknowns ($ρ_0$ and H); it is therefore possible to solve for the two unknowns and determine H in this way. To this end, it is useful to take the ratio between Equation (11) and the equivalent for ¹³CO, which yields a relatively simple equation for H:

$$ \frac{^{12}\text{C} \ N_{\text{crit},13}}{^{13}\text{C} \ N_{\text{crit},12}} \text{erfc}\left(\frac{z_{\text{em},12}}{\sqrt{2}H}\right) = \text{erfc}\left(\frac{z_{\text{em},13}}{\sqrt{2}H}\right) \quad (16) $$

that we can apply directly to the data, if we have measurements of both emitting layers. The equation has a straightforward interpretation as it is merely saying that the column at the ¹³CO emitting layer must be a factor $\frac{^{12}\text{C} \ N_{\text{crit},13}}{^{13}\text{C} \ N_{\text{crit},12}}$ higher than the column at the ¹²CO emitting layer. Interestingly, all other quantities cancel out, and therefore, H depends only on the isotopic ratio. The equation needs to be solved numerically, which can be done using a standard root finding algorithm (we use the Brent solver as implemented in the scipy function brentq). Once H is known, we can apply Equation (14) to retrieve the density in the midplane.

### 2.3.2. The Vertically Stratified Case

We first note that this case can be solved only if we know the vertical temperature profile $T(z)$. This is because the temperature potentially takes a different value at every vertical coordinate z, and for this reason, we can no longer describe the temperature with a single parameter. In the context of exoALMA, the radial and vertical temperature profiles are fitted to the observed brightness temperatures in the work of M. Galloway-Sprietsma et al. (2025). While we refer to that

work for the details of the temperatures we employ, for the remainder of this section, we will simply assume the temperature profile is known.

A fixed temperature profile implies, solving the hydrostatic equilibrium equation, that the density profile is also fixed, save for a normalization factor. Determining this factor is our goal. It is convenient to write the vertical density profile as

$$ \rho(z) = \rho_0 f(z) \quad (17) $$

introducing the value of the density in the midplane $ρ_0$ and encoding the vertical dependency of the density in the function $f(z)$. It is trivial to define along the same lines the CO number density in the midplane $n_0$. With this definition, the photodissociation condition (Equation (7)) becomes

$$ -n_0 \int_{+\infty}^{z_{\text{ph}}} f(z) dz = N_{\text{ph}} \quad (18) $$

while the optical depth condition (Equation (9), which we rewrite using the notation of Equation (10)) becomes

$$ -n_0 \int_{z_{\text{ph}}}^{z_{\text{em}}} \sigma(T(z)) f(z) \frac{dz}{\cos i} = 2/3. \quad (19) $$

For computational reasons, we rewrite this condition as

$$ -\frac{n_0}{\cos i} \left( \int_{+\infty}^{z_{\text{em}}} \sigma(T) f(z) dz - \int_{+\infty}^{z_{\text{ph}}} \sigma(T) f(z) dz \right) = 2/3. \quad (20) $$

Practically speaking, the two conditions contain integrals that we can tabulate once for each radius, irrespective of the value of $n_0$; this is thus a problem that can be solved very efficiently, with typical run times shorter than 1 s on a standard desktop CPU to solve the inverse problem for a given disk. Operationally, we use a root finding algorithm (we use the Brent solver as for the isothermal case) to find the root of Equation (20) in terms of $n_0$. This requires first, for a given choice of $n_0$, to solve Equation (18) for $z_{ph}$, since it is needed in Equation (20). Since we tabulate the integral appearing in Equation (18), this can be done efficiently using linear interpolation. Finally, we note that the Brent solver requires specifying a lower and upper limit bracketing the root. In order to do this, we initially neglect the photodissociation column $N_{ph}$ (i.e., we assume $z_{ph} = \infty$), which provides a reasonable first guess for $ρ_0$ by inverting Equation (20). We use values a factor 1000 above and below this initial estimate to provide the search range; we have verified this range is sufficiently broad to always encompass the actual root.<fnref n="31" />

### 2.4. The Role of the Unknown CO Abundance

We remark that strictly speaking the constraints we derive in this work only concern the CO molecule. Our procedure amounts to turn a measurement of the height at which CO emission is optically thick into a measurement of the total CO column. Only under the assumption that there is a constant scaling factor between the total number density and the CO number density, namely, the $x_{CO}$ factor, can we turn this constraint into a constraint on the total surface density when solving the “inverse” problem. It is thus not a surprise that the surface densities we derive in this work are inversely

<footnote n="30">³⁰ Clearly, this observational fact is in contradiction with the assumption that the disk is vertically isothermal. We will assess in the rest of the paper how well the isothermal assumption holds.</footnote>

<footnote n="31">³¹ As expected, the root is always higher than the initial estimate.</footnote>
proportional to $x_{CO}$, as it is clear for example from Equation (14). Unless otherwise specified, from now on, we will use the standard value $x_{CO} = 10^{-4}$, as derived from the interstellar matter (ISM). There is however a lot of evidence that CO is depleted in protoplanetary disks (see, e.g., A. Miotello et al. 2023 for a review). We will discuss in Section 5.3 whether this value is justified for protoplanetary disks.

### 3. Calibrating the Photodissociation Column $N_{ph}$

As mentioned in the previous section, before we can solve the “inverse” problem and derive the disk surface density given the observed emitting height, there remains a free parameter to be estimated: the photodissociation column $N_{ph}$. Determining this parameter is the aim of this section. In order to accomplish this goal, we compare our model, solving the “direct” problem, to the grid of thermochemical models presented in T. Paneque-Carreño et al. (2025).

#### 3.1. DALI Thermochemical Models

While we refer to T. Paneque-Carreño et al. (2025) for details, we briefly describe here the models we employ, which correspond to the fiducial parameters in T. Paneque-Carreño et al. (2025). Our motivation for only considering fiducial parameters is their conclusion that the parameters mostly affecting the emitting height are the disk mass and CO abundance, while the other parameters have a more limited effect. The models have been run using the DALI thermochemical code (S. Bruderer 2013), which solves self-consistently for the (dust and gas) temperature and the chemical abundances of the various molecular species. We particularly stress that for this work we use the chemical network from A. Miotello et al. (2014), which crucially accounts for CO isotope-selective effects (isotope-selective photodissociation and fractionation reactions) in addition to self-shielding and freeze-out. This is of fundamental importance for our work because, due to isotope-selective effects, the relative abundance of $^{12}$CO and $^{13}$CO can in principle differ from the overall isotopic ratio. This is a particular concern close to the photodissociation region, since photodissociation is an isotope-selective process.

The surface density of the disk is assumed to follow a self-similar solution (D. Lynden-Bell & J. E. Pringle 1974):

$$ \Sigma(R) = \Sigma_c \left(\frac{R}{R_c}\right)^{-p} \exp\left[-\left(\frac{R}{R_c}\right)^{2-p}\right], \quad (21) $$

where $\Sigma_c$ is a normalization factor, $R_c$ is a scale radius, and $p$ is a free parameter describing the steepness of the surface density. In this work, we restrict to the case $p=1$ while we consider $R_c = 100$ au, and we consider four different values for $\Sigma_c$, setting it so that the total mass is [0.0001, 0.001, 0.01, 0.1]$M_{\odot}$. We include both a population of small and large grains; the large grains dominate the mass, containing 95% of the dust mass. The total dust-to-gas ratio is 0.01. The small grain population is assumed to be well coupled to the gas, while the large grain population is vertically settled, with a settling parameter (which should be interpreted as the ratio between the dust and gas scale-height, see S. Bruderer 2013) $\chi = 0.2$. As illuminating spectra, we consider both a T Tauri star with $1 L_{\odot}$ and a Herbig star with $17 L_{\odot}$. The disk structure is computed

solving for hydrostatic equilibrium; an iterative procedure (L. Trapman et al. 2017; T. Paneque-Carreño et al. 2025) is performed to ensure that the density and temperature structures are consistent with each other. The typical run time for the complete iteration process varies between 6 and 10 days for each full disk model, depending on the disk mass, underscoring why a simpler, but still effective, model is desirable to interpret observations. Finally, for each model, we ray trace the $^{12}$CO and $^{13}$CO emission and compute the height where the vertical optical depth at the line center is 2/3.

#### 3.2. Comparison between the Emitting Heights from DALI and from the “Direct” Problem

We now compare the emitting heights between DALI and the solution of our direct problem. To make sure this is a fair comparison, we read the temperature computed by the DALI model and use it in our approach. We also use the DALI total mass density structure, which already accounts for hydrostatic equilibrium, rather than solving for it explicitly.32 Note this can correspond to a different CO volume density, since in our approach the CO abundance is parameterized, while in DALI it is computed self-consistently from the chemical network.

We show in the left panel of Figure 2 an example with the results of the comparison for the case with $M_{disk} = 0.01 M_{\odot}$ for $^{12}$CO around a T Tauri star. The figure shows the emitting height as a function of the radius for the DALI model and for our model for different values of the free parameter $N_{ph}$. The figure shows that in general our model matches well the results of the more accurate, but also significantly more computationally expensive, DALI model, and therefore serves also as a validation that our model, while simplified, correctly captures the most important aspects of the problem. However, the exact result depends on the assumed value of the self-shielding column. As expected, increasing this value lowers the emitting height, since the photodissociation layer becomes thicker. A simple inspection of the figure leads to a value of $N_{ph} \sim 3-4 \times 10^{16} \text{ cm}^{-2}$ as the one minimizing the difference between our model and DALI. To be more quantitative, we computed the mean squared deviation between our model and DALI model and looked for the value of $N_{ph}$ that minimizes it. For this particular case, we settle on a value of $N_{ph} = 4 \times 10^{16} \text{ cm}^{-2}$. We repeated the analysis for the other values of the disk mass, finding in general similar values, although we do find that $N_{ph} = 3 \times 10^{16} \text{ cm}^{-2}$ is a better fit for the case with $M_{disk} = 0.0001 M_{\odot}$. In a related problem, determining the radius beyond which CO is photodissociated, L. Trapman et al. (2023) also found that the critical gas column to shield against photodissociation gets smaller for lower disk masses. Given the relatively modest change in this case, we ignore this difference in what follows.

The case of $^{13}$CO deserves a separate discussion. Simple tests show that using the same value we used for $^{12}$CO leads to severely underestimating the emitting height. We show this case in the middle panel of Figure 2 as the blue line; the model is the same as in the left panel, but we now consider the $^{13}$CO emitting layer. This assumption corresponds to considering

<footnote n="32">32 We verified that using our solver for hydrostatic equilibrium we get a similar density profile to DALI. Some minor difference arises because of the DALI iterative procedure, which means that the density structure is computed solving hydrostatic equilibrium with the temperature computed at the previous step. This effect however is small and negligible for our purposes.</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Left panel: Emitting height as a function of radius for ¹²CO, assuming a disk mass of 0.01M☉. We show the values from the DALI thermochemical model, as well as from our semianalytical model, for different values of the free parameter $N_{ph}$. Central panel: Emitting height as a function of radius for ¹³CO, assuming a disk mass of 0.01M☉ for different values of $N_{ph}$. We compare the DALI data (black line) with the emitting layer predicted assuming $N_{ph,13} = N_{ph,12}$ (blue line) and $N_{ph,13} = N_{ph,12C}/^{12}\text{C}$ (orange line). Right panel: Emitting height as a function of radius for ¹²CO, assuming a disk mass of 0.01M☉ orbiting a Herbig star. We show the values from the DALI thermochemical model, as well as from our semianalytical model, for different values of the free parameter $N_{ph}$.</figcaption></figure>

only shielding by ¹³CO itself. We find instead that, for ¹³CO rescaling $N_{ph}$ with the isotopic abundance, employing a value $N_{ph,13} = N_{ph,12C}^{13}/^{12}\text{C}$ leads to a much better agreement with the DALI models, as shown by the orange line in the middle panel of Figure 2. This implies that in practice isotope-selective photodissociation is not particularly effective for ¹³CO, a finding in line with previous dedicated studies (e.g., R. Visser et al. 2009; A. Miotello et al. 2014), which they attributed to isotope-exchange reactions.

Since we have DALI models with a Herbig star as the central source, we repeated the test to verify whether the optimal value of $N_{ph}$ changes significantly compared to the T Tauri case, as it could be expected from the fact that Herbig stars have a significantly more intense UV field, potentially leading to enhanced CO photodissociation. We show our results in the right panel of Figure 2. We find that the value of $N_{ph}$ minimizing the discrepancy between our model and DALI remains consistent with the T Tauri scenario, with a slight increase to $6 \times 10^{16} \text{cm}^{-2}$. We use this value for consistency, although we note that in practice the difference is small.

We expect that other parameters that we did not explore here, such as for example $f_{large}$, the fraction of large grains comprising the dust, could have an effect on $N_{ph}$ comparable to the difference between T Tauris and Herbigs, since they affect the amount of UV extinction by the dust. Although it is a slightly different problem, we note that also L. Trapman et al. (2023) found a weak (a factor ~2) dependence of the CO dissociation column on this parameter. A detailed exploration of how $N_{ph}$ depends on the microphysics is however outside the scope of this paper.

## 4. Application to Observations

In this section, we describe how we apply the model we presented in the previous sections to interpret measurements of the emission height of protoplanetary disks. As part of the exoALMA paper series, in this paper, we gather most of our sample from the exoALMA sample. In addition, we consider also the disks included in the MAPS (K. I. Öberg et al. 2021) program since they also have high quality kinematical data available. For what concerns the exoALMA sample (see R. Teague et al. 2025 for an overview), measurements of the emitting height are available using two different techniques, DISCMINER(A. Izquierdo et al. 2025) and *disksurf*

(M. Galloway-Sprietsma et al. 2025). The latter paper presents a comparison between the two measurements and a benchmark based on radiative transfer models. In general, the two techniques yield similar results, but nevertheless, the differences are notable. Based on the results of the comparison with radiative transfer models, which show that DISCMINER typically reconstructs the emitting height with higher accuracy at high inclinations due to less confusion with the backside, we opt to use the DISCMINER surfaces. This has the disadvantage, though, that we must work with parameterized surfaces. We show in Appendix A how the results change when using *disksurf*.

For what concerns exoALMA, we are not able to apply our model to the whole sample. We discarded in particular sources that show strong nonaxisymmetric features (MWC 758, CQ Tau), since in our model we implicitly assume azimuthal symmetry, and low inclination ones (HD 135344B, HD 143006, J1604), for which the emitting layer cannot be reliably extracted. Finally, we exclude AA Tau for which the emitting layer extraction is problematic, most likely due to the high inclination. Therefore, the sample used in this work includes DM Tau, HD 34282, J1615, J1842, J1852, LkCa15, PDS66, SY Cha, and V4046. Our requirements are very similar to the analysis of rotation curves presented in C. Longarini et al. (2025), and our samples mostly overlap, with the exception of AA Tau for which they still present their results (although they caveat they may be inaccurate). For what concerns the MAPS sample (AS 209, MWC 480, HD163296, IM Lup, and GM Aur), for consistency, we use the measurements of the emitting heights described in A. F. Izquierdo et al. (2023) that also used DISCMINER, rather than the surfaces originally reported by the MAPS collaboration (C. J. Law et al. 2021) that were derived with *disksurf*. Combining the exoALMA and MAPS sample, we have in total 14 disks.

As previously mentioned, we need to rely on an independent measurement of the temperature to apply our method. For the MAPS sample, we adopt the profiles reported by C. J. Law et al. (2021), while for the exoALMA sample we adopt the profiles reported by M. Galloway-Sprietsma et al. (2025)—note that the two works assumed a different functional form for the temperature profile. For what concerns ¹³CO, M. Galloway-Sprietsma et al. (2025) reports that in the outer parts of the disk the brightness
temperature can be lower than 20 K. Since at face value this points to the emission no longer being optically thick, we rescale in this case the critical optical depth traced by the emitting layer (see Equation (9)) by the ratio between the brightness temperature and 20 K. We find however this is only a small correction (a 10% factor on average on the total mass). In terms of radial range, we only consider radii larger than 0.3 since the inner part of the disk is affected by beam smearing. For the MAPS sample, C. J. Law et al. (2021) reported radial ranges where their fits are valid, and we therefore apply the further constraint that the radius must be contained inside the reported range. For the exoALMA sample, we use the parameter $R_{\text{out}}$ reported by the DISCMINER fit to set the outer radius. Since $^{13}\text{CO}$ emission is typically less extended than $^{12}\text{CO}$, we are in general able to derive a surface density on a wider radial extent for $^{12}\text{CO}$ than we are for $^{13}\text{CO}$. We find that in most cases (see also Section 5.2) there is negligible mass in the region where we get constraints only from $^{12}\text{CO}$.

A final point to discuss when dealing with observations is which errors to report for our estimates. Since our method is computationally cheap, we can afford to bootstrap the errors on parameters inputs. For the temperatures, we simply use the uncertainties reported by M. Galloway-Sprietsma et al. (2025) and C. J. Law et al. (2021) for their temperature fits. For the emission surfaces, however, we cannot directly use the errors reported by DISCMINER; they are likely to be severely underestimated as a result of the strong systematic choice of a parameterized emission surface, and using them directly, we would obtain a severely underestimated error. Instead, we choose to use a representative uncertainty based on benchmark performed on radiative transfer models (see the appendix of M. Galloway-Sprietsma et al. 2025). We find that 25% is a typical uncertainty, and we therefore use this value in the rest of the paper. Because of this choice, we remark that the uncertainties we show should not be considered as true statistical uncertainties but should instead be considered as a representative. Operationally, we assume the error is normally distributed, and we bootstrap 1000 times to quantify uncertainties. Since the results can vary by orders of magnitude, the values we report in this paper are computed from the average and standard deviation in log space of the resulting 1000 samples.

The two measurements of $^{12}\text{CO}$ and $^{13}\text{CO}$ emitting heights correspond to two different measurements of the disk surface density when using the stratified model —we will discuss their (dis)agreement in Section 5.2. For what concerns the isothermal model, as discussed in Section 2, there are in principle multiple options depending on which temperature we use. Combined with the stratified case, this leads potentially to a large number of surface density measurements, and it needs to be discussed which of these are the most reliable. To simplify the illustration of our results, we now focus mostly on the stratified model, which we consider our best estimate. We discuss in Section 5.1 the results of the isothermal model when we use the line brightness temperature. We anticipate the results in this case are relatively comparable to the stratified case, while we find instead that solving the system of two equations to derive $H$ cannot be always reliably applied. For this reason, we discuss this case in Appendix B.

## 4.1. A Worked Example: LkCa 15

Given the large number of surface density measurements when applying the various combinations in our method to the whole sample, we now choose to focus on a specific example

<figure><img src="image_6.png" /><figcaption>**Figure 3.** Surface density as a function of radius derived with the stratified model for LkCa 15. We show the results for both $^{12}\text{CO}$ and $^{13}\text{CO}$; shaded regions show the uncertainties (see text for a discussion of the radial dependence of the uncertainty). For comparison, we also show the surface density implied by the fit to the rotation curve presented in this paper series in C. Longarini et al. (2025) and a constraint from optically thin $C^{18}\text{O}$ and $C^{17}\text{O}$ lines presented in J. A. Sturm et al. (2023).</figcaption></figure>

that we illustrate in detail. As a representative example, we choose LkCa 15, a large and bright disk that thus constitutes one of the best case scenarios for applying our method; in addition, it also has independent constraints from optically thin lines ($C^{18}\text{O}$ and $C^{17}\text{O}$), presented in J. A. Sturm et al. (2023), that we can compare to. Figure 3 shows the fitted surface density. We also show the surface density profile resulting from the disk mass reported in the accompanying paper by C. Longarini et al. (2025), obtained by fitting the disk rotation curve. We remark that, while the method we employ here makes no assumption about the shape of the surface density profile, C. Longarini et al. (2025) and J. A. Sturm et al. (2023) assumed a self-similar profile.

We first note that the surface density profiles we derive with our method appear to be in a reasonable range—their values are bracketed in most of the radial range by the other two constraints. The implied surface density profile decreases with the radius, as one may naively expect. As already mentioned, our method can be essentially considered as an extrapolation of the density from the emitting height down to the midplane, and as such, it yields large uncertainties; for reference, at a radius of 100 au, the uncertainty is 0.7 dex. The advantage of our method is that it can be applied to optically thick emission, but this comes at the cost of large uncertainties. We also note that the uncertainty decreases with the radius. This is driven by the decreasing trend we described of the surface density with the radius; a higher surface density implies that the emitting height is in a region where the density gradient with the height is steep, causing larger uncertainties in the extrapolation. Finally, we find that the uncertainty is entirely dominated by the uncertainty in the emitting height determination; including only the error on the temperature produces an uncertainty that is not visible on the scale of the plot. This is because the errors reported by M. Galloway-Sprietsma et al. (2025) are in general small (e.g., fractions of a kelvin for the normalizations).

Comparing $^{12}\text{CO}$ and $^{13}\text{CO}$, we note that the values are close to each other and are compatible within the uncertainties for most of the radial range. More quantitatively, the total disk mass derived from $^{12}\text{CO}$ is $0.02M_{\odot} \pm 0.7$ dex (i.e., $0.02^{+0.08}_{-0.08}$)
and $0.06M_{\odot} \pm 1$ dex for $^{13}$CO. Given the uncertainties, the two estimates are in good agreement. The value reported in J. A. Sturm et al. (2023) is slightly lower at $0.01M_{\odot}$. From the figure, it could seem that this estimate should be even lower, but this is because the mass in this case is integrated over the whole domain; integrating over only the radial extent where we get a surface density measurement from $^{12}$CO yields a factor 2 lower mass.

A simple fit with a power-law surface density for $^{12}$CO returns a value of -2.3 for the exponent. This value is steep, but it should be interpreted with caution since a comparable fit to the optically thin constraint yields a similar value of -2.4; this is due to the exponential drop off of the disk surface density in the outer regions.

Our estimates are lower than the dynamical mass ($\sim 0.1M_{\odot}$), and this is significant for $^{12}$CO. Bearing in mind that in order to derive our estimates we had to assume a functional form of the temperature structure, which in principle could affect our results, we interpret the difference as due to CO depletion. Since our constraints are derived assuming the standard ISM CO abundance of $10^{-4}$, the implication is that in this source CO must have a factor $\sim 10$ lower abundance. We will show this conclusion is not an isolated case but holds for the whole sample (and in fact, LkCa 15 is the source in the sample with the smallest amount of depletion).

## 4.2. Surface Density Profiles for the Whole Sample

After describing LkCa 15 in detail, we now show the results for the whole sample in Figure 4. As for LkCa 15, we find that uncertainties are always dominated by the error we attribute to the emitting height determination, rather than the temperature uncertainty. While we will not describe each disk individually, we note that many of the trends we described before still hold. To start with, the surface densities we derive decrease with the radius.33 The two estimates from $^{12}$CO and $^{13}$CO generally agree within the uncertainties, but they are always lower than the estimate from the rotation curve (reported in C. Longarini et al. 2025 for exoALMA; and in P. Martire et al. 2024 for MAPS). In the interpretation we discussed above, this points to CO depletion being a widespread phenomenon. We will quantify the level of depletion in the next section. As for LkCa 15, we remark that the estimates derived from the rotation curve assume a self-similar profile with a fixed slope (but see S. M. Andrews et al. 2024 for MWC 480, where the slope is let free to vary). We also added as visual aids black dashed lines showing $r^{-1}$ and $r^{-2}$ profiles in order to gauge the steepness of the surface density slope. Most of the surface densities we derive follow more closely the steeper profile (see also the discussion in M. Galloway-Sprietsma et al. 2025).

Finally, for the MAPS sample, we also included as the red line the surface density constraints from (partially) optically thin $C^{18}$O (K. Zhang et al. 2021). We take their result on the CO column density and rescale it to a total density with the same assumption we employ to derive our constraints, i.e., a constant CO abundance of $10^{-4}$. We note this is different from the surface density reported in K. Zhang et al. (2021), who also employed information from the spectral energy distribution to derive the total disk mass and in this way derive a spatially resolved CO depletion profile. We choose to compare our

results in this way since the raw output of our method is a CO column, and comparing directly the surface density reported by K. Zhang et al. (2021) would not be a fair comparison.

# 5. Discussion

## 5.1. Comparison with the Isothermal Model

As mentioned, we consider the stratified model our best estimate since it uses the full information on how the temperature profile depends on the vertical coordinate. However, this approach requires having data for at least two emission lines. It is conceivable that in some circumstances high quality observations could be available instead for a single emission line, for example because one of the lines (in this case, likely the most abundant, $^{12}$CO) is cloud contaminated, or because a line may be particularly faint (in this case, likely the least abundant, $^{13}$CO), or simply because of requirements in the spectral setup that did not permit to include two. It is interesting in this case to quantify how well the surface density can be reconstructed.

We show the comparison in Figure 5 for the reference case of LkCa 15. For the isothermal model, solid lines are for $^{12}$CO while dashed lines are for $^{13}$CO. As discussed in Section 2.3.1, we show estimates from the isothermal model using the brightness temperature at the emission surface to set the scale-height. In addition, for reference, we include an additional estimate where we set the temperature using the temperature fit of M. Galloway-Sprietsma et al. (2025) in the midplane. This estimate cannot be computed when observations of a single line are available, and it is therefore of limited practical use. However, it is useful for illustrative reasons. First of all, we note that using the emitting height temperature gives an estimate that compares reasonably well with the more complex stratified case—for example, for $^{12}$CO, the difference is roughly a factor 2. When using instead the midplane temperature, we see that it produces a much higher value than the stratified case (note the scale of the y-axis)—to be more quantitative, integrating the surface density for $^{13}$CO leads to a problematic disk mass of $\sim 1M_{\odot}$ when using the midplane temperature, compared to $0.07M_{\odot}$ for the stratified case, and in a similar way, for $^{12}$CO, the difference is a factor 6. We discuss more in detail the reasons why the midplane temperature tends to lead to an overestimate of the surface density, despite naively one expecting it to work better since most of the mass mass is at the midplane, in Appendix C.

Moving to the whole sample, we show in Figure 6 the results for all the disk masses. We do not show error bars on this plot since the various data points come from the same measurements, and the goal is to compare the different techniques with each other. We can see that the trends discussed for LkCa 15 remain also for the whole population: the mass estimated using the midplane temperature always overestimates the temperature of the stratified model, sometimes severely so (the average overestimate is a factor $\sim 6$, but the plot shows in some cases it can be larger than 1 order of magnitude). On the other hand, using the emitting height temperature returns a value that most of the time is close to the stratified case, and in general underestimated. To quantify the discrepancy, we compute that the average ratio between the two masses is 0.8 (the median has a similar value) with a spread of $\sim 0.2$ dex, which we can consider as an empirical uncertainty of the approximation. We conclude that, when only one line is available, using the isothermal model and using the brightness temperature at the emitting surface to set the temperature are satisfactory approximations.

<footnote n="33">33 Although some of our sources are transition disks and potentially have inner gas cavities, the cavity radius is inside the two central beams.</footnote>
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Like Figure 3, but for all the disks in our sample. We plot the constraint from the rotation curve (P. Martire et al. 2024; C. Longarini et al. 2025) as a dashed line when they only constitute upper limits. For the MAPS sample, the red line shows the surface density constraints from (partially) optically thin C$^{18}$O (K. Zhang et al. 2021). Black dashed lines are visual aids showing $r^{-1}$ and $r^{-2}$ profiles.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Like Figure 3, but we now show the results of the stratified model in comparison with the results of the isothermal model. We consider two cases, described in detail in the main text, in which we set the temperature based either on the brightness temperature of the line or on the fitted disk temperature profile in the midplane. For the isothermal model, solid lines are for $^{12}$CO while dashed lines are for $^{13}$CO.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 6.** Masses derived from the isothermal model, assuming either the temperature in the midplane or the temperature at the emitting height, in comparison with the stratified case for the whole sample.</figcaption></figure>

## 5.2. Difference between $^{12}$CO and $^{13}$CO Estimates

To make comparisons across the whole sample, we now refer only to the disk mass. In this subsection, we investigate whether there is a systematic difference between the mass derived from $^{12}$CO and $^{13}$CO. We show the result of the comparison in Figure 7, in which we plot the mass derived from $^{13}$CO as a function of the mass derived from $^{12}$CO. For consistency, the $^{12}$CO mass reported in this plot has been computed only over the radial range over which the $^{13}$CO estimate is defined, although this is a minor difference (a factor 2 on average, with the notable exception of IM Lup for which the difference is significant). The gray solid line represents equal values. The first thing to note is that there are a few outliers, which are significantly far from the line where the two masses are the same. Two of these are HD163296 and GM Aur, for which the $^{12}$CO mass is higher than the $^{13}$CO. An inspection of the surface density reveals that the high $^{12}$CO mass

<figure><img src="image_10.png" /><figcaption>**Figure 7.** Comparison between the mass inferred using the $^{12}$CO emitting height and the mass inferred using the $^{13}$CO emitting height. The gray solid line marks where the two masses are the same. The red lines are the results of the power-law fit described in the main text.</figcaption></figure>

is dominated by the inner regions of the disk where the $^{12}$CO surface density goes to very high values. LkCa 15 is instead in the opposite situation where the $^{13}$CO mass is significantly higher than the $^{12}$CO one. It is likely that in these disks there is an issue in the inner disk with one emission surface; however, the high surface density also implies large error bars, and these disks remain compatible with lying on the equality line. For this reason, we decide to keep them in the sample.

For what concerns the other disks, there is overall good agreement between the two mass estimates, although it appears that on average the $^{13}$CO mass is slightly higher than the $^{12}$CO one. To quantify this, we fitted the values using a power-law relation with the method<fnref n="34" /> developed by B. C. Kelly (2007), which allows us to include errors on both axes. We show the results of the fit with the red lines. The fit confirms that there is a good correlation between the two values, with a slope of $1.1 \pm 0.2$, and also that on average the $^{13}$CO mass is higher than $^{12}$CO, with a value of the normalization $2.5 \pm 0.8$ dex. We note however that this value has a large error, meaning that the two mass estimates remain compatible within the uncertainties. We remark that a discrepancy could point to an isotopic ratio $^{12}$C/$^{13}$C different from the value of ~77 we have assumed here. Indeed, recently, T. C. Yoshida et al. (2022) found evidence in TW Hya, where the almost face-on geometry allows to measure the column of material in the optically thin line wings, of a reduced $^{12}$C/$^{13}$C isotopic ratio when comparing $^{13}$CO to $^{12}$CO (see also E. A. Bergin et al. 2024 for a discussion of C fractionation in this source). They measured a value of 21 at radii between 70 and 110 au. The value of the normalization we derive would correspond to an isotopic ratio of ~30, but, considering the uncertainties in our estimates, we are unable to confirm whether we have evidence

<footnote n="34">³⁴ An open-source implementation is available at the following link: https://github.com/jmeyers314/linmix.</footnote>
<table><caption>**Table 1**

Masses Derived from ¹²CO and ¹³CO in This Work, with Their Uncertainties</caption><thead><tr><th></th><th>¹²CO mass<br>(M<sub>☉</sub>)</th><th>Uncertainty<br>(dex)</th><th>¹³CO Mass<br>(M<sub>☉</sub>)</th><th>Uncertainty<br>(dex)</th></tr></thead><tbody><tr><td>DM Tau</td><td>3.0 × 10<sup>-3</sup></td><td>0.2</td><td>1.9 × 10<sup>-3</sup></td><td>0.1</td></tr><tr><td>HD 34282</td><td>2.2 × 10<sup>-2</sup></td><td>0.6</td><td>1.1 × 10<sup>-2</sup></td><td>0.5</td></tr><tr><td>J1615</td><td>7.2 × 10<sup>-3</sup></td><td>0.7</td><td>3.9 × 10<sup>-3</sup></td><td>0.3</td></tr><tr><td>J1842</td><td>5.7 × 10<sup>-4</sup></td><td>0.3</td><td>5.6 × 10<sup>-4</sup></td><td>0.1</td></tr><tr><td>J1852</td><td>2.7 × 10<sup>-4</sup></td><td>0.2</td><td>4.3 × 10<sup>-4</sup></td><td>0.1</td></tr><tr><td>LkCa 15</td><td>1.7 × 10<sup>-2</sup></td><td>0.7</td><td>6.4 × 10<sup>-2</sup></td><td>1.1</td></tr><tr><td>PDS 66</td><td>8.1 × 10<sup>-5</sup></td><td>0.2</td><td>5.2 × 10<sup>-5</sup></td><td>0.0</td></tr><tr><td>SY Cha</td><td>2.2 × 10<sup>-3</sup></td><td>0.3</td><td>2.0 × 10<sup>-3</sup></td><td>0.2</td></tr><tr><td>V4046 Sgr</td><td>5.5 × 10<sup>-4</sup></td><td>0.2</td><td>9.5 × 10<sup>-4</sup></td><td>0.3</td></tr><tr><td>MWC 480</td><td>2.8 × 10<sup>-3</sup></td><td>0.2</td><td>3.3 × 10<sup>-3</sup></td><td>0.2</td></tr><tr><td>IM Lup</td><td>4.6 × 10<sup>-2</sup></td><td>1.2</td><td>5.2 × 10<sup>-3</sup></td><td>0.4</td></tr><tr><td>GM Aur</td><td>2.5 × 10<sup>-2</sup></td><td>0.9</td><td>2.9 × 10<sup>-3</sup></td><td>0.4</td></tr><tr><td>HD 163296</td><td>1.9 × 10<sup>-1</sup></td><td>1.3</td><td>9.7 × 10<sup>-3</sup></td><td>0.5</td></tr><tr><td>AS 209</td><td>4.3 × 10<sup>-4</sup></td><td>0.3</td><td>1.5 × 10<sup>-4</sup></td><td>0.1</td></tr></tbody></table>

of the same behavior in our sample. It is worth remembering that our uncertainties are not true statistical uncertainties and could be overestimated or underestimated, but since we regard them as typical, we will not speculate further on the isotopic ratio. Further work is needed to either reduce the error bars or increase the sample to come to firmer conclusions.

Overall, the fact that the masses derived from the two isotopologues are in such a good agreement is remarkable considering that our method relies on knowing the disk temperature structure and that we are estimating the mass from optically thick emission, although as we highlighted this comes at the cost of large error bars. We also note that the measurements of the two emitting heights are separate, so there is no guarantee a priori that the two values be similar.

## 5.3. Comparison with Other Mass Estimates and CO Depletion

We now compare the masses resulting from our method, which we list in Table 1, with the dynamical mass estimates from the accompanying work of C. Longarini et al. (2025). In addition, we also consider the estimates based on C¹⁸O for the MAPS disks (K. Zhang et al. 2021), which we already showed in Section 4.2. The result of the comparison is shown in Figure 8.

The first thing to note is that our mass estimates are systematically lower than the dynamical masses, pointing strongly toward CO depletion as a likely explanation. In addition, it is reassuring that the sources where C. Longarini et al. (2025) reports only upper limits on the disk mass (denoted as arrows in the figure) correspond to the lowest masses in the sample also with our methodology.

For what concerns the comparison with the optically thin derived mass, we note first that the strongest discrepancy is MWC480. This discrepancy is driven by radii smaller than 100 au, where the surface density reported by K. Zhang et al. (2021) increases steeply—outside this radius, there is reasonable agreement with our results. The origin of this discrepancy is unclear. A possible explanation could be a high level of CO depletion spatially localized in the upper layers of the inner disk, explaining why other regions of the disk are in better agreement. It is worth noting that in the inner disk the rotation curve estimate (P. Martire et al. 2024; and see also S. M. Andrews et al. 2024 for a more flexible parameterization that gives similar results) is comparable to the results of K. Zhang et al. (2021), possibly implying no CO depletion at those radii.

Even considering MWC 480, in the MAPS sample, the ¹³CO mass is on average a factor 1.25 higher than the C¹⁸O estimate. Such an agreement is remarkable considering our estimates are based on optically thick emission. For what concerns ¹²CO, when taking into account uncertainties, there is also good agreement at face value, with an average of 1.3. Looking however at individual sources tells a different story, with GM Aur and HD163296 (which as mentioned in Section 5.2 are among the largest discrepancies between ¹²CO and ¹³CO) being severely overestimated (but with large error bars) and MWC 480 that we already discussed. As mentioned in the previous section, for the first two, the discrepancy is driven by the ¹²CO surface density reaching very high values in the inner (≤100 au) disk.

We now consider more in detail what levels of CO depletion are implied by our results. We show in Figure 9 the ratio between the dynamical mass estimates and the ¹³CO masses (using this or the ¹²CO one is arbitrary since we showed they are statistically the same). Most of the data points cluster around a value of a few tens. Excluding the disks that have only upper limits from the rotation curve, the median CO depletion factor is ~20, and the average (taking into account error bars) CO depletion factor is ~50. This value is in line with what was found recently by T. Paneque-Carreño et al. (2025) using a similar methodology to what we employ here. More broadly, our results confirm the mounting evidence (e.g., M. K. McClure et al. 2016; E. A. Bergin & J. P. Williams 2017; A. Miotello et al. 2017) that CO appears heavily depleted in protoplanetary disks. The sample we present here is currently the largest sample for which the mass is robustly constrained from the rotation curve.

The depletion of CO is probed observationally as part of this paper series also by another study (L. Trapman et al. 2025). In that case, the emission lines used probe deeper into the disk midplane, and are therefore more sensitive to the levels of depletion in the midplane rather than in the disk upper layers. We refer to the paper for an extensive discussion of their findings.

## 5.4. Caveats and Future Prospects

As repeatedly stated, our surface density estimates are extrapolations from the emitting height down to the midplane, and as such, they require knowledge of the disk temperature profile. In this work, we have employed a parameterization of the disk temperature motivated from previous works, but it is likely that real disks have a more complex temperature structure. Our error analysis has shown that errors in the emitting height dominate over errors in temperature, but we reached this conclusion using only the statistical errors on the temperature; it is considerably more complex to estimate the uncertainty due to the chosen temperature parameterization. In addition, while we showed that the isothermal model produces similar results to our best estimate when using the line brightness temperature, this is only a validation of the former against the latter, and in principle, we cannot exclude that both estimates could be off from the correct value. If in the future more flexible, or even nonparametric, temperature structures become available, it would be important to assess their impact on our results.

Another aspect where our estimates could be improved consists in the determination of the emitting heights, since the two methods employed in exoALMA can sometimes have notable differences (see the appendix of M. Galloway-Sprietsma et al. 2025; and Appendix A). As shown, our uncertainties are currently dominated by the emitting height determination, and improvement
<figure><img src="image_11.png" /><figcaption>**Figure 8.** Comparison between the masses measured with various methods, namely, the dynamical masses derived from the rotation curve (P. Martire et al. 2024; C. Longarini et al. 2025), the emitting heights (this work), and (moderately) optically thin $C^{18}O$ (K. Zhang et al. 2021) for the disks in our sample.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 9.** CO depletion factors needed with respect to the ISM value $x_{CO} = 10^{-4}$ to match the masses derived in this work from $^{13}$CO with the kinematical masses derived in C. Longarini et al. (2025) and P. Martire et al. (2024).</figcaption></figure>

in the accuracy of this estimate is therefore an obvious point where to make progress.

In this work, we use the $J=3-2$ transition from the exoALMA sample and the $J=2-1$ transition from the MAPS sample for $^{12}$CO and $^{13}$CO. Those two samples currently constitute the best available observations, but we note that in the literature there are already estimates (C. J. Law et al. 2022, 2023, 2024; T. Paneque-Carreño et al. 2023; L. M. Stapper et al. 2023) of the emitting height for additional sources in CO molecular emission, albeit at lower data quality. Only in some of these cases, there are also temperature fits available, and the extraction procedure of the emitting height is not homogeneous across the literature — further work would be required before we can apply our method to these sources, but it should be possible. Considering there are ongoing ALMA Large Programs to survey gas emission in large samples (e.g., DECO and CHEER) and that they should be able to measure the emitting height at least for the brightest and largest disks, it is likely the observational sample our technique can be applied to will expand in the future. In addition, more isotopologues and transitions are observable with ALMA, and in the future, this work could be extended in this direction. In particular, multiple transitions would allow a more precise measurement of the vertical temperature profile, in a similar fashion to what was done in V. Pezzotta et al. (2025) for the measurement of the rotation curve, rather than using only two as we employ here. We also remark that it should be especially useful to consider high-$J$ transitions that have a significantly different opacity from those we considered here. Different transitions of the

same molecule, with an accurate measurement of the temperature profile, would also allow us to study how the abundance of the molecule changes at different heights, along similar lines to what we did here when comparing the surface density inferred from $^{12}$CO and $^{13}$CO to investigate the isotopic ratio.

In this work, we have derived the CO depletion factor comparing our results with the dynamical mass. In the last few years, chemical methods have been devised (D. E. Anderson et al. 2019, 2022; L. Trapman et al. 2022) to infer the same factor using observations of the same disk in CO isotopologues and $N_2H^+$. A comparison between the chemical methods and the dynamical masses is presented as part of this paper series in L. Trapman et al. (2025). They find lower levels of depletion than what we report here, and suggest this may be because their method is probing deeper layers in the disk, where CO may be less depleted.

Finally, in this work, we focused on the CO molecule. We note however that our model could be easily extended to other optically thick molecules. The only requirement is to have a simple model describing how the abundance of the molecule changes throughout the disk. If this is not available, thermochemical models remain the only way to link the observed emitting heights with the disk's physical structure.

## 6. Conclusions

In this paper, we introduced a semianalytical model linking the disk density and temperature structure with the emitting height of the CO molecule. The model depends on a free
parameter describing photodissociation that we calibrated
against thermochemical calculations.

We then applied our model to the MAPS and exoALMA
disk samples where the CO emitting height has been measured.
Our model turns these measurements into measurements of the
CO column and therefore, assuming a CO abundance, into a
disk surface density. By integrating this value, we also provide
measurements of the disk mass.

The key findings of our work are as follows:

1. In this work, we provide a method to measure the disk surface density from optically thick, rather than thin as commonly done, emission. This comes at the cost of increased uncertainty. We find that the dominant source of error in our measurements is the error on the emitting height of the transition.

2. In the sample, we consider the temperature has been determined using two different, optically thick transitions. We show that when only one transition is available the use of a vertically isothermal model still yields satisfactory constraints on the disk surface density if using the line brightness temperature to set the disk temperature.

3. Our model yields two different constraints for the ¹²CO and the ¹³CO molecule. In general, they are in relatively good agreement assuming a standard isotopic ratio, which is a confirmation that our results are robust. The ¹³CO masses are systematically higher by a factor ~2.5, in line with recent findings in other works of variations of the isotopic ratio in protoplanetary disks, but given the large uncertainties, we are unable to further investigate this difference.

4. We compared the masses we obtain with the dynamical fits to the disk rotation curve. Our values are in general lower, implying that the CO abundance needs to be reduced with respect to the standard ISM value we assumed. A median CO depletion of a factor ~20 is needed to reconcile our mass measurements with the dynamical values, in line with other results showing that CO is depleted in protoplanetary disks.

Acknowledgments

We thank an anonymous referee for comments that improved
the clarity of the paper. This paper makes use of the following
ALMA data: ADS/JAO.ALMA#2021.1.01123.L. ALMA is a
partnership of ESO (representing its member states), NSF
(USA), and NINS (Japan), together with NRC (Canada),
MOST and ASIAA (Taiwan), and KASI (Republic of Korea),
in cooperation with the Republic of Chile. The Joint ALMA
Observatory is operated by ESO, AUI/NRAO, and NAOJ. The
National Radio Astronomy Observatory is a facility of the
National Science Foundation operated under a cooperative
agreement by Associated Universities, Inc. We thank the North
American ALMA Science Center (NAASC) for their generous
support including providing computing facilities and financial
support for student attendance at workshops and publications.

J.B. acknowledges support from NASA XRP grant No. 80NSSC23K1312. M.B., D.F., and J.S. have received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (PROTOPLANETS, grant agreement No. 101002188). Computations by J.S. have been performed on the "Mesocentre SIGAMM" machine, hosted by Observatoire de la Cote d'Azur. P.C.

acknowledges support by the Italian Ministero dell'Istruzione,
Università e Ricerca through the grant Progetti Premiali 2012 -
iALMA (CUP C52113000140001) and by the ANID BASAL
project FB210003. S.F. is funded by the European Union (ERC,
UNVEIL, 101076613), and acknowledges financial contribution
from PRIN-MUR 2022YP5ACE. M.F. is supported by a Grant-in-
Aid from the Japan Society for the Promotion of Science
(KAKENHI: No. JP22H01274). J.D.I. acknowledges support from
an STFC Ernest Rutherford Fellowship (ST/W004119/1) and a
University Academic Fellowship from the University of Leeds.
Support for A.F.I. was provided by NASA through the NASA
Hubble Fellowship grant No. HST-HF2-51532.001-A awarded by
the Space Telescope Science Institute, which is operated by the
Association of Universities for Research in Astronomy, Inc., for
NASA, under contract NAS5-26555. C.L. has received funding
from the European Union's Horizon 2020 research and innovation
program under the Marie Skłodowska-Curie grant agreement No.
823823 (DUSTBUSTERS) and by the UK Science and
Technology research Council (STFC) via the consolidated grant
ST/W000997/1. C.P. acknowledges Australian Research Council
funding via FT170100040, DP18010423, DP220103767, and
DP240103290. D.P. acknowledges Australian Research Council
funding via DP18010423, DP220103767, and DP240103290. G.
R. acknowledges funding from the Fondazione Cariplo, grant No.
2022-1217, and the European Research Council (ERC) under the
European Union's Horizon Europe Research & Innovation
Programme under grant agreement No. 101039651 (DiscEvol).
H.-W.Y. acknowledges support from National Science and
Technology Council (NSTC) in Taiwan through grant NSTC
113-2112-M-001-035- and from the Academia Sinica Career
Development Award (AS-CDA-111-M03). G.W.F. acknowledges
support from the European Research Council (ERC) under the
European Union Horizon 2020 research and innovation program
(grant agreement No. 815559 (MHDiscs)). G.W.F. was granted
access to the HPC resources of IDRIS under the allocation
A0120402231 made by GENCI. Support for B.Z. was provided by
The Brinson Foundation. C.H. acknowledges support from NSF
AAG grant No. 2407679. G.L. has received funding from the
European Union's Horizon 2020 research and innovation program
under the Marie Skłodowska-Curie grant agreement No. 823823
(DUSTBUSTERS). T.C.Y. acknowledges support by Grant-in-Aid
for JSPS Fellows JP23KJ1008. A.J.W. has received funding from
the European Union's Horizon 2020 research and innovation
program under the Marie Skłodowska-Curie grant agreement No
101104656. Views and opinions expressed by ERC-funded
scientists are however those of the author(s) only and do not
necessarily reflect those of the European Union or the European
Research Council. Neither the European Union nor the granting
authority can be held responsible for them. T.P.-C. thanks the
Heising-Simons Foundation for their funding through the 51
Pegasi b Postdoctoral Fellowship. The National Radio Astronomy
Observatory and Green Bank Observatory are facilities of the U.S.
National Science Foundation operated under cooperative agree-
ment by Associated Universities, Inc.

Facility: ALMA.

Software: DALI.

**Appendix A**

Comparison between *disksurf* and DISCMINER Emitting
Heights

We show in Figure 10 a comparison between the disk masses
obtained using the emitting heights from *disksurf* and those
<figure><img src="image_13.png" /><figcaption>**Figure 10.** Comparison between the disk masses obtained using the emitting heights from *disksurf* and those from *DISCMINER*.</figcaption></figure>

from DISCMINER. The two estimates agree across the sample, with an average ratio of 0.9 for ¹²CO and of 1.8 for ¹³CO, meaning there is no particular bias between the two methods, but the estimates can vary significantly from source to source. Indeed, the standard deviation of the mass ratio is 0.6 dex for ¹²CO and 0.4 dex for ¹³CO. The plot confirms visually that the discrepancy can be very significant for ¹²CO (particularly recognizable, as expected from the comparison shown in M. Galloway-Sprietsma et al. 2025, are HD 34282 and J1615), although when this happens the error bars are also large.

Considering that the disk mass uncertainties we report in this paper are comparable to the standard deviation we reported above, our results will not significantly change due to the choice of the way in which the emitting heights are retrieved. On the other hand, the differences that do emerge also show that there is still work to do in finding and developing the best method to retrieve emission heights.

## Appendix B
### Isothermal Approach Using ¹²CO and ¹³CO at the Same Time

We describe in the main text in Section 5.1 the results of the isothermal model when using the temperature at the emitting height. In this appendix, we describe the results when using the other approach we introduced in Section 2.3.1, namely, using two different isotopologues at the same time to jointly derive the scale-height *H* and the disk surface density.

For the example case of LkCa 15, we show the results of the isothermal model in comparison with the stratified one in the first panel of Figure 11. While in the outer disk the retrieved surface density is broadly comparable to the stratified case, we note that in the inner disk the isothermal profile is significantly denser. To illustrate why, we need to consider the second panel, where we show a comparison between the temperature returned by the isothermal model and the 2D temperature fit (M. Galloway-Sprietsma et al. 2025). Since in the latter case the temperature is a function of the vertical coordinate, we show the temperatures at three different layers: the extrapolation to the midplane and the values at the emitting heights of ¹²CO and ¹³CO. It can be seen how the temperature returned by the isothermal model becomes unphysically low in the inner disk. Eventually, for radii smaller than ≈100 au, the isothermal

profile is not defined, and no value for the surface density and temperature can be found that solves the system of equations defined in Section 2.3.1.

We attribute the behavior in the inner disk to the fact that the emitting heights of ¹²CO and ¹³CO become similar within 100 au, and within ~60 au, the emitting height of ¹³CO is higher than that of ¹²CO, as can be seen in the third panel of the figure. In our model, by construction, ¹³CO must come from a deeper layer than ¹²CO, explaining why no solution can be found inside ~60 au. At the radii where two heights are similar, Equation (16) implies that the solution must have a small *H*, so that both emitting heights are far from the center of the Gaussian, in the region where the density rapidly drops. Such a high ¹³CO emitting surface also explains why even for the vertically stratified case the surface density measured from ¹³CO is significantly higher in the inner disk than when measuring it from ¹²CO. Considering that ¹²CO is significantly more abundant than ¹³CO, this behavior is clearly unphysical, and it is likely driven by errors in the retrieval of the emitting height, in particular for ¹³CO. Because we use a parameterized surface in DISCMINER, it is possible that a more flexible parameterization would alleviate the issue.

This behavior is not unique of LkCa 15—we find no solution for the isothermal model also in the inner disk of V4046 Sgr. In addition, MWC 480 and SY Cha (in this case, in the outer disk) are also severely affected by the issue. Formally, we are still able to find solutions because the ¹³CO emitting height is smaller than the ¹²CO one, but because they become very similar, the isothermal surface density estimate becomes orders of magnitude higher than the stratified case.

Even when excluding such extreme cases, we find that using the two emitting height measurements at the same time often performs significantly less well than when using the brightness temperature at the emission surface; in some cases, it yields disk masses that can differ up to 1 order of magnitude from either of the stratified estimates. We attribute this failure to the fact that ultimately ¹²CO and ¹³CO come from layers with different temperatures. Because the model we discuss in this appendix uses their relative ratio to set the scale-height, it is more sensitive to the vertical variations of the temperature than the models using a single emission line to estimate the surface
<figure><img src="image_14.png" /><figcaption>**Figure 11.** Left panel: surface density $\Sigma$ in LkCa 15 of the stratified models in comparison with the isothermal model when solving for the scale-height using both the $^{12}$CO and $^{13}$CO emitting height at the same time. Right panel: comparison between the temperatures returned by the solution shown in the left panel and the temperature fit reported in M. Galloway-Sprietsma et al. (2025; at three different layers, the extrapolation to the midplane and the values at the emitting heights of $^{12}$CO and $^{13}$CO). Bottom panel: comparison between the emitting heights of $^{12}$CO and $^{13}$CO, showing that in the inner disk the former is smaller than the latter.</figcaption></figure>

density. Because of these reasons, we will not consider further using both emitting heights at the same time.

## Appendix C
### Density Profile Comparison between Isothermal and Stratified Models

Our method can essentially be understood as an extrapolation from the conditions at the emitting height to the midplane. The extrapolation assumes hydrostatic equilibrium and as such relies on knowledge of the temperature structure. To visualize in practice how the extrapolation happens, we show in the upper panel of Figure 12 the vertical density profile for LkCa 15 at a reference radius of 80 au, comparing the stratified model where the temperature varies with the vertical coordinate, as shown in the bottom panel, with the two isothermal models described in the text, one where we set the temperature to the value at the emitting height and one where we set the temperature to the midplane value. We also mark the $^{12}$CO emitting height with the dashed vertical line.

We normalize the profiles so that they all pass by the same point at the emitting height, with a value set by the stratified case with respect to the midplane (i.e., for the stratified case, the density in the midplane is normalized to 1). This is justified by the fact that, to the first order, the models constrain the

density to be similar (but not mathematically equal) at the emitting height, since they all require the line emission to become optically thick, and the column integrals are dominated by the region close to the emitting height. The plot shows visually that the value returned by the extrapolation for the midplane density depends on the steepness of the density profile; for example, the midplane temperature is the lowest, implying a steeper profile that returns the highest midplane density. The opposite happens when using the emitting height temperature, which therefore produces the lowest midplane density; the stratified case falls somewhere in between. We find this ordering to apply to most cases, although there are exceptions. We interpret this as due to the fact that the optical depth depends on the whole vertical integral up to infinity. The steepness of the density profile has therefore an impact in the density retrieved at the emitting height: for example, both isothermal cases fall off more steeply than the stratified case since the temperature will still increase in the upper layers, yielding for both a slightly higher density at the emitting point. This further exacerbates the overestimate for the midplane temperature case, while it partially compensates the underestimate for the emitting height temperature one, explaining why the latter generally fares better than the former. In some rare cases, we even find that the latter effect dominates and that
<figure><img src="image_15.png" /><figcaption>**Figure 12.** Top panel: density profiles in LkCa 15 at 80 au. We show the density profiles for the stratified case and for the two isothermal cases described in Section 5.1. The stratified profile is normalized to 1 in the midplane, and the other two profiles are normalized so that at the ¹²CO emitting height (marked by the dashed vertical line) they intersect with the stratified one. Bottom panel: temperature profile derived from the M. Galloway-Sprietsma et al. (2025) fit.</figcaption></figure>

the emitting height temperature model overestimates the stratified value. In practice, we note that the exact magnitude of the underestimate/overestimate depends on the exact temperature profile, and even qualitatively, there can be a large variety of cases: for example, the emitting height could be located at a height where the temperature profile is close to the midplane value, still increasing (as it is here) or where it has plateaued to the surface value. We can easily imagine that in the first case (temperature at the emitting height close to the midplane value) using the midplane temperature should give an estimate very close to the isothermal value.

The arguments above justify why in general the midplane temperature gives a worse estimate than the emitting height temperature, but we remark that the truth of the statement depends somehow on the temperature profile. Bearing in mind that both isothermal models are approximations, we consider this result as empirical, and it can only be justified a posteriori, as we have done here.

## ORCID iDs

Giovanni P. Rosotti https://orcid.org/0000-0003-4853-5736
Cristiano Longarini https://orcid.org/0000-0003-4663-0318
Teresa Paneque-Carreño https://orcid.org/0000-0002-4044-8016

Gianni Cataldi https://orcid.org/0000-0002-2700-9676

Maria Galloway-Sprietsma https://orcid.org/0000-0002-5503-5476

Sean M. Andrews https://orcid.org/0000-0003-2253-2270

Jaehan Bae https://orcid.org/0000-0001-7258-770X

Marcelo Barraza-Alfaro https://orcid.org/0000-0001-6378-7873

Myriam Benisty https://orcid.org/0000-0002-7695-7605

Pietro Curone https://orcid.org/0000-0003-2045-2154

Ian Czekala https://orcid.org/0000-0002-1483-8811

Stefano Facchini https://orcid.org/0000-0003-4689-2684

Daniele Fasano https://orcid.org/0000-0003-4679-4072

Mario Flock https://orcid.org/0000-0002-9298-3029

Misato Fukagawa https://orcid.org/0000-0003-1117-9213

Himanshi Garg https://orcid.org/0000-0002-5910-4598

Cassandra Hall https://orcid.org/0000-0002-8138-0425

Jane Huang https://orcid.org/0000-0001-6947-6072

John D. Ilee https://orcid.org/0000-0003-1008-1142

Andrés F. Izquierdo https://orcid.org/0000-0001-8446-3026

Kazuhiro Kanagawa https://orcid.org/0000-0001-7235-2417

Geoffroy Lesur https://orcid.org/0000-0002-8896-9435

Giuseppe Lodato https://orcid.org/0000-0002-2357-7692

Ryan A. Loomis https://orcid.org/0000-0002-8932-1219

Ryuta Orihara https://orcid.org/0000-0003-4039-8933

Christophe Pinte https://orcid.org/0000-0001-5907-5179

Daniel J. Price https://orcid.org/0000-0002-4716-4235

Jochen Stadler https://orcid.org/0000-0002-0491-143X

Richard Teague https://orcid.org/0000-0003-1534-5186

Gaylor Wafflard-Fernandez https://orcid.org/0000-0002-3468-9577

Andrew J. Winter https://orcid.org/0000-0002-7501-9801

Lisa Wölfer https://orcid.org/0000-0002-7212-2416

Hsi-Wei Yen https://orcid.org/0000-0003-1412-893X

Tomohiro C. Yoshida https://orcid.org/0000-0001-8002-8473

Brianna Zawadzki https://orcid.org/0000-0001-9319-1296

## References

Aikawa, Y., van Zadelhoff, G. J., van Dishoeck, E. F., & Herbst, E. 2002, A&A, 386, 622

Anderson, D. E., Blake, G. A., Bergin, E. A., et al. 2019, ApJ, 881, 127

Anderson, D. E., Cleeves, L. I., Blake, G. A., et al. 2022, ApJ, 927, 229

Andrews, S. M., Teague, R., Wirth, C. P., Huang, J., & Zhu, Z. 2024, ApJ, 970, 153

Barbier, D. 1943, AnAp, 6, 113

Bergin, E. A., Bosman, A., Teague, R., et al. 2024, ApJ, 965, 147

Bergin, E. A., & Williams, J. P. 2017, in Formation, Evolution, and Dynamics of Young Solar Systems, Astrophysics and Space Science Library, ed. M. Pessah & O. Gressel, Vol. 445 (Cham: Springer), 1

Bruderer, S. 2013, A&A, 559, A46

Casassus, S., & Pérez, S. 2019, ApJL, 883, L41

Condon, J. J., & Ransom, S. M. 2016, Essential Radio Astronomy (Princeton, NJ: Princeton Univ. Press)

Flaherty, K. M., Hughes, A. M., Rosenfeld, K. A., et al. 2015, ApJ, 813, 99

Flaherty, K., Hughes, A. M., Simon, J. B., et al. 2020, ApJ, 895, 109

Flaherty, K., Hughes, A. M., Simon, J. B., et al. 2024, MNRAS, 532, 363

Flock, M., Fromang, S., González, M., & Commerçon, B. 2013, A&A, 560, A43

Galloway-Sprietsma, M., Bae, J., Izquierdo, A., et al. 2025, ApJL, 984, L10

Izquierdo, A., Stadler, J., Galloway-Sprietsma, M., et al. 2025, ApJL, 984, L8

Izquierdo, A. F., Testi, L., Facchini, S., et al. 2023, A&A, 674, A113

Izquierdo, A. F., Testi, L., Facchini, S., Rosotti, G. P., & van Dishoeck, E. F. 2021, A&A, 650, A179

Kelly, B. C. 2007, ApJ, 665, 1489

Law, C. J., Benisty, M., Facchini, S., et al. 2024, ApJ, 964, 190

Law, C. J., Crystian, S., Teague, R., et al. 2022, ApJ, 932, 114

Law, C. J., Teague, R., Loomis, R. A., et al. 2021, ApJS, 257, 4

Law, C. J., Teague, R., Öberg, K. I., et al. 2023, ApJ, 948, 60

Longarini, C., Lodato, G., Rosotti, G., et al. 2025, ApJL, 984, L17

Lynden-Bell, D., & Pringle, J. E. 1974, MNRAS, 168, 603

Martire, P., Longarini, C., Lodato, G., et al. 2024, A&A, 686, A9

McClure, M. K., Bergin, E. A., Cleeves, L. I., et al. 2016, ApJ, 831, 167

Miotello, A., Bruderer, S., & van Dishoeck, E. F. 2014, A&A, 572, A96

Miotello, A., Kamp, I., Birnstiel, T., Cleeves, L. C., & Kataoka, A. 2023, in ASP Conf. Ser. 534, Protostars and Planets VII, ed. S. Inutsuka et al. (San Francisco, CA: ASP), 501

Miotello, A., van Dishoeck, E. F., Williams, J. P., et al. 2017, A&A, 599, A113

Öberg, K. I., Guzmán, V. V., Walsh, C., et al. 2021, ApJS, 257, 1

Paletou, F. 2018, OAst, 27, 76

Paneque-Carreño, T., Izquierdo, A. F., Teague, R., et al. 2024, A&A, 684, A174

Paneque-Carreño, T., Miotello, A., van Dishoeck, E. F., et al. 2023, A&A, 669, A126
Paneque-Carreño, T., Miotello, A., van Dishoeck, E. F., Rosotti, G., & Tabone, B. 2025, arXiv:2501.08294

Pezzotta, V., Facchini, S., Longarini, C., Lodato, G., & Martire, P. 2025, A&A, 694, A108

Pinte, C., Ménard, F., Duchêne, G., et al. 2018, A&A, 609, A47

Pinte, C., Teague, R., Flaherty, K., et al. 2023, in ASP Conf. Ser. 534, Protostars and Planets VII, ed. S. Inutsuka et al. (San Francisco, CA: ASP), 645

Qi, C., D'Alessio, P., Öberg, K. I., et al. 2011, ApJ, 740, 84

Rosenfeld, K. A., Andrews, S. M., Hughes, A. M., Wilner, D. J., & Qi, C. 2013, ApJ, 774, 16

Rosotti, G. P. 2023, NewAR, 96, 101674

Rosotti, G. P., Benisty, M., Juhász, A., et al. 2020b, MNRAS, 491, 1335

Rosotti, G. P., Teague, R., Dullemond, C., Booth, R. A., & Clarke, C. J. 2020a, MNRAS, 495, 173

Schöier, F. L., van der Tak, F. F. S., van Dishoeck, E. F., & Black, J. H. 2005, A&A, 432, 369

Stapper, L. M., Hogerheijde, M. R., van Dishoeck, E. F., & Paneque-Carreño, T. 2023, A&A, 669, A158

Sturm, J. A., Booth, A. S., McClure, M. K., Leemker, M., & van Dishoeck, E. F. 2023, A&A, 670, A12

Takeuchi, T., & Lin, D. N. C. 2002, ApJ, 581, 1344

Teague, R. 2019, JOSS, 4, 1220

Teague, R., Benisty, M., Facchini, S., et al. 2025, ApJL, 984, L6

Toci, C., Rosotti, G., Lodato, G., Testi, L., & Trapman, L. 2021, MNRAS, 507, 818

Trapman, L., Longarini, C., Rosotti, G., et al. 2025, ApJL, 984, L18

Trapman, L., Miotello, A., Kama, M., van Dishoeck, E. F., & Bruderer, S. 2017, A&A, 605, A69

Trapman, L., Rosotti, G., Zhang, K., & Tabone, B. 2023, ApJ, 954, 41

Trapman, L., Zhang, K., van't Hoff, M. L. R., Hogerheijde, M. R., & Bergin, E. A. 2022, ApJL, 926, L2

Urbina, F., Miley, J., Kama, M., & Keyte, L. 2024, A&A, 686, A120

van der Tak, F. F. S., Black, J. H., Schöier, F. L., Jansen, D. J., & van Dishoeck, E. F. 2007, A&A, 468, 627

van Dishoeck, E. F., & Black, J. H. 1988, ApJ, 334, 771

Visser, R., van Dishoeck, E. F., & Black, J. H. 2009, A&A, 503, 323

Williams, J. P., & Best, W. M. J. 2014, ApJ, 788, 59

Wilson, T. L., & Rood, R. 1994, ARA&A, 32, 191

Woitke, P., Min, M., Pinte, C., et al. 2016, A&A, 586, A103

Yoshida, C. T., Curone, P., Stadler, J., et al. 2025, ApJL, 984, L19

Yoshida, T. C., Nomura, H., Furuya, K., Tsukagoshi, T., & Lee, S. 2022, ApJ, 932, 126

Zhang, K., Booth, A. S., Law, C. J., et al. 2021, ApJS, 257, 5