exoALMA. V. Gaseous Emission Surfaces and
Temperature Structures

Maria Galloway-Sprietsma, Jaehan Bae, Andrés F. Izquierdo, Jochen Stadler,
Cristiano Longarini, Richard Teague, Sean M. Andrews, Andrew J. Winter,
Myriam Benisty, Stefano Facchini, et al.

► To cite this version:

Maria Galloway-Sprietsma, Jaehan Bae, Andrés F. Izquierdo, Jochen Stadler, Cristiano Longarini, et al.. exoALMA. V. Gaseous Emission Surfaces and Temperature Structures. The Astrophysical Journal Letters, 2025, 984, 10.3847/2041-8213/adc437 . insu-05099881

HAL Id: insu-05099881

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

Submitted on 5 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
exoALMA. V. Gaseous Emission Surfaces and Temperature Structures

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

**Abstract**

An analysis of the gaseous component in protoplanetary disks can inform us about their thermal and physical structure, chemical composition, and kinematic properties, all of which are crucial for understanding various processes within the disks. By exploiting the asymmetry of the line emission, or via line profile analysis, we can locate the emitting surfaces. Here, we present the emission surfaces of the exoALMA sources in $^{12}$CO J = 3–2, $^{13}$CO J = 3–2, and CS J = 7–6. We find that $^{12}$CO traces the upper disk atmosphere, with mean $\langle z/r \rangle$ values of ≈0.28, while $^{13}$CO and CS trace lower regions of the disk with mean $\langle z/r \rangle$ values of ≈0.16 and ≈0.18, respectively. We find that $^{12}$CO $\langle z/r \rangle$ and the disk mass are positively correlated with each other; this relationship offers a straightforward way to infer the disk mass. We derive 2D $r-z$ temperature distributions of the disks. Additionally, we search for substructure in the surfaces and radial intensity profiles; we find evidence of localized substructure in the emission surfaces and peak intensity profiles of nearly every disk, with this substructure often being coincident between molecular tracers, intensity profiles, and kinematic perturbations. Four disks display evidence of potential photodesorption, implying that this effect may be common even in low far-ultraviolet star-forming regions. For most disks, we find that the physical and thermal structure is more complex than analytical models can account for, highlighting a need for more theoretical work and a better understanding of the role of projection effects on our observations.

*Unified Astronomy Thesaurus concepts: Protoplanetary disks (1300); Planet formation (1241); CO line emission (262); High angular resolution (2167)*

**1. Introduction**

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

Planets are born within the gaseous and dusty environments of circumstellar disks. Despite the many observations of fully formed planets, our knowledge of the preceding planet formation process is still largely incomplete. In order to best understand this

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

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

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

<footnote n="4">⁴Université Côte d'Azur, Observatoire de la Côte d'Azur, CNRS, Laboratoire Lagrange, France</footnote>

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

<footnote n="6">⁶Dipartimento di Fisica, Università degli Studi di Milano, Via Celoria 16, 20133 Milano, Italy</footnote>

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

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

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

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

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

<footnote n="13">¹³School of Physics and Astronomy, Monash University, VIC 3800, Australia</footnote>

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

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

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

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

<footnote n="18">¹⁸Department of Physics and Astronomy, Rice University, 6100 Main St., Houston, TX 77005, USA</footnote>

<footnote n="19">¹⁹Los Alamos National Laboratory, Los Alamos, NM 87545, USA</footnote>

<footnote n="20">²⁰School of Physics and Astronomy, Monash University, Clayton VIC 3800, Australia</footnote>

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

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

<footnote n="23">²³Institute for Artificial Intelligence, The University of Georgia, Athens, GA 30602, USA</footnote>

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

<footnote n="25">²⁵School of Physics and Astronomy, University of Leeds, Leeds, LS2 9JT, UK</footnote>

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

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

<footnote n="28">²⁸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="29">²⁹Department of Astronomical Science, The Graduate University for Advanced Studies, SOKENDAI, 2-21-1 Osawa, Mitaka, Tokyo 181-8588, Japan
Received 2024 November 25; revised 2025 March 14; accepted 2025 March 15; published 2025 April 28</footnote>

<footnote n="30">³⁰ 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>
process, we must first characterize the physical and chemical
properties of protoplanetary disks. In recent years, the high
angular and spectral resolution offered by the Atacama Large
Millimeter/submillimeter Array (ALMA) has allowed us to map
planetary nurseries in deeper detail than ever before. In
particular, molecular line observations of protoplanetary disks
can inform us about the gas structure, which holds a wealth of
information about the disk and potential planets forming within
(A. Miotello et al. 2023).

Protoplanetary disks are composed of two main components: dust and gas. The millimeter and larger dust grains generally settle toward the midplane due to gravity exerted by the central star, aerodynamical gas drag, and a vertical pressure gradient (L. Barrière-Fouchet et al. 2005). In contrast, the bulk of the gaseous disk component extends beyond the dust both radially and vertically (e.g., M. Ansdell et al. 2018; C. J. Law et al. 2021a; F. Long et al. 2022a). Although the gas is primarily composed of radio-quiet molecular hydrogen (H₂), there are trace amounts of other molecules that we can readily observe with ALMA. Because these molecules have differing opacities, observations targeting various molecular tracers allow us to probe the 2D structure of the disk, from the upper atmospheres down to the midplane. The vertical and radial distribution of these molecules is dependent on the physical conditions found within the disk, particularly the temperature, density, and radiation (G. J. van Zadelhoff et al. 2001; C. Walsh et al. 2015; K. I. Öberg et al. 2023). We can directly measure the vertical emission surfaces for the frontside and backside of the disk via several methods. The vertical extent of the surfaces can be directly probed in edge-on disks, as in A. Dutrey et al. (2017), M. Villenave et al. (2020), G. Duchêne et al. (2024), and M. Villenave et al. (2024). For midinclination disks (generally 30°–60°), the geometry of the disk can be exploited to determine isovelocity contours in each channel. The frontside and backside emission will be well separated for disks of these intermediate inclinations, thus allowing us access to the full 3D structure covering the entire azimuth; see E. Dartois et al. (2003), C. Pinte et al. (2018), R. Teague et al. (2021a), and T. Paneque-Carreño et al. (2023). A line-profile analysis also allows for a measurement of the emission surfaces by identifying the multiple components of the line from the frontside and backside of the disk (A. F. Izquierdo et al. 2021, 2023). These techniques can identify CO snow lines (C. Pinte et al. 2018), map the 2D structure of disks (C. J. Law et al. 2021a; E. A. Rich et al. 2021; C. J. Law et al. 2022; T. Paneque-Carreño et al. 2023; L. M. Stapper et al. 2023; C. Hernández-Vera et al. 2024), identify substructures in emission surfaces that relate to disk physical processes (R. Teague et al. 2021b; T. Paneque-Carreño et al. 2022; M. Galloway-Sprietsma et al. 2023; A. F. Izquierdo et al. 2023; F. Urbina et al. 2024; T. Paneque-Carreño et al. 2025), and infer turbulence levels at multiple vertical locations (T. Paneque-Carreño et al. 2024). Importantly, knowing the radial and vertical location of molecular line emission is critical to interpreting the kinematic signatures in the gas (see review by C. Pinte et al. 2023).

The ALMA Large Program exoALMA aims to detect plants
forming in situ at submillimeter wavelengths, utilizing
kinematic planet-detection techniques. To achieve those goals,
15 protoplanetary disks were observed at high angular and
spectral resolution to target the molecular line emission of ¹²CO
J = 3–2, ¹³CO J = 3–2, and CS J = 7–6. In this Letter, we

present the emission surfaces, intensity radial profiles, and 2D
temperature profiles of the 15 sources. In Section 2, we
describe the calibration and imaging of the data utilized in this
Letter. In Section 3, we outline the methods used to derive the
emission surfaces and the results. In Section 4, we detail
the methods used to obtain radial temperature profiles and the
subsequent 2D temperature profile calculation. In Section 5, we
identify substructures in the surfaces and radial profiles, and
compare these to continuum structures and kinematic perturba-
tions. We discuss the results and implications of this study in
Section 6, and summarize our work in Section 7.

## 2. Data

The data used for this analysis were obtained for exoALMA,<fnref n="31" />
an ALMA Large Program (2021.1.01123.L; R. Teague et al.
2025). This program observed 15 protoplanetary disks, most
with previously observed continuum substructures indicative of
ongoing planet formation processes. Observations of ¹²CO
J = 3–2, ¹³CO J = 3–2, and CS J = 7–6 at 0.′1 angular
resolution and a 35.5 kHz (26 m s⁻¹) spectral resolution with
a sensitivity of 3 K over 150 m s⁻¹ were requested to ensure a
detection of line emission in the outer disk where we expect the
gas temperature, T<sub>B</sub>, to be ≈20 K for all molecules. A complete
description of observations, source selection, and imaging
techniques used are described in R. Teague et al. (2025) and
R. A. Loomis et al. (2025). Figure 1 presents the exoALMA
disk sample in ¹²CO J = 3–2 moment eight maps, depicting the
peak intensity. These moment maps were made using
bettermoments (R. Teague & D. Foreman-Mackey 2018),
with the quadratic method and assuming the full Planck
function. The fiducial images are shown in the blue contour
(0.′15, 100 m s⁻¹ resolution) and the high surface brightness
sensitivity images are in the background in gray (0.′30,
100 m s⁻¹ resolution). A 5σ signal-to-noise ratio (SNR) clip
has been applied.

## 3. Emission Surfaces

In this section, we discuss the method used to constrain the
emission surfaces of the 15 target disks and then present the
results. The nonparametric surfaces discussed represent the radial
and vertical points, and the parametric surfaces represent the
tapered power-law fits. We report surfaces for the three target
molecules: ¹²CO J = 3–2, ¹³CO J = 3–2, and CS J = 7–6.

3.1. *Nonparametric Surfaces*

To derive the nonparametric surfaces, we employ
disksurf<fnref n="32" /> (R. Teague et al. 2021a), which utilizes the
method outlined in C. Pinte et al. (2018). This method relies on
the asymmetry of the line emission above the disk midplane to
infer an emission height, which allows us to locate emission
arising from specific locations in the disk. Using this method, we
obtain the deprojected radius *r* and emission height *z*, along with
the brightness temperature *T*<sub>B</sub> and the SNR for each point, which
will be used to derive 2D temperature structure later in Section 4.

Unless otherwise specified, all of the nonparametric surfaces
were derived using the non-continuum-subtracted cubes with
100 m s⁻¹ resolution and 0.′15 beam size for ¹²CO and ¹³CO,
and the continuum-subtracted cubes with 200 m s⁻¹ resolution

<footnote n="31">³¹ https://www.exoalma.com</footnote>

<footnote n="32">³² https://github.com/richteague/disksurf</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Peak intensity maps of ¹²CO J = 3–2 for all exoALMA sources made with the fiducial images. Parametric tapered power-law surfaces from DISCMINER are overlaid in white, with each contour being separated by the beam size (0.′′15); the bold lines represent the major and minor axes. The synthesized beam is shown in the lower left corner of each plot. The 0.′′30 beam cubes are plotted in the background in gray.</figcaption></figure>

and 0.′′15 beam size for CS. One exception to this is DM Tau, where we use the 0.′′30 beam size cubes for both ¹²CO and ¹³CO surface extraction to retrieve the emission surface out to a larger radial extent due to the low SNR in the 0.′′15 cubes (see DM Tau in Figure 1).

We opt for the non-continuum-subtracted cubes for ¹²CO and ¹³CO to ensure we retrieve accurate temperature measurements (discussed more in Section 4). We found no significant differences between the surface structure retrieved from continuum-subtracted versus non-continuum-subtracted cubes for ¹²CO and ¹³CO. However, the CS emission tends to have comparable intensities to the continuum, affecting the surface retrievals. In some cases, the continuum emission was confused for the gaseous emission; therefore, we opt to use the continuum-subtracted cubes for the CS surfaces. We use geometric parameters ($M_*$, $i$, PA, $v_{LSR}$, $x_0$, $y_0$) derived from DISCMINER (see R. Teague et al. 2025; A. F. Izquierdo et al. 2025).

We employ several methods to improve emission surfaces obtained with *disksurf* by using its updated functionalities. First, we use the iterative surface-finding approach, which relies on a previously derived surface and rotation curve, $z(r)$ and $v_\phi(r)$. This surface and rotation curve is then used to draw a mask for the subsequent iteration, using a specified beam size and SNR. The SNR for each point is measured as the intensity divided by the cube rms. We repeat the iteration five times for each surface, allowing us to decrease the maximum beam size from two beams for the first iteration to one beam on the subsequent iterations, and increase minimum SNR linearly

from 1.0 to 5.0 for ¹²CO. For ¹³CO, we increased the SNR from 1.0 to 3.5, and for CS, we increased the SNR from 1.0 to 3.0. Therefore, the final SNR clips for the surfaces are 5.0 for ¹²CO, 3.5 for ¹³CO, and 3.0 for CS. This iterative approach does not allow for surface points to fall below $z = 0$; therefore, we do not use it for the back surfaces, and all front surfaces are masked below this point. Despite this, when comparing the surfaces retrieved from the noniterative versus iterative, there are not significant differences in the height of the surfaces. We apply additional masks to cut out visually identified noise, as well as points that belong to more diffuse CO emission in the outer disk, which does not clearly belong to either the front or back surfaces (see Section 6.3). These steps were done on a per-disk basis, and were visually justified by looking at the channel maps with the detected points overlaid. The individual $r-z$ points derived via this process comprise the nonparametric surfaces. These surfaces are presented in Figure 2 as both raw and binned data points; all points are for the front surfaces only. The points are binned by a quarter of a beam size. For individual molecular surfaces, see Figures 9, 10, and 11 in Appendix A. We have excluded J1604 and HD 135344B from this analysis, as they are too face-on to reliably separate the front and back surfaces using the *disksurf* method.

## 3.2. Parametric Fitting

We find a parametric surface using both *disksurf* and DISCMINER<fnref n="33" /> (see A. F. Izquierdo et al. 2021, 2023).

<footnote n="33">³³ https://github.com/andizq/discminer</footnote>
<figure><img src="image_4.png" /><figcaption>Figure 2. Raw and binned front surfaces for disks with viable emission surfaces in $^{12}$CO J = 3–2 (blue), $^{13}$CO J = 3–2 (orange), and CS J = 7–6 (purple). The binned diffuse backside points of AA Tau, SY Cha, and LkCa 15 are shown as “x” points. Lines of constant $\langle z/r \rangle$ from 0.1 to 0.6 are shown as gray background lines. The beam size is shown in the upper right corner of each plot; the gray shaded region corresponds to everything interior to one beam size.</figcaption></figure>

DISCMINER fits a Keplerian model to each channel, and has been used to derive the geometric parameters for each exoALMA source (A. F. Izquierdo et al. 2025). We assume a tapered power-law analytic prescription as the surface model, given by

$$z(r) = z_0 \left( \frac{r}{100 \text{ au}} \right)^\psi \exp \left( - \left[ \frac{r}{r_{\text{taper}}} \right]^{q_{\text{taper}}} \right), \quad (1)$$

where $z_0$ is the normalization factor of the surface at 100 au, $\psi$ is the power-law exponent, $r_{\text{taper}}$ is the characteristic radius for the exponential taper, and $q_{\text{taper}}$ is the exponent of the taper term.

This tapered power-law model is fit to the nonparametric surfaces using the Markov Chain Monte Carlo method from disksurf, which uses emcee (D. Foreman-Mackey et al. 2013), a Python implementation of an affine invariant ensemble sampler. We employ 256 walkers, 100 burn-in steps, and 1000 steps to obtain a fit to the nonparametric surface data points. The best-fit surfaces represent the 50th percentile parameters

from the posterior distributions, and the errors are obtained from the 16th and 84th percentiles.

We present $z_0$, $r_{\text{taper}}$, $\psi$, and $q_{\text{taper}}$, derived from disksurf and DISCMINER, in Table 1. Figure 1 presents the DISCMINER tapered power-law fits overlaid on peak intensity maps of each disk to demonstrate the 3D geometry of the surfaces. Figures 9, 10, and 11 in Appendix A compare the parametric surfaces from disksurf and DISCMINER, overlaid on the nonparametric surfaces. The differences between the parametric surfaces are due to the underlying methodological differences between disksurf and DISCMINER; the disksurf parametric surfaces are fit to the data points extracted, which are shown in the plot, while DISCMINER fits a model to the line profiles. In some cases, the difference may be due to backside contamination; DISCMINER can readily distinguish between emission stemming from the front versus back surface via line profile modeling, but disksurf relies on finding peaks along isovelocity contours. Additionally, if there is a large deviation from Keplerian rotation, DISCMINER takes this into account during
<table><caption>**Table 1**
Exponentially Tapered Power-law Fits for All Disks and All Molecules from *disksurf* and DISCMINER; Characteristic $\langle z/r \rangle$ Values Determined from *disksurf* Data Points</caption><thead><tr><th rowspan="2">Source</th><th rowspan="2">Line</th><th rowspan="2">z / r</th><th colspan="4">disksurf</th><th colspan="4">DISCMINER</th></tr><tr><th>z<sub>0</sub><br>(au)</th><th>ψ</th><th>r<sub>taper</sub><br>(au)</th><th>q<sub>taper</sub></th><th>z<sub>0</sub><br>(au)</th><th>ψ</th><th>r<sub>taper</sub><br>(au)</th><th>q<sub>taper</sub></th></tr></thead><tbody><tr><td rowspan="3">DM Tau</td><td><sup>12</sup>CO J = 3–2</td><td>0.34</td><td>25.2</td><td>1.8</td><td>366</td><td>3.2</td><td>87</td><td>1.9</td><td>79.6</td><td>0.48</td></tr><tr><td><sup>13</sup>CO J = 3–2</td><td>0.14</td><td>21.5</td><td>1.9</td><td>262</td><td>0.35</td><td>20</td><td>2.3</td><td>242</td><td>0.93</td></tr><tr><td>CS J = 7–6</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>6.5</td><td>3.2</td><td>94.1</td><td>0.92</td></tr><tr><td rowspan="3">AA Tau</td><td><sup>12</sup>CO J = 3–2</td><td>0.41</td><td>42.5</td><td>0.79</td><td>289</td><td>4.8</td><td>50</td><td>1.2</td><td>240</td><td>1.3</td></tr><tr><td><sup>13</sup>CO J = 3–2</td><td>0.32</td><td>33.4</td><td>0.92</td><td>309</td><td>2.6</td><td>52</td><td>1.4</td><td>151</td><td>1.4</td></tr><tr><td>CS J = 7–6</td><td>0.32</td><td>70.9</td><td>1</td><td>164</td><td>0.28</td><td>33</td><td>0.83</td><td>232</td><td>2.2</td></tr><tr><td rowspan="3">LkCa 15</td><td><sup>12</sup>CO J = 3–2</td><td>0.34</td><td>32.2</td><td>1.1</td><td>649</td><td>4.8</td><td>29</td><td>1.1</td><td>795</td><td>3.2</td></tr><tr><td><sup>13</sup>CO J = 3–2</td><td>0.21</td><td>36.4</td><td>1.2</td><td>399</td><td>0.65</td><td>27</td><td>0.87</td><td>511</td><td>3.5</td></tr><tr><td>CS J = 7–6</td><td>0.21</td><td>28.3</td><td>0.94</td><td>244</td><td>1.5</td><td>29</td><td>0.72</td><td>303</td><td>4.6</td></tr><tr><td rowspan="3">HD 34282</td><td><sup>12</sup>CO J = 3–2</td><td>0.41</td><td>41.3</td><td>1</td><td>480</td><td>4.9</td><td>34</td><td>1.2</td><td>512</td><td>3.2</td></tr><tr><td><sup>13</sup>CO J = 3–2</td><td>0.15</td><td>12.8</td><td>2</td><td>347</td><td>4.6</td><td>27</td><td>0.79</td><td>510</td><td>4.4</td></tr><tr><td>CS J = 7–6</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>16</td><td>1.7</td><td>268</td><td>1.8</td></tr><tr><td rowspan="3">MWC 758</td><td><sup>12</sup>CO J = 3–2</td><td>0.29</td><td>64.9</td><td>0.96</td><td>176</td><td>0.24</td><td>16</td><td>0.97</td><td>254</td><td>5.3</td></tr><tr><td><sup>13</sup>CO J = 3–2</td><td>...</td><td>34.1</td><td>0.89</td><td>204</td><td>1.1</td><td>67</td><td>3.2</td><td>12.5</td><td>0.71</td></tr><tr><td>CS J = 7–6</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>6.1</td><td>5</td><td>120</td><td>3.3</td></tr><tr><td rowspan="3">CQ Tau</td><td><sup>12</sup>CO J = 3–2</td><td>0.17</td><td>21.7</td><td>1.1</td><td>187</td><td>1.9</td><td>42</td><td>1.2</td><td>346</td><td>0.087</td></tr><tr><td><sup>13</sup>CO J = 3–2</td><td>0.14</td><td>53.8</td><td>1.8</td><td>73.8</td><td>2</td><td>39</td><td>1.1</td><td>16.6</td><td>0.41</td></tr><tr><td>CS J = 7–6</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>28</td><td>4.6</td><td>80.2</td><td>1.5</td></tr><tr><td rowspan="3">SY Cha</td><td><sup>12</sup>CO J = 3–2</td><td>0.29</td><td>40.1</td><td>1.6</td><td>335</td><td>1</td><td>43</td><td>1.8</td><td>210</td><td>1</td></tr><tr><td><sup>13</sup>CO J = 3–2</td><td>0.16</td><td>21.3</td><td>2.3</td><td>185</td><td>2.4</td><td>73</td><td>2.4</td><td>66.1</td><td>0.7</td></tr><tr><td>CS J = 7–6</td><td>0.16</td><td>34.1</td><td>3.1</td><td>128</td><td>3.2</td><td>50</td><td>1.9</td><td>124</td><td>1.7</td></tr><tr><td rowspan="3">PDS 66</td><td><sup>12</sup>CO J = 3–2</td><td>0.12</td><td>12.4</td><td>0.81</td><td>143</td><td>1.4</td><td>17</td><td>1.8</td><td>127</td><td>4.5</td></tr><tr><td><sup>13</sup>CO J = 3–2</td><td>...</td><td>13.5</td><td>0.14</td><td>43.2</td><td>0.79</td><td>7.5</td><td>1.2</td><td>29</td><td>1.5</td></tr><tr><td>CS J = 7–6</td><td>...</td><td>38.1</td><td>2.5</td><td>62.6</td><td>1.7</td><td>1.2</td><td>2.9</td><td>91.7</td><td>8.4</td></tr><tr><td rowspan="3">HD 135344B</td><td><sup>12</sup>CO J = 3–2</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>14</td><td>1.4</td><td>226</td><td>10</td></tr><tr><td><sup>13</sup>CO J = 3–2</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>CS J = 7–6</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td rowspan="3">HD 143006</td><td><sup>12</sup>CO J = 3–2</td><td>0.4</td><td>49.6</td><td>1.1</td><td>146</td><td>3.4</td><td>40</td><td>1.9</td><td>161</td><td>6</td></tr><tr><td><sup>13</sup>CO J = 3–2</td><td>0.29</td><td>43.4</td><td>0.91</td><td>147</td><td>1.8</td><td>24</td><td>2.2</td><td>146</td><td>13</td></tr><tr><td>CS J = 7–6</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>17</td><td>0.72</td><td>103</td><td>4.5</td></tr><tr><td rowspan="3">RXJ1604.3-2130 A</td><td><sup>12</sup>CO J = 3–2</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td><sup>13</sup>CO J = 3–2</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>CS J = 7–6</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td rowspan="3">RXJ1615.3-3255</td><td><sup>12</sup>CO J = 3–2</td><td>0.34</td><td>27.9</td><td>1.5</td><td>403</td><td>5</td><td>26</td><td>1</td><td>530</td><td>6.9</td></tr><tr><td><sup>13</sup>CO J = 3–2</td><td>0.17</td><td>17.7</td><td>1.2</td><td>381</td><td>4.9</td><td>19</td><td>1</td><td>425</td><td>5.9</td></tr><tr><td>CS J = 7–6</td><td>0.17</td><td>14.9</td><td>2.3</td><td>185</td><td>1.3</td><td>38</td><td>2.2</td><td>121</td><td>0.8</td></tr><tr><td rowspan="3">V4046 Sgr</td><td><sup>12</sup>CO J = 3–2</td><td>0.13</td><td>26.4</td><td>1.5</td><td>196</td><td>0.94</td><td>26</td><td>1.8</td><td>151</td><td>1.2</td></tr><tr><td><sup>13</sup>CO J = 3–2</td><td>0.052</td><td>12.9</td><td>0.54</td><td>257</td><td>0.093</td><td>33</td><td>1.6</td><td>65.6</td><td>1.1</td></tr><tr><td>CS J = 7–6</td><td>0.052</td><td>5.52</td><td>1</td><td>153</td><td>0.86</td><td>0.086</td><td>0.0039</td><td>85.9</td><td>7.1</td></tr><tr><td rowspan="3">RXJ1842.9-3532</td><td><sup>12</sup>CO J = 3–2</td><td>0.2</td><td>47.4</td><td>2.2</td><td>124</td><td>1.5</td><td>26</td><td>1.5</td><td>211</td><td>1.9</td></tr><tr><td><sup>13</sup>CO J = 3–2</td><td>0.091</td><td>15</td><td>0.61</td><td>207</td><td>0.98</td><td>17</td><td>1.7</td><td>143</td><td>2</td></tr><tr><td>CS J = 7–6</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>29</td><td>1.4</td><td>145</td><td>4.5</td></tr><tr><td rowspan="3">RXJ1852.3-3700</td><td><sup>12</sup>CO J = 3–2</td><td>0.15</td><td>63.1</td><td>2.2</td><td>82.1</td><td>1.1</td><td>75</td><td>1.8</td><td>60.9</td><td>0.84</td></tr><tr><td><sup>13</sup>CO J = 3–2</td><td>0.079</td><td>15.7</td><td>1.4</td><td>155</td><td>1.3</td><td>31</td><td>2.7</td><td>90.3</td><td>1.3</td></tr><tr><td>CS J = 7–6</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>1</td><td>3.6</td><td>108</td><td>1.4</td></tr></tbody></table>

modeling, but *disksurf* assumes a purely circular motion. However, the definite underlying cause of the differences is still unknown; a comparison of the two codes is discussed in Appendix D. Despite the slight parametric differences, both surfaces are broadly consistent.

### 3.3. Emission Surface Results

We obtain nonparametric emission surfaces in $^{12}$CO and $^{13}$CO J = 3–2 for 13 out of the 15 exoALMA sources, and six nonparametric emission surfaces in CS J = 7–6. Although
emission was detected for each molecule from each disk, we
are only able to characterize a surface for those with high
enough SNR. The two disk surfaces we do not report belong to
HD 135344B and J1604, both of which are too face-on to
retrieve accurate emission heights using the *disksurf* and
DISCMINER methods. The raw and binned surfaces for each
molecule are presented in Figure 2. The surfaces plotted
individually by line are shown in Figures 9, 10, and 11 in
Appendix A. Figure 12, also in Appendix A, plots all rolling
surfaces together for comparison, where the rolling surface
represents a moving average of the binned data points. We find
that, in most disks, the ¹²CO surface lies above that of the ¹³CO
surface, and the ¹³CO surface generally lies above that of CS
surface, which is expected as reaching an optical depth of ~1
will happen at deeper layers for less abundant molecules, and
different molecules will be distributed differently throughout
the disk. However, this is not true of all disks; in particular, the
¹³CO and CS molecular surfaces of AA Tau, J1615, LkCa 15,
SY Cha, and V4046 Sgr overlap and appear to be stemming
from similar heights, at least in the inner disk regions (e.g.,
~150 au). In nearly all disks, the ¹²CO surfaces extend radially
farther out than the ¹³CO and CS surfaces, which is due to the
higher SNR achieved with ¹²CO in those regions, as well as the
intrinsic distribution of molecules due to the temperature and
density. DM Tau hosts the largest ¹²CO surface in both height
and radius, extending out to ~850 au and reaching a height of
~250 au, although we note most of the emission in the outer
disk is low SNR and scattered. PDS 66 has the lowest surface
among the molecules detected; it is consistent with a mostly flat
surface (e.g., consistent with an origin in the disk midplane) in
all molecules and extends out to a radius of ~150 au. While
nearly all surfaces follow the expected power-law or tapered
power-law shape, V4046 Sgr is a notable exception with a
double-peaked emission surface in ¹²CO. This will be
discussed later in Section 6.

We calculate a representative $\langle z/r \rangle$ for each surface by taking an average of the binned surface points. We cut points past the peak of the emission surface, thus excluding those that belong to the tapered portion of the surfaces. This differs slightly from the method employed in previous works such as C. J. Law et al. (2021a) or T. Paneque-Carreño et al. (2023) who use $0.8 \times r_{\text{taper}}$ to make the radial cut. We opt to use the radial location of the peak of the surface due to a few surface fits having $r_{\text{taper}}$ values interior to most of the rising emission; these two methods do not give discrepant results, and the method employed here still captures only the rising portion of the emission surfaces, which we verified by eye. We additionally exclude any data interior to one beam size in the $\langle z/r \rangle$ calculation. We have verified that the $\langle z/r \rangle$ value is consistent when measured at the same radial location in each disk (125 au), and is not inherently influenced by the radial extent of the disk.

We find $\langle z/r \rangle$ values for ¹²CO emission ranging from ~0.13 to ~0.4, with an overall average $\langle z/r \rangle$ of 0.28. AA Tau hosts the most elevated ¹²CO surface with a $\langle z/r \rangle$ of 0.41, whereas PDS 66 hosts the flattest one, with a $\langle z/r \rangle$ of 0.12. The ¹³CO surfaces range in $\langle z/r \rangle$ from 0.05 (excluding PDS 66, which exhibits a completely flat surface in ¹³CO) to 0.32, with an average of 0.16. The ¹³CO surfaces generally do not exhibit a taper in the outer disk, likely because the SNR is too low in this portion of the disks. We retrieve CS emission surfaces for six disks. The remaining disks had too low SNR to retrieve a

viable surface. Of these six, two are consistent with a flat
surface close to z = 0 (PDS 66 and V4046 Sgr). The remaining
CS surfaces are elevated, and that of SY Cha even exhibits the
characteristic tapered emission surface. The average $\langle z/r \rangle$ for
CS is 0.18, with AA Tau having the largest at ~0.32. In
general, the CS surfaces that we can detect are elevated, but
span smaller radii due to low abundance and/or low SNR in
the outer disk. We note that this result is dictated by the SNR
we are able to achieve.

The error on each individual *r*, *z* point is dependent on
several factors, mostly stemming from the SNR and the
geometric parameters *M*<sub>*</sub>, *i*, PA, *v*<sub>LSR</sub>, *x*<sub>0</sub>, and *y*<sub>0</sub> (see
S. M. Andrews et al. 2024), derived in A. F. Izquierdo et al.
(2025). Of course, a lower SNR will lead to higher errors, due
to the increased difficulty of correctly identifying isovelocity
contours. The effect of low SNR will lead to scatter in the
surface, or nondetection entirely. Of the geometric parameters,
inclination is one of the largest contributors to the error of the
emission surfaces. By nature, disks that are more face-on will
have a larger scatter in the *z*-direction, due to the fact that the
front and back surfaces become less separated from our
viewing angle. The molecular layers of low-inclination sources
will therefore also be less distinguishable, causing layers to
appear to be overlapping. Disks with a large-scale non-
Keplerian motion also introduce an error, as this emission does
not follow the traditional “butterfly” pattern seen in channel
maps, and can lead to emission being incorrectly vertically
classified. While we do not report explicit errors on each *r*, *z*
point, the amount of scatter in the surface reflects the
confidence of the emission surfaces. This scatter can be seen
as the gray background points in Figure 2 and Figures 9, 10,
and 11 in Appendix A. This is also reflected in the upper and
lower shaded regions of the rolling surfaces, which represent
the standard deviation within each bin, shown in Figure 12 in
Appendix A. In particular, the surfaces of HD 143006, J1852,
MWC 758, and PDS 66 exhibit a large vertical scatter due to
either a face-on orientation or a large-scale non-Keplerian
motion, or a combination of the both.

4. Temperature Radial Profiles and 2D Structure

In Section 4.1, we describe how the gas disk size was
determined. In Section 4.2, we outline the methods used to
obtain the azimuthally averaged temperature profiles and the
gas temperatures for each emission surface point. In
Section 4.3, we discuss how we fit an analytical model to the
temperatures to obtain a 2D temperature distribution for
each disk.

4.1. Gas Disk Size

We compute the size of the gaseous disk utilizing the integrated intensity radial profiles of ¹²CO via a procedure that sums the emission so that the total flux can be determined at various radii, broadly consistent with literature methods (M. Ansdell et al. 2018; L. Trapman et al. 2023). When first calculating the radial profiles, we do not define an outer cutoff radius and instead rely on a hand-selected field of view that encompasses all the emission, plus an extent of a few arcseconds with just noise. We then take a cumulative sum of the intensity, $F_{\text{cuml}}$, which is calculated as a function of the radius (see Equation (A.1) in E. Sanchis et al. 2020). This cumulative sum is then normalized by its maximum flux, $F_{\text{max}}$.
We define an outer edge, $R_{edge}$, as the first peak of $F_{cuml}/F_{max}$, found using scipy with the peak.finder function (P. Virtanen et al. 2020). We make a radial cut of all the intensity past $R_{edge}$; we use the remaining intensity to calculate the total flux, $F_{tot}$. We additionally calculate $R_{95}$, $R_{90}$, and $R_{68}$ by taking the 95th, 90th, and 68th percentiles of the total flux. We report $R_{edge}$ and $R_{90}$ in Table 2. The flux here differs slightly from that reported in R. Teague et al. (2025) due to differences in methods, but our results are broadly consistent. Figure 3 demonstrates this procedure for LkCa 15.

## 4.2. Brightness Temperatures

We use the radial_profile function from GoFish (R. Teague 2019) to obtain the integrated and peak intensity radial profiles from the 3D data cubes in Jy beam⁻¹ km s⁻¹ and kelvins, respectively, assuming the full Planck function. This method shifts and stacks the spectra, utilizing the known rotation of the disk, to better resolve the line profile (see also R. Teague et al. 2016; H.-W. Yen et al. 2016). We assume the geometric parameters and parametric surfaces as determined by DISCMINER for each fit; the differences between using the DISCMINER surface versus the disksurf surface are negligible. For the integrated intensity radial profiles, we use the continuum-subtracted cubes for all three molecules. For the peak intensity radial profiles, we use the non-continuum-subtracted fiducial image cubes for ¹²CO and ¹³CO. E. Weaver et al. (2018) finds that continuum subtraction leads to an underestimation of the gas temperature when the line is optically thick and absorbs the underlying continuum. Therefore, we opt to use the non-continuum-subtracted cubes as to not underestimate the temperatures. The radial profiles from GoFish are determined by azimuthally averaging radial bins that are binned by one-quarter of the FWHM of the synthesized beam. We report uncertainties of the radial profiles that account for the noise in each channel as described in H.-W. Yen et al. (2016).

We also obtain peak intensity radial profiles from the emission surface data. Each r, z point found with disksurf has an associated emission intensity, $I$, in Jy beam⁻¹. We obtain the intensity for each of these r-z points, and convert that to a brightness temperature using the full Planck function:

$$T_B = \frac{h\nu}{k_B} \left[ \ln\left(\frac{2h\nu^3}{c^2I_\nu} + 1\right) \right]^{-1}. \quad (2)$$

In Equation (2), h is the Planck constant, $\nu$ is the frequency, $k_B$ is the Stefan-Boltzmann constant, c is the speed of light, and $I_\nu$ is the emission intensity. Again, we opt to use the non-continuum-subtracted cubes for surface extraction to properly estimate the temperatures. We have repeated the emission surface analysis with the continuum-subtracted cubes and do not find statistically significant differences between them. For a more in-depth discussion on the imaging procedure of continuum subtraction on the exoALMA sources, see R. Teague et al. (2025) and R. A. Loomis et al. (2025).

The resulting peak intensity radial profiles for each molecule are shown as solid lines (for GoFish values) and scatter points (for disksurf values) in Figure 4. The error bars on these points represent the standard deviation of the temperature within each bin. The integrated intensity profiles are shown in Appendix C in Figures 15, 16, and 17. We also compare radial profiles derived directly from DISCMINER and find they are

broadly consistent with those derived with go_fish and disksurf; Figure 18 in Appendix C shows the comparison.

The brightness temperatures obtained via GoFish and disksurf are in good agreement (Figure 4). Nearly all of the binned data points from disksurf are within the error bounds of the peak intensities from GoFish. Some disks have slight discrepancies between GoFish and disksurf peak intensities in the outer regions of the disks (e.g., AA Tau, DM Tau, HD 34282, and J1842.) This is likely due to disksurf tracing discrete vertical and radial locations within the disk, whereas the shifted and staked radial profiles could face contamination from the backside more easily in the outer disk. Even so, these differences are slight and generally under 5 K.

Figure 5 shows the emission surfaces with each respective point colored by its brightness temperature. The background contours of this plot will be discussed in Section 4.3. The smaller background points are the raw r-z surfaces. All three molecules have been binned separately to distinguish between the layers, and the color of the binned points represents the temperature of that specific radial bin. Points below 20 K are denoted by a triangle. From Figure 5, we can see how the temperature decreases as we probe deeper vertical layers and further radial extents. In each disk, we see that the surface layer probed by ¹²CO is the hottest, whereas the ¹³CO surfaces are generally below 30 K, and the CS surface points are almost all below 20 K. For DM Tau, the 0.730 beam size cubes were used for the emission surface extraction, causing the measured temperatures to be colder due to beam dilution.

The mean ¹²CO brightness temperature at 100 au across all disks is ~40 K ± 7 K. For ¹³CO, the mean brightness temperature at 100 au is ~26 K ± 8 K, and for CS, this value is ~12 K ± 5 K. The brightness temperatures of ¹³CO are generally colder than that of ¹²CO, and the same is true for ¹³CO and CS. The average ¹³CO temperature is ~22 K, with much of the emission falling below 20 K. If we assume a CO freeze-out temperature of 20 K, this implies that at least some of the ¹³CO emission is optically thin. Nearly all of the observed CS $J=7-6$ emission is colder than that of ¹³CO, falling below 20 K. The one exception is the outer CS disk of PDS 66, which becomes slightly hotter than ¹³CO at ~75 au. The freeze-out temperature of CS in disks is ≈50-60 K (R. T. Garrod & E. Herbst 2006; G. van der Plas et al. 2014), which implies that nearly all of the CS is optically thin.

## 4.3. 2D Temperature Structure

With multiple molecular tracers probing different heights across the disk, we can construct the 2D temperature structure of the disk in an r-z plane (see E. Dartois et al. 2003; K. A. Rosenfeld et al. 2013; C. P. Dullemond et al. 2020; C. J. Law et al. 2021a; M. Leemker et al. 2022; C. J. Law et al. 2023, 2024). If we assume that the CO lines are optically thick and in local thermodynamic equilibrium, we can assume that the measured brightness temperature $T_B$ is equal to the gas temperature. In this analysis, we adopt the two-layer model connecting the atmosphere temperature to the midplane temperature with a sinusoidal equation originally used in E. Dartois et al. (2003). We tested other functional forms, including the hyperbolic tangent function used in C. P. Dullemond et al. (2020) and C. J. Law et al. (2021a). We find that the sinusoidal function described below best reproduces the observations and gives more reasonable estimates of the midplane temperatures with fewer discrepancies between the
<table><caption>**Table 2**
Outer Disk Radii, R<sub>edge</sub>, and the Radii, Which Encompass 95%, 90%, and 68% of the Total Flux</caption><thead><tr><th>Source</th><th>Line</th><th>R<sub>edge</sub><br>(au)</th><th>R<sub>95</sub><br>(au)</th><th>R<sub>90</sub><br>(au)</th><th>R<sub>68</sub><br>(au)</th></tr></thead><tbody><tr><td rowspan="3">DM Tau</td><td><sup>12</sup>CO J = 3-2</td><td>991 ± 92</td><td>856<sup>+54</sup><sub>-49</sub></td><td>786<sup>+32</sup><sub>-27</sub></td><td>580<sup>+27</sup><sub>-32</sub></td></tr><tr><td><sup>13</sup>CO J = 3-2</td><td>532 ± 22</td><td>483<sup>+27</sup><sub>-22</sub></td><td>462<sup>+16</sup><sub>-22</sub></td><td>310<sup>+22</sup><sub>-22</sub></td></tr><tr><td>CS J = 7-6</td><td>143 ± 22</td><td>121<sup>+16</sup><sub>-32</sub></td><td>111<sup>+27</sup><sub>-27</sub></td><td>78<sup>+5</sup><sub>-11</sub></td></tr><tr><td rowspan="3">AA Tau</td><td><sup>12</sup>CO J = 3-2</td><td>644 ± 10</td><td>498<sup>+141</sup><sub>-56</sub></td><td>442<sup>+56</sup><sub>-51</sub></td><td>265<sup>+20</sup><sub>-20</sub></td></tr><tr><td><sup>13</sup>CO J = 3-2</td><td>331 ± 20</td><td>265<sup>+61</sup><sub>-30</sub></td><td>235<sup>+30</sup><sub>-20</sub></td><td>164<sup>+10</sup><sub>-5</sub></td></tr><tr><td>CS J = 7-6</td><td>336 ± 141</td><td>285<sup>+45</sup><sub>-20</sub></td><td>265<sup>+20</sup><sub>-45</sub></td><td>154<sup>+5</sup><sub>-10</sub></td></tr><tr><td rowspan="3">LkCa 15</td><td><sup>12</sup>CO J = 3-2</td><td>964 ± 12</td><td>799<sup>+53</sup><sub>-59</sub></td><td>699<sup>+41</sup><sub>-29</sub></td><td>456<sup>+18</sup><sub>-18</sub></td></tr><tr><td><sup>13</sup>CO J = 3-2</td><td>946 ± 18</td><td>634<sup>+65</sup><sub>-83</sub></td><td>533<sup>+35</sup><sub>-24</sub></td><td>386<sup>+12</sup><sub>-24</sub></td></tr><tr><td>CS J = 7-6</td><td>339 ± 206</td><td>309<sup>+24</sup><sub>-35</sub></td><td>292<sup>+29</sup><sub>-88</sub></td><td>127<sup>+12</sup><sub>-12</sub></td></tr><tr><td rowspan="3">HD 34282</td><td><sup>12</sup>CO J = 3-2</td><td>840 ± 0</td><td>678<sup>+151</sup><sub>-81</sub></td><td>597<sup>+81</sup><sub>-46</sub></td><td>422<sup>+23</sup><sub>-23</sub></td></tr><tr><td><sup>13</sup>CO J = 3-2</td><td>886 ± 23</td><td>701<sup>+174</sup><sub>-185</sub></td><td>562<sup>+174</sup><sub>-104</sub></td><td>342<sup>+35</sup><sub>-23</sub></td></tr><tr><td>CS J = 7-6</td><td>388 ± 162</td><td>342<sup>+35</sup><sub>-35</sub></td><td>319<sup>+35</sup><sub>-23</sub></td><td>261<sup>+12</sup><sub>-23</sub></td></tr><tr><td rowspan="3">MWC 758</td><td><sup>12</sup>CO J = 3-2</td><td>406 ± 0</td><td>301<sup>+99</sup><sub>-47</sub></td><td>266<sup>+53</sup><sub>-29</sub></td><td>172<sup>+18</sup><sub>-18</sub></td></tr><tr><td><sup>13</sup>CO J = 3-2</td><td>278 ± 12</td><td>219<sup>+53</sup><sub>-58</sub></td><td>172<sup>+58</sup><sub>-23</sub></td><td>114<sup>+12</sup><sub>-12</sub></td></tr><tr><td>CS J = 7-6</td><td>403 ± 58</td><td>357<sup>+35</sup><sub>-70</sub></td><td>333<sup>+58</sup><sub>-94</sub></td><td>146<sup>+105</sup><sub>-35</sub></td></tr><tr><td rowspan="3">CQ Tau</td><td><sup>12</sup>CO J = 3-2</td><td>277 ± 0</td><td>204<sup>+67</sup><sub>-78</sub></td><td>170<sup>+101</sup><sub>-56</sub></td><td>109<sup>+34</sup><sub>-34</sub></td></tr><tr><td><sup>13</sup>CO J = 3-2</td><td>277 ± 11</td><td>221<sup>+50</sup><sub>-95</sub></td><td>165<sup>+106</sup><sub>-50</sub></td><td>92<sup>+11</sup><sub>-11</sub></td></tr><tr><td>CS J = 7-6</td><td>131 ± 129</td><td>109<sup>+17</sup><sub>-17</sub></td><td>103<sup>+22</sup><sub>-17</sub></td><td>81<sup>+6</sup><sub>-11</sub></td></tr><tr><td rowspan="3">SY Cha</td><td><sup>12</sup>CO J = 3-2</td><td>734 ± 7</td><td>604<sup>+123</sup><sub>-61</sub></td><td>543<sup>+61</sup><sub>-55</sub></td><td>378<sup>+20</sup><sub>-27</sub></td></tr><tr><td><sup>13</sup>CO J = 3-2</td><td>529 ± 96</td><td>461<sup>+61</sup><sub>-130</sub></td><td>358<sup>+109</sup><sub>-48</sub></td><td>242<sup>+20</sup><sub>-20</sub></td></tr><tr><td>CS J = 7-6</td><td>222 ± 34</td><td>188<sup>+27</sup><sub>-14</sub></td><td>174<sup>+27</sup><sub>-7</sub></td><td>147<sup>+7</sup><sub>-7</sub></td></tr><tr><td rowspan="3">PDS 66</td><td><sup>12</sup>CO J = 3-2</td><td>149 ± 0</td><td>119<sup>+26</sup><sub>-22</sub></td><td>112<sup>+33</sup><sub>-32</sub></td><td>78<sup>+15</sup><sub>-11</sub></td></tr><tr><td><sup>13</sup>CO J = 3-2</td><td>123 ± 4</td><td>97<sup>+22</sup><sub>-44</sub></td><td>83<sup>+37</sup><sub>-33</sub></td><td>50<sup>+18</sup><sub>-15</sub></td></tr><tr><td>CS J = 7-6</td><td>123 ± 22</td><td>112<sup>+7</sup><sub>-11</sub></td><td>105<sup>+11</sup><sub>-7</sub></td><td>86<sup>+4</sup><sub>-4</sub></td></tr><tr><td rowspan="3">HD 135344B</td><td><sup>12</sup>CO J = 3-2</td><td>291 ± 0</td><td>215<sup>+71</sup><sub>-25</sub></td><td>195<sup>+30</sup><sub>-15</sub></td><td>144<sup>+15</sup><sub>-10</sub></td></tr><tr><td><sup>13</sup>CO J = 3-2</td><td>251 ± 5</td><td>175<sup>+71</sup><sub>-20</sub></td><td>159<sup>+30</sup><sub>-15</sub></td><td>114<sup>+10</sup><sub>-15</sub></td></tr><tr><td>CS J = 7-6</td><td>170 ± 20</td><td>154<sup>+10</sup><sub>-35</sub></td><td>139<sup>+25</sup><sub>-30</sub></td><td>99<sup>+5</sup><sub>-10</sub></td></tr><tr><td rowspan="3">HD 143006</td><td><sup>12</sup>CO J = 3-2</td><td>204 ± 0</td><td>160<sup>+38</sup><sub>-44</sub></td><td>147<sup>+50</sup><sub>-38</sub></td><td>103<sup>+19</sup><sub>-25</sub></td></tr><tr><td><sup>13</sup>CO J = 3-2</td><td>173 ± 13</td><td>135<sup>+31</sup><sub>-19</sub></td><td>129<sup>+38</sup><sub>-19</sub></td><td>97<sup>+13</sup><sub>-13</sub></td></tr><tr><td>CS J = 7-6</td><td>135 ± 100</td><td>122<sup>+6</sup><sub>-13</sub></td><td>116<sup>+13</sup><sub>-13</sub></td><td>97<sup>+6</sup><sub>-0</sub></td></tr><tr><td rowspan="3">RXJ1604.3-2130 A</td><td><sup>12</sup>CO J = 3-2</td><td>306 ± 11</td><td>247<sup>+22</sup><sub>-11</sub></td><td>236<sup>+11</sup><sub>-11</sub></td><td>176<sup>+5</sup><sub>-5</sub></td></tr><tr><td><sup>13</sup>CO J = 3-2</td><td>241 ± 11</td><td>214<sup>+22</sup><sub>-16</sub></td><td>192<sup>+22</sup><sub>-11</sub></td><td>138<sup>+11</sup><sub>-5</sub></td></tr><tr><td>CS J = 7-6</td><td>187 ± 195</td><td>149<sup>+33</sup><sub>-11</sub></td><td>144<sup>+11</sup><sub>-11</sub></td><td>122<sup>+0</sup><sub>-5</sub></td></tr><tr><td rowspan="3">RXJ1615.3-3255</td><td><sup>12</sup>CO J = 3-2</td><td>722 ± 0</td><td>553<sup>+99</sup><sub>-35</sub></td><td>506<sup>+29</sup><sub>-23</sub></td><td>365<sup>+18</sup><sub>-12</sub></td></tr><tr><td><sup>13</sup>CO J = 3-2</td><td>477 ± 12</td><td>395<sup>+64</sup><sub>-41</sub></td><td>348<sup>+35</sup><sub>-23</sub></td><td>249<sup>+12</sup><sub>-12</sub></td></tr><tr><td>CS J = 7-6</td><td>325 ± 64</td><td>290<sup>+29</sup><sub>-41</sub></td><td>260<sup>+58</sup><sub>-41</sub></td><td>155<sup>+41</sup><sub>-23</sub></td></tr><tr><td rowspan="3">V4046 Sgr</td><td><sup>12</sup>CO J = 3-2</td><td>481 ± 0</td><td>369<sup>+48</sup><sub>-27</sub></td><td>334<sup>+21</sup><sub>-21</sub></td><td>215<sup>+13</sup><sub>-11</sub></td></tr><tr><td><sup>13</sup>CO J = 3-2</td><td>278 ± 62</td><td>240<sup>+35</sup><sub>-29</sub></td><td>210<sup>+29</sup><sub>-21</sub></td><td>130<sup>+8</sup><sub>-8</sub></td></tr><tr><td>CS J = 7-6</td><td>154 ± 13</td><td>138<sup>+8</sup><sub>-8</sub></td><td>130<sup>+5</sup><sub>-8</sub></td><td>101<sup>+3</sup><sub>-3</sub></td></tr><tr><td rowspan="3">RXJ1842.9-3532</td><td><sup>12</sup>CO J = 3-2</td><td>388 ± 0</td><td>314<sup>+68</sup><sub>-40</sub></td><td>280<sup>+45</sup><sub>-28</sub></td><td>195<sup>+17</sup><sub>-17</sub></td></tr><tr><td><sup>13</sup>CO J = 3-2</td><td>286 ± 11</td><td>258<sup>+23</sup><sub>-51</sub></td><td>224<sup>+45</sup><sub>-40</sub></td><td>133<sup>+17</sup><sub>-11</sub></td></tr><tr><td>CS J = 7-6</td><td>144 ± 23</td><td>127<sup>+11</sup><sub>-28</sub></td><td>110<sup>+28</sup><sub>-17</sub></td><td>88<sup>+11</sup><sub>-11</sub></td></tr><tr><td rowspan="3">RXJ1852.3-3700</td><td><sup>12</sup>CO J = 3-2</td><td>322 ± 0</td><td>256<sup>+61</sup><sub>-44</sub></td><td>223<sup>+44</sup><sub>-28</sub></td><td>140<sup>+11</sup><sub>-11</sub></td></tr><tr><td><sup>13</sup>CO J = 3-2</td><td>240 ± 33</td><td>179<sup>+35</sup><sub>-22</sub></td><td>157<sup>+28</sup><sub>-17</sub></td><td>113<sup>+6</sup><sub>-6</sub></td></tr><tr><td>CS J = 7-6</td><td>157 ± 22</td><td>141<sup>+11</sup><sub>-17</sub></td><td>130<sup>+17</sup><sub>-11</sub></td><td>96<sup>+6</sup><sub>-6</sub></td></tr></tbody></table>

observed temperatures and the 2D fit (see Appendix B and Figure 13). With this functional form, the midplane temperature ($T_{\text{mid}}$) and atmospheric temperature ($T_{\text{atm}}$) are described

by power laws:

$$T_{\text{mid}}(r) = T_{\text{mid},0}(r/r_0)^{q_{\text{mid}}} \quad (3)$$
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Normalized integrated intensity for LkCa 15. The dashed vertical lines show R<sub>edge</sub>, R<sub>95</sub>, R<sub>90</sub>, and R<sub>68</sub>. Only intensity inward of R<sub>edge</sub> is used in the calculation of the total flux and respective 95th, 90th, and 68th percentile radii.</figcaption></figure>

$$T_{\text{atm}}(r) = T_{\text{atm},0}(r/r_0)^{q_{\text{atm}}} \quad (4)$$

where T<sub>mid,0</sub> and T<sub>atm,0</sub> represent the atmospheric and midplane temperature at r<sub>0</sub> = 100 au, and q<sub>mid</sub> and q<sub>atm</sub> are the slopes to the power laws, respectively. At each radius, the temperature is smoothly connected between the midplane and atmosphere using a sinusoidal function:

$$T(z) = \begin{cases} T_{\text{atm}} + (T_{\text{mid}} - T_{\text{atm}})\cos^2\left(\frac{\pi z}{2 z_q}\right) & \text{if } z < z_q \\ T_{\text{atm}} & \text{if } z \ge z_q \end{cases}. \quad (5)$$

Here, z<sub>q</sub> is given by a power law of the form z<sub>q</sub>(r) = z<sub>0</sub>(r/r<sub>0</sub>)<sup>β</sup> where z<sub>0</sub> is the height in the disk at r = r<sub>0</sub>, and β is the power-law exponent. To determine the 2D temperature structure, we fit these equations to the nonparametric surface data for six parameters: T<sub>atm,0</sub>, T<sub>mid,0</sub>, q<sub>atm</sub>, q<sub>mid</sub>, z<sub>0</sub>, and β. In this fit, we discard the data points with brightness temperatures below 20 K because those likely probe optically thin emission and thus do not reflect the true gas temperature. All CS temperatures used in the temperature fitting fall below 20 K, and are therefore not used.

We have constrained the parameter space with priors as follows: T<sub>mid,0</sub> < T<sub>atm,0</sub> < 150 K, -2 < q<sub>atm</sub> < 0, -2 < q<sub>mid</sub> < 0, 0 < z<sub>0</sub> < r<sub>max</sub> where r<sub>max</sub> is the largest radial point, and 0 < β < 2. These priors were informed by previous observations of similar disks and modeling efforts; see E. Dartois et al. 2003, K. Zhang et al. 2021, J. K. Calahan et al. 2021, and C. J. Law et al. 2021a. We attempted to keep the priors listed above as broad as possible as to not impart inherent assumptions on the temperature profiles, but the convergence of the fits was contingent on these priors being in place. The lower bounds of T<sub>mid</sub> are set by the continuum temperatures T<sub>cont</sub> derived in P. Curone et al. (2025). At radial locations where there is no continuum information, we set the lower bound for T<sub>mid</sub> to the cosmic microwave background (CMB) temperature T<sub>CMB</sub> = 2.73 K. We additionally place an upper bound on T<sub>mid</sub> using the expected irradiation temperature given by

$$T_{\text{irr}}(r) = \left( \frac{\phi L_{*}}{8\pi r^2 \sigma_{\text{SB}}} \right)^{1/4} \quad (6)$$

where L<sub>*</sub> is the stellar luminosity, r is the distance from the central star, σ<sub>SB</sub> is the Stefan-Boltzmann constant, and φ is an

irradiation factor that is dependent on the angle that the stellar radiation impinges onto the disk (P. D’Alessio et al. 2001). The midplane temperature bounds are then max(T<sub>cont</sub>, T<sub>CMB</sub>) < T<sub>mid</sub> < T<sub>irr</sub>. We opt for φ = 1 as a conservative upper bound for stellar irradiation-dominated disk; in practice, we find that the posterior midplane temperature is usually much lower than T<sub>irr</sub>.

To fit Equations (4), (3), and (5), we use emcee (D. Foreman-Mackey et al. 2013), adopting 256 walkers, 1000 steps, and 100 burn-in steps. We start with a randomized initial guess within the bounds of the prior conditions to obtain the first fit. We then use the results of this first fit to inform the initial guess of a second round with the same number of walkers, steps, and burn-in steps. We found that this procedure yielded more converged results due to walkers getting stuck in the first round of fitting, but we note that differences between the first and second round of fitting are not generally large. The results of the 2D temperature structure are reported in Table 3 where we report the 50th percentiles as the best-fit values and the 16th and 84th percentiles as upper and lower errors. These represent only the statistical errors and do not include systematic errors, which can be larger than the statistical errors. Because of this, some of the errors reported in Table 3 are zero, particularly for q<sub>mid</sub> and β. For a discussion of the systematic errors, see Appendix B.

The derived 2D temperature profiles are shown in Figure 5 as filled contours. In the figure, the larger circles show the data points binned by 1/4 of the beam size, and the triangles show the binned data that were not utilized in the fitting procedure because it was below the CO freeze-out temperature (<20 K). Lines of constant temperature are drawn at 50, 40, 30, and 20 K. From the fitting procedure, we obtain an estimate of the atmospheric temperature as a function of the radius, T<sub>atm</sub>(r), and the midplane temperature, T<sub>mid</sub>(r). We find atmospheric temperatures (at 100 au) range from ~67 K (HD 34282) to ~34 K (J1615), with a mean atmospheric temperature (T<sub>atm</sub>) at 100 au across all disks of 43.4 K. We find midplane temperatures ranging from 13 K (AA Tau) up to 32 K (HD 34282). The average midplane temperature at 100 au across the disk sample is 24 K. Because there is no optically thick data at the midplane, this value is an extrapolation informed by the priors. Due to this, we stipulate that the reported T<sub>mid</sub> is closer to an upper bound to the true midplane temperature. In J1852 and in V4046 Sgr, the midplane temperature inverts, as shown by the contour at 20 K in Figure 5. Both disks have low inclinations around 35°, which causes highly scattered data and difficulty in separating the molecular layers, particularly in the outer disk. This leads to larger uncertainties in the temperature profiles.

## 5. Substructure

Vertical substructures have been observed in many emission surfaces (see C. J. Law et al. 2021a; C. J. Law et al. 2023; T. Paneque-Carreño et al. 2023). These localized structures have been found to be colocated with substructures identified in the continuum (C. J. Law et al. 2021a; T. Paneque-Carreño et al. 2023), intensity radial profiles (C. J. Law et al. 2021b), or velocity perturbations (R. Teague et al. 2019a, 2021b; M. Galloway-Sprietsma et al. 2023; T. Paneque-Carreño et al. 2023). Their identification can provide insight into finer disk structures; surface modulations are expected to arise during planet formation, which is expected to open gaseous gaps (R. Dong et al. 2015), or produce changes in the
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Peak intensities of ¹²CO *J* = 3–2 (blue), ¹³CO *J* = 3–2 (orange), and CS*J* = 7–6 (pink) in kelvins for all exoALMA sources. The solid black lines show the peak intensities found with GoFish, and the points are the temperature points found using *disksurf*, binned by a quarter of the beam size. The error bars on the points represent the standard deviation of the temperature within each bin. The radial profiles are reported out to the molecular emission *R*<sub>edge</sub> radius.</figcaption></figure>

molecular abundance and chemical stratification (Y. Aikawa & E. Herbst 1999). In the following section, we detail the process of substructure identification for both the emission surfaces and radial intensity profiles, and compare these with those seen in the continuum and the velocity profiles.

## 5.1. Emission Surface Substructures

Some of the most prominent emission surface features easily identified by eye are the sharp cutoff from rising emission seen

in AA Tau, SY Cha, and LkCa 15. In ¹²CO, both surfaces of AA Tau and SY Cha sharply drop at 237 and 359 au, respectively. In ¹³CO, LkCa 15 also exhibits a sharp drop in the emission surface at 274 au. These emission surface morphologies differ from the traditional tapered power law in that the falling portion of the emission is distinctly separated from the rising portion of the surface, and that this outer emission is scattered with no well-defined shape. While *disksurf* has identified these regions as part of the upper emission surface, a
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Emission surface points with associated brightness temperatures; raw surface data points are shown as smaller background points. The raw data is binned by a quarter of a beam size for molecules above 20 K (larger circles) and below 20 K (triangles). The contour background shows the 2D temperature profiles calculated with Equation (5). Lines of constant temperatures at 50, 40, 30, and 20 K are shown in black.</figcaption></figure>

<table><caption>**Table 3**
2D Temperature Structure Fits</caption><thead><tr><th>Source</th><th>T<sub>atm,0</sub><br>(K)</th><th>T<sub>mid,0</sub><br>(K)</th><th>q<sub>atm</sub></th><th>q<sub>mid</sub></th><th>z<sub>0</sub><br>(au)</th><th>β</th></tr></thead><tbody><tr><td>DM Tau</td><td>36<sup>+0.22</sup><sub>-0.3</sub></td><td>23<sup>+0.15</sup><sub>-0.9</sub></td><td>-0.13<sup>+0.015</sup><sub>-0.1</sub></td><td>-0.5<sup>+0.015</sup><sub>-0.1</sub></td><td>25<sup>+0.79</sup><sub>-0.19</sub></td><td>2<sup>+0.0057</sup><sub>-0.93</sub></td></tr><tr><td>AA Tau</td><td>41<sup>+0.43</sup><sub>-0.39</sub></td><td>13<sup>+0.31</sup><sub>-0.33</sub></td><td>-0.51<sup>+0.017</sup><sub>-0.018</sub></td><td>-0.21<sup>+0.017</sup><sub>-0.018</sub></td><td>60<sup>+0.87</sup><sub>-0.8</sub></td><td>0.075<sup>+0.015</sup><sub>-0.015</sub></td></tr><tr><td>LkCa 15</td><td>48<sup>+0.43</sup><sub>-0.4</sub></td><td>20<sup>+0.15</sup><sub>-0.15</sub></td><td>-0.55<sup>+0.0075</sup><sub>-0.0077</sub></td><td>-0.23<sup>+0.0075</sup><sub>-0.0077</sub></td><td>55<sup>+0.69</sup><sub>-0.67</sub></td><td>0.59<sup>+0.0073</sup><sub>-0.0074</sub></td></tr><tr><td>HD 34282</td><td>67<sup>+1.3</sup><sub>-0.78</sub></td><td>32<sup>+0.14</sup><sub>-0.14</sub></td><td>-0.00093<sup>+0.00071</sup><sub>-0.0025</sub></td><td>-0.25<sup>+0.00071</sup><sub>-0.0025</sub></td><td>87<sup>+2.5</sup><sub>-1.6</sub></td><td>0.69<sup>+0.0059</sup><sub>-0.009</sub></td></tr><tr><td>SY Cha</td><td>45<sup>+0.37</sup><sub>-0.35</sub></td><td>24<sup>+0.13</sup><sub>-0.15</sub></td><td>-0.58<sup>+0.0096</sup><sub>-0.01</sub></td><td>-0.3<sup>+0.0096</sup><sub>-0.01</sub></td><td>56<sup>+0.83</sup><sub>-0.82</sub></td><td>0.0069<sup>+0.0096</sup><sub>-0.005</sub></td></tr><tr><td>RXJ1615.3-3255</td><td>34<sup>+0.051</sup><sub>-0.044</sub></td><td>24<sup>+0.12</sup><sub>-0.15</sub></td><td>-0.098<sup>+0.0019</sup><sub>-0.0032</sub></td><td>-0.25<sup>+0.0019</sup><sub>-0.0032</sub></td><td>32<sup>+0.22</sup><sub>-0.23</sub></td><td>1.1<sup>+0.011</sup><sub>-0.018</sub></td></tr><tr><td>V4046 Sgr</td><td>37<sup>+0.036</sup><sub>-0.039</sub></td><td>28<sup>+0.067</sup><sub>-0.067</sub></td><td>-0.62<sup>+0.0022</sup><sub>-0.0021</sub></td><td>-0.35<sup>+0.0022</sup><sub>-0.0021</sub></td><td>10<sup>+0.06</sup><sub>-0.062</sub></td><td>0.00025<sup>+0.00051</sup><sub>-0.00019</sub></td></tr><tr><td>RXJ1842.9-3532</td><td>43<sup>+0.14</sup><sub>-0.13</sub></td><td>25<sup>+0.08</sup><sub>-0.094</sub></td><td>-0.45<sup>+0.0062</sup><sub>-0.0062</sub></td><td>-0.23<sup>+0.0062</sup><sub>-0.0062</sub></td><td>27<sup>+0.25</sup><sub>-0.25</sub></td><td>0.00077<sup>+0.0015</sup><sub>-0.00058</sub></td></tr><tr><td>RXJ1852.3-3700</td><td>40<sup>+0.092</sup><sub>-0.088</sub></td><td>30<sup>+0.15</sup><sub>-0.16</sub></td><td>-0.87<sup>+0.0072</sup><sub>-0.0071</sub></td><td>-0.36<sup>+0.0072</sup><sub>-0.0071</sub></td><td>16<sup>+0.2</sup><sub>-0.2</sub></td><td>0.00081<sup>+0.0019</sup><sub>-0.00061</sub></td></tr></tbody></table>

visual inspection and line profile analysis has shown that they are likely a combination of emission from the upper and lower surfaces, as well as the midplane. For this reason, we exclude these points when fitting the surface with *disksurf*. The physical reasons for this ambiguous “diffuse” emission are unclear. It may be a projection effect, given that all three disks are highly inclined at around 50°. It may imply that the disks have different sizes between the upper and lower surfaces, or that the outer region of the upper surface becomes colder than that of the backside. It could also be due to photodesorption in

the outer regions of the disks. Sections 6.3 and 6.4.1 discuss these possibilities in more detail.

For smaller-scale substructure identification, we rely on the rolling surfaces derived using *disksurf* shown in Figure 12 in Appendix A. The rolling surfaces shown in this figure have been smoothed with a Savitzky–Golay filter for visual clarity, with a window length of 1" and a polynomial order of 2; surface perturbations smaller than the beam size are smoothed out. To identify surface structure, we utilize the *findPeaks* function from *scipy* (P. Virtanen et al. 2020). This function
finds maxima or minima in 1D data sets by comparison of neighboring values. To ensure the substructures are physical, we set the minimum width of each substructure and the minimum distance between substructures to be at least the beam size (0.715), up to 2 times the beam size (0.73). We scale the maximum height, or prominence, of the peak to the beam size or larger. After running *findPeaks*, we visually confirm by eye the validity of the substructures, as this method can be sensitive to small perturbations, even with the imposed limits. We report each vertical substructure dip in Table 4 and in Figure 12 with dashed lines. Each substructure is labeled with “Z,” in the nomenclature of C. J. Law et al. (2021a), followed by the radial location. For visual clarity, we show only the structures identified as dips, and not the peaks.

We find evidence of localized substructures in 76% of disks with viable emission surfaces. Of these, seven of the substructures are in $^{12}\text{CO}$ emission. Nine of the disks also have vertical substructure in $^{13}\text{CO}$ emission, which is generally colocated with that seen in $^{12}\text{CO}$. Only two definite substructures were identified in CS, although several slope changes in the surfaces are observed. The emission surfaces extracted using the iterative methodology described in Section 3 are generally smoother than the initial surfaces found with just one iteration; we have verified that all surface substructures presented in Figure 12 are consistent between the iterative and noniterative surfaces. We do not find localized surface perturbations large enough to warrant classification in HD 34282 or PDS 66.

## 5.2. Radial Peak Intensity Substructures

We repeat the structure identification process outlined in the previous section, but for the radial peak intensity profiles, smoothed using a Savitzky-Golay filter with a window length of 5 and a polynomial order of 3. To identify structure, we again utilize the peak-finding function from *scipy* (P. Virtanen et al. 2020). We again set the minimum radial width of each substructure and the minimum distance between substructures to be at least the beam size. We scale the maximum height, or prominence, of the peak to the temperature, ranging from 7 to 15 K for some disks. We find peak intensity structure in the disks of AA Tau, CQ Tau, HD 34282, J1615, J1842, LkCa 15, and SY Cha, visualized in Figure 14 in Appendix C. The profiles of HD 135344B, HD 143006, J1604, J1852, MWC 758, and PDS 66 are generally smooth, and we do not characterize any substructure for these disks. We do not characterize the first peak in the intensity profiles, as this is an expected feature present in all disks due to beam dilution in the inner radial regions causing a drop in intensity. We also do not characterize any structure interior to this first peak.

AA Tau, LkCa 15, and SY Cha disks exhibit a wavelike pattern in their $^{12}\text{CO}$ intensity profiles. In AA Tau, there are two slight bright bumps at ~120 au and ~195 au. The wave pattern in LkCa 15 is distinct; there is an initial bright spot at 174 au, also seen in $^{13}\text{CO}$, with subsequent dips and bright spots at nearly evenly spaced intervals, ending at a dip identified at 781 au. The wave patterns are spaced roughly 175 au apart. SY Cha exhibits two bright peaks in a distinct pattern as well, at 276 and 447 au, which are observed to a lesser extent in the $^{13}\text{CO}$ peak intensities. The potential causes of this structure are discussed in Section 6.4.

CQ Tau, HD 34282, J1615, J1842, and V4046 Sgr exhibit singular bright spots in multiple molecules. In CQ Tau, a bright

“shoulder” is seen at 125 au; the shape of the $^{13}\text{CO}$ profile is similar, but no definite peak was identified via our process. HD 34282 displays a prominent bright shoulder in $^{12}\text{CO}$ as well, at 295 au, but the remaining molecules do not show any structure. J1615 exhibits a bright shoulder in $^{12}\text{CO}$ at 366 au, and there is a colocated bright spot radially inward at 132 au for $^{13}\text{CO}$ and CS. In J1842, we identified a bright shoulder at 88 au that is colocated between all three molecules. Finally, V4046 Sgr exhibits a bright shoulder only in CS emission, at 98 au.

We do not find a significant correlation between substructures identified in the peak intensity and those identified in the emission surfaces. One exception is LkCa 15, whose surface dips line up within a few astronomical units of some of the “wiggles” seen in the peak intensity.

## 5.3. Comparison with Dust Continuum

We compare the emission surface substructure radial locations with the rings and gaps in the continuum reported by P. Curone et al. (2025). Unlike the results of the MAPS emission surface study (C. J. Law et al. 2021a), we do not find significant correlation between the molecular emission surfaces and the dust millimeter continuum rings and gaps for our larger sample size. We find that the disks of AA Tau, J1615, and MWC 758 are the only disks whose molecular emission surface substructure potentially correlates with that of the continuum. We present these structures, along with their continuum counterparts, in Figure 6, with areas of interest highlighted in yellow. For AA Tau and J1615, this correlation is only present in the CS $J=7-6$ emission, while MWC 758 exhibits a substructure match in $^{12}\text{CO}$ and $^{13}\text{CO}$.

In AA Tau, the dust gap at 80 au, and the subsequent rings surrounding it, align with a slight dip in the CS surface, visualized in the leftmost panel of Figure 6. We do not see this same correlation in $^{12}\text{CO}$ or $^{13}\text{CO}$; there is no gap structure in either line near 80 au. The correlation between the dust substructure and CS emission in J1615 is slight. J1615 exhibits a dust gap 125.5 au, and at this point, the CS surface dips slightly and then experiences a slope change, visualized in the middle panel of Figure 6.

MWC 758 has the clearest correlation between emission surface structure and dust rings/gaps. MWC 758 displays two clear dips in the $^{12}\text{CO}$ and $^{13}\text{CO}$ surfaces at 56 and 62 au, respectively. The emission surface dip in $^{13}\text{CO}$ aligns nearly perfectly with a continuum gap at 61.3 au. Additionally, the $^{12}\text{CO}$ surface appears to increase in height slightly at the innermost dust ring, around 50 au. The continuum rings and gaps overlaid with the surfaces that have a substructure match are shown in the rightmost panel of Figure 6. These surface dips may be due to the continuum subtraction, but a physical mechanism may be at play given that we do not see them at the locations of the other continuum gaps.

C. J. Law et al. (2024) note a slope change in $^{12}\text{CO}$ $J=3-2$ emission present in PDS 70 at the location of the dust cavity wall, which is reminiscent of what we see in AA Tau, J1615, and MWC 758. Of these disks, AA Tau and MWC 758 are transition disks, and the observed slope changes may be due to the decreased gas surface density within the cavity, as seen in C. J. Law et al. (2024). A more thorough investigation of the inner radial regions of gas in exoALMA sources and their relation to dust cavity walls is reserved for future work.

The substructures identified in the peak intensity align with several continuum rings and gaps. Figure 14 in Appendix C
<table><caption>**Table 4**
Emission Surface Substructures (Gaps Only) and Peak Intensity Substructures (Peaks Only)</caption><thead><tr><th>Source</th><th>Line</th><th>Emission Surface Dips<br>(au)</th><th>Peak Intensity Bumps<br>(au)</th><th>Peak Intensity Dips<br>(au)</th></tr></thead><tbody><tr><td rowspan="3">DM Tau</td><td><sup>12</sup>CO J = 3–2</td><td>120 ± 22, 183 ± 22, 293 ± 22</td><td>...</td><td>...</td></tr><tr><td><sup>13</sup>CO J = 3–2</td><td>194 ± 22</td><td>127 ± 22</td><td>89 ± 22</td></tr><tr><td>CS J = 7–6</td><td>...</td><td>...</td><td>...</td></tr><tr><td rowspan="3">AA Tau</td><td><sup>12</sup>CO J = 3–2</td><td><b>235 ± 20</b></td><td>93 ± 20, 187 ± 20</td><td>86 ± 20</td></tr><tr><td><sup>13</sup>CO J = 3–2</td><td>...</td><td>119 ± 20</td><td>86 ± 20</td></tr><tr><td>CS J = 7–6</td><td>90 ± 20</td><td>...</td><td>...</td></tr><tr><td rowspan="3">LkCa 15</td><td><sup>12</sup>CO J = 3–2</td><td>200 ± 24, 342 ± 24</td><td>176 ± 24, 333 ± 24, 504 ± 24, 651 ± 24</td><td>302 ± 24, 451 ± 24, 604 ± 24, 781 ± 24</td></tr><tr><td><sup>13</sup>CO J = 3–2</td><td>119 ± 24, <b>282 ± 24</b></td><td>203 ± 24</td><td>174 ± 24</td></tr><tr><td>CS J = 7–6</td><td>177 ± 24</td><td>192 ± 24</td><td>...</td></tr><tr><td rowspan="3">HD 34282</td><td><sup>12</sup>CO J = 3–2</td><td>...</td><td>295 ± 46</td><td>214 ± 46</td></tr><tr><td><sup>13</sup>CO J = 3–2</td><td>...</td><td>...</td><td>...</td></tr><tr><td>CS J = 7–6</td><td>...</td><td>...</td><td>...</td></tr><tr><td rowspan="3">MWC 758</td><td><sup>12</sup>CO J = 3–2</td><td>56 ± 23</td><td>...</td><td>...</td></tr><tr><td><sup>13</sup>CO J = 3–2</td><td>62 ± 23</td><td>...</td><td>...</td></tr><tr><td>CS J = 7–6</td><td>62 ± 23</td><td>...</td><td>...</td></tr><tr><td rowspan="3">CQ Tau</td><td><sup>12</sup>CO J = 3–2</td><td>101 ± 22</td><td>125 ± 22</td><td>...</td></tr><tr><td><sup>13</sup>CO J = 3–2</td><td>90 ± 22</td><td>...</td><td>...</td></tr><tr><td>CS J = 7–6</td><td>...</td><td>...</td><td>...</td></tr><tr><td rowspan="3">SY Cha</td><td><sup>12</sup>CO J = 3–2</td><td>220 ± 27, <b>354 ± 27</b></td><td>276 ± 27, 447 ± 27</td><td>242 ± 27, 406 ± 27</td></tr><tr><td><sup>13</sup>CO J = 3–2</td><td>...</td><td>150 ± 27</td><td>113 ± 27, 222 ± 27</td></tr><tr><td>CS J = 7–6</td><td>...</td><td>125 ± 27</td><td>...</td></tr><tr><td rowspan="3">PDS 66</td><td><sup>12</sup>CO J = 3–2</td><td>...</td><td>...</td><td>...</td></tr><tr><td><sup>13</sup>CO J = 3–2</td><td>...</td><td>...</td><td>...</td></tr><tr><td>CS J = 7–6</td><td>...</td><td>...</td><td>...</td></tr><tr><td rowspan="3">HD 135344B</td><td><sup>12</sup>CO J = 3–2</td><td>...</td><td>...</td><td>...</td></tr><tr><td><sup>13</sup>CO J = 3–2</td><td>...</td><td>...</td><td>...</td></tr><tr><td>CS J = 7–6</td><td>...</td><td>...</td><td>...</td></tr><tr><td rowspan="3">HD 143006</td><td><sup>12</sup>CO J = 3–2</td><td>...</td><td>...</td><td>...</td></tr><tr><td><sup>13</sup>CO J = 3–2</td><td>112 ± 25</td><td>...</td><td>...</td></tr><tr><td>CS J = 7–6</td><td>...</td><td>...</td><td>...</td></tr><tr><td rowspan="3">RXJ1604.3-2130 A</td><td><sup>12</sup>CO J = 3–2</td><td>...</td><td>...</td><td>...</td></tr><tr><td><sup>13</sup>CO J = 3–2</td><td>...</td><td>...</td><td>...</td></tr><tr><td>CS J = 7–6</td><td>...</td><td>...</td><td>...</td></tr><tr><td rowspan="3">RXJ1615.3-3255</td><td><sup>12</sup>CO J = 3–2</td><td>472 ± 23</td><td>366 ± 23</td><td>...</td></tr><tr><td><sup>13</sup>CO J = 3–2</td><td>250 ± 23</td><td>132 ± 23</td><td>108 ± 23</td></tr><tr><td>CS J = 7–6</td><td>113 ± 23</td><td>132 ± 23</td><td>120 ± 23</td></tr><tr><td rowspan="3">V4046 Sgr</td><td><sup>12</sup>CO J = 3–2</td><td>173 ± 11, 318 ± 11</td><td>...</td><td>...</td></tr><tr><td><sup>13</sup>CO J = 3–2</td><td>174 ± 11</td><td>...</td><td>...</td></tr><tr><td>CS J = 7–6</td><td>91 ± 11</td><td>98 ± 11</td><td>42 ± 11</td></tr><tr><td rowspan="3">RXJ1842.9-3532</td><td><sup>12</sup>CO J = 3–2</td><td>138 ± 23</td><td>88 ± 23</td><td>59 ± 23</td></tr><tr><td><sup>13</sup>CO J = 3–2</td><td>83 ± 23, 133 ± 23</td><td>88 ± 23</td><td>65 ± 23</td></tr><tr><td>CS J = 7–6</td><td>...</td><td>...</td><td>...</td></tr><tr><td rowspan="3">RXJ1852.3-3700</td><td><sup>12</sup>CO J = 3–2</td><td>117 ± 22, 161 ± 22</td><td>...</td><td>...</td></tr><tr><td><sup>13</sup>CO J = 3–2</td><td>124 ± 22</td><td>...</td><td>...</td></tr><tr><td>CS J = 7–6</td><td>...</td><td>...</td><td>...</td></tr></tbody></table>

**Note.** The bolded substructures represent the radial location of the larger emission surface drops seen in AA Tau, SyCha, and LkCa 15, discussed in Section 6.3. Errors are generically defined as the beam size (0.′/15). We exclude HD 135344 B and J1604, which we do not derive emission surfaces for.

shows the radial profiles of the peak intensities, with the dust gaps as black dashed lines and the dust rings as black solid lines. In AA Tau, we see a dip in $^{13}$CO that aligns within a beam size of the dust gap at 79 au. We see similar structure in

J1615, where gaps identified in $^{13}$CO and CS align with the dust gap at 125 au. This is also seen in J1842; dips in the intensity in $^{12}$CO and $^{13}$CO align with the dust gap at 63 au. There are several instances where dust rings align with the
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Rolling emission surfaces (CS *J* = 7–6, pink; $^{12}$CO *J* = 3–2, blue; $^{13}$CO *J* = 3–2, orange) with continuum substructures identified by P. Curone et al. (2025). The dashed lines show the continuum gaps, and the solid lines show the continuum rings. The yellow shaded regions highlight the surface structures that overlap with continuum and or gas substructures.</figcaption></figure>

brightest intensity peak, as well as substructures interior to this
peak, but many of these are within one beam, and we exclude
them from this analysis.

5.4. *Comparison with Velocity Signatures*

Several of the exoALMA sources display kinematic
perturbations in the individual channel maps, which are
discussed and characterized in C. Pinte et al. (2025);
azimuthally averaged velocity perturbations are derived in
J. Stadler et al. (2025). The signatures characterized as
“velocity kinks” may be indicative of embedded protoplanets
(C. Pinte et al. 2018, 2023), although other explanations are
possible (G. Lodato et al. 2023). We find that surface
modulations in AA Tau, J1842, and LkCa 15 are radially
coincident within a beam size of the proposed planet locations
listed in Table 2 of C. Pinte et al. (2025). The radial locations of
the proposed planets are plotted as green lines in Figure 12 in
Appendix A. In AA Tau, the velocity kink aligns with a dip in
the CS surface at 95 au, which is also within a beam size of the
dust continuum feature at 80 au shown more clearly in
Figure 6. In J1842, the proposed planet location of 150 au
aligns with two dips in the emission surfaces, coincident
between $^{12}$CO and $^{13}$CO. In LkCa 15, the most prominent
surface modulation in $^{12}$CO aligns with the proposed planet
location of 370 au.

In AA Tau, J1615, LkCa 15, and SY Cha, peak intensity
structures align with the proposed planet locations. The dip
identified in $^{13}$CO at 86 au aligns with the proposed planet
location of 85 au, which also coaligns with a continuum gap.
For J1615, the bright peak intensity shoulder at 336 au aligns
with the proposed planet at 370 au. In LkCa 15, the second
“wave” peak is within a beam size of the proposed planet
location. In SY Cha, a $^{13}$CO peak at 150 au aligns with the
proposed planet location.

6. Discussion

6.1. *Comparison with Source Properties*

Numerous studies have shown that emission height is
correlated with some disk and stellar properties. C. J. Law
et al. (2021a) showed that the disk-averaged ⟨z/r⟩ is positively
correlated with the gas disk size. This trend was further verified
in C. J. Law et al. (2022, 2023) and T. Paneque-Carreño et al.
(2023). Relationships between ⟨z/r⟩, stellar mass, and gas
temperature were also seen. Most recently, T. Paneque-Carreño
et al. (2025) demonstrated the correlation between emission
height, CO depletion, and disk mass. Here, we build off of the

data compiled in C. J. Law et al. (2022, 2023) and T. Paneque-Carreño et al. (2025), adding the exoALMA sources, to further explore these correlations.

Plots showing $^{12}$CO emission surface $\langle z/r \rangle$ versus the stellar mass, average gas temperature, and the gas disk size $R_{\text{edge}}$ are shown in the top panels of Figure 7. Here, we focus our discussion to $^{12}$CO emission surface and briefly comment on $^{13}$CO emission surface. The $\langle z/r \rangle$ calculation is described in Section 3.3; we use methods consistent with those used for the compiled background points. We find consistent results with those observed in C. J. Law et al. (2022, 2023). The emission surface $\langle z/r \rangle$ weakly decreases with both the stellar mass and average $^{12}$CO temperature. To explain these trends, C. J. Law et al. (2022) invokes the equation relating the gas pressure scale height with the stellar mass and disk temperature (see Equation (3) in C. J. Law et al. 2022). Assuming that the emission surfaces scale with gas pressure scale heights and that the stellar luminosity and mass are related following $L_* \propto M_*^3$, as in C. J. Law et al. (2022), $\langle z/r \rangle \sim M_*^{-1/8}$ and $\langle z/r \rangle \sim T^{-1/6}$. We plot these relations on the upper left and upper middle panels in Figure 7. We quantify the correlation between these variables using Spearman's rank correlation coefficient ($r_s$), which is used for nonlinear relations. For $\langle z/r \rangle$ versus the stellar mass, we find $r_s = -0.16$. For $\langle z/r \rangle$ versus average brightness temperature, we find $r_s = -0.13$. Overall, we find a trend that $\langle z/r \rangle$ decreases with the stellar mass and temperature although we note that these trends are highly scattered, and not reliably correlated, similar to what C. J. Law et al. (2023) found (see also T. Paneque-Carreño et al. 2023).

Next, we look at the correlation between $\langle z/r \rangle$ and the gas
disk size. We find a linear correlation between them as shown
in the top rightmost panel of Figure 7. To keep consistency
with previous studies (e.g., C. J. Law et al. 2022, 2023), we use
*linmix* (B. C. Kelly 2007) to fit a line to the data. We
find a linear relation quantified by $\langle z/r \rangle = 0.090^{+0.030}_{-0.030} +
0.00033^{+0.000053}_{-0.000052} \times R_{\text{edge}}$. This is steeper than previous relations
found in C. J. Law et al. (2022, 2023), but overall consistent.
We use the Pearson correlation coefficient (PCC, or $r_p$) to
quantify the linear correlation between $\langle z/r \rangle$ and $R_{\text{CO}}$. We find
a PCC of $r_p = 0.67$, indicating a strong positive correlation. We
repeated the $\langle z/r \rangle$ measurement, but only taking emission
surface points at 125 au for each disk, $\pm$ the beam size of 0.9'15,
essentially measuring the $\langle z/r \rangle$ at the same radius for each disk.
We found that the measured $\langle z/r \rangle$ values did not change
significantly, and still resulted in a correlation with $R_{\text{CO}}$. We
note that HD 143006 (labeled as point 5) is an outlier, likely
due to the face-on nature of this disk and subsequent associated
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Top panel: Average $z/r$ in $^{12}$CO for disks with moderate inclination, plotted against $M_*$ (leftmost plot), average $^{12}$CO brightness temperature (middle plot), and the outer radius of the gaseous disk (rightmost plot). The points are colored by brightness temperature at 100 au. Background gray points are literature values from previous studies. Bottom panel: Average $z/r$ in $^{12}$CO for disks with moderate inclination, plotted against disk mass calculated via the rotation curve from C. Longarini et al. (2025; leftmost plot), disk mass calculated via the dust mass from P. Curone et al. (2025; middle plot), and disk mass calculated using $N_2H^+$ from L. Trapman et al. (2025; rightmost plot). Points that represent exoALMA sources are colored by the $^{12}$CO $R_{\text{edge}}$. Smaller gray points are MAPS sources with mass measurements from P. Martire et al. (2024), S. M. Andrews et al. (2018), and L. Trapman et al. (2025). The Pearson correlation coefficient, $r_p$, or the Spearman correlation coefficient, $r_s$, is shown in the top right of each figure. When a linear function is fit, the results are plotted, and the equation is shown in the top left of each figure. Disk labels: AA Tau = 1, CQ Tau = 2, DM Tau = 3, HD 143006 = 5, HD 34282 = 6, J1615 = 8, J1842 = 9, J1852 = 10, LkCa 15 = 11, MWC 758 = 12, PDS 66 = 13, SY Cha = 14, V4046 Sgr = 15. Disks without retrievable surfaces are excluded.</figcaption></figure>

$\langle z/r \rangle$ error. With this point removed, the PCC becomes $r_p = 0.76$.

We have also explored trends in the stellar mass, average temperature, and outer edge for $^{13}$CO emission. In general, we find $\langle z/r \rangle_{13\text{CO}}$ decreases with mass and mean brightness temperature, consistent with C. J. Law et al. (2023). Similar to what we found for $^{12}$CO, we retrieve a linear relationship (albeit much steeper) between $\langle z/r \rangle_{13\text{CO}}$ and $R_{\text{CO}}$ compared to previous studies, although this linear relation suffers from a larger amount of scatter.

G. Rosotti et al. (2025) predicts the emission surface height should positively correlate with the gas surface density and hence the disk mass. T. Paneque-Carreño et al. (2025) also recently presented a relation between $\langle z/r \rangle$ and the disk mass via chemical modeling; they find that the disk mass and volatile carbon abundance are primarily responsible for setting the observed $\langle z/r \rangle$. Here, we check whether the observations are in agreement with the predictions. The bottom left panel of Figure 7 shows $\langle z/r \rangle$ versus the disk mass derived from the rotation curves in C. Longarini et al. (2025; left panel), the disk mass derived from dust mass measurements from P. Curone et al. (2025; middle panel), and the disk mass derived using $N_2H^+$ from L. Trapman et al. (2025; rightmost panel), along with literature values shown in gray. As shown in the figure, we find linear relations between $\langle z/r \rangle$ and the disk mass derived

from multiple methods. For $\langle z/r \rangle$ versus the rotation disk mass, we find the strongest linear relation, with a PCC of $r_p = 0.64$. The emission height versus the disk mass via the dust mass has a PCC of $r_p = 0.38$, and $\langle z/r \rangle$ versus the disk mass via $N_2H^+$ has the weakest linear relation with a PCC of $r_p = 0.20$. In the former two plots, HD 34282 (6) and J1615 (8) stand as outliers, respectively. This may hint that the $\langle z/r \rangle$/disk mass relation flattens out at higher masses, or that the relation is better quantified by a nonlinear equation. T. Paneque-Carreño et al. (2025) finds that the $\langle z/r \rangle$/disk mass relation differs between T Tauri and Herbig stars, which may at least partially contribute to HD 34282 being an outlier, although we note that MWC 758 and a few literature sources are also Herbig stars and do not exhibit outlier behavior. The strong correlation between $\langle z/r \rangle$ and the disk mass measured with rotation curves may be in part due to the fact that this method uses the emission surface as a prior. However, independent of the mass measurements, we find a strong linear relation between $\langle z/r \rangle$ and $R_{\text{edge}}$. This implies that the relation between $\langle z/r \rangle$ and the dynamical disk mass is not entirely due to the use of emission surface. The results presented in Figure 7 imply that the disk mass and $\langle z/r \rangle$ are physically related, as predicted in G. Rosotti et al. (2025) and modeled in T. Paneque-Carreño et al. (2025).

While various methods have been proposed to measure the mass and surface density of protoplanetary disks (e.g.,
S. M. Andrews et al. 2024; C. Longarini et al. 2025; G. Rosotti et al. 2025; L. Trapman et al. 2025; see also review by A. Miotello et al. 2023), these methods require observations of multiple molecular lines (L. Trapman et al. 2025) or sophisticated velocity analyses (S. M. Andrews et al. 2024; C. Longarini et al. 2025). Compared with other approaches, the extraction of ¹²CO emission surface is relatively straightforward, offering an alternative way to estimate disk mass and surface density, as in T. Paneque-Carreño et al. (2025). We also calculate $\langle z/r \rangle$ versus the gas surface density $\Sigma$ at 100 au, inferred based on the rotational disk masses assuming the surface density profile used in C. Longarini et al. (2025), where $\gamma = 1$ (see their Equation (A1)). We find a positive relationship between the emission height and surface density, with a PCC of $r_p = 0.75$, in agreement with the prediction by G. Rosotti et al. (2025). To this end, we provide fits between $\langle z/r \rangle$ and the disk mass and surface density, with the caveat that this relation is specific to dynamical masses, and other methods will not follow this relation as closely:

$$ \frac{M_{\text{disk}}}{M_{\odot}} = 0.020^{+0.047}_{-0.056} + 0.26^{+0.23}_{-0.089} \langle z/r \rangle \quad (7) $$

$$ \frac{\Sigma(100 \text{ au})}{\text{g cm}^{-2}} = 1.010^{+2.26}_{-0.29} + 11.41^{+14.19}_{-7.22} \langle z/r \rangle. \quad (8) $$

The remaining linear relations derived using the disk mass via dust and disk mass via N₂H⁺ are shown in the upper right corners of Figure 7. We use emcee (D. Foreman-Mackey et al. 2013) to obtain these linear fits with 100 walkers and 5000 iterations. This allows us to place priors such that the model cannot be below zero. The errors are propagated from $\langle z/r \rangle$ and $M_{\text{disk}}$.

Additionally, we explore the relationship between $\langle z/r \rangle$ and disk mass for ¹³CO. We find similar correlations to that seen with ¹²CO. For $\langle z/r \rangle$ versus the dynamical disk mass, $r_p = 0.88$; for $\langle z/r \rangle$ versus the disk mass via the dust mass, $r_p = 0.27$; for $\langle z/r \rangle$ versus the disk mass via N₂H⁺, $r_p = 0.34$. We note that we do not include archival estimates here.

## 6.2. Comparison to Previous Results

Several disks from the exoALMA sample have previous gas component observations at a different spatial and angular resolution, and with other molecular transitions. In the following sections, we compare our results to previous studies of the same disks when available.

### 6.2.1. Emission Surfaces

Previous work by C. J. Law et al. (2023) derived the emission heights, temperature profiles, and outer radii for DM Tau (¹²CO J = 3–2 and ¹³CO J = 2–1 emission), HD 34282 (¹²CO J = 2–1 and ¹²CO J = 3–2 emission), and LkCa 15 (¹²CO J = 2–1, ¹²CO J = 3–2, ¹³CO J = 2–1, and ¹³CO J = 3–2 emission). Additionally, C. J. Law et al. (2022) derived the emission surface of V4046 Sgr (¹²CO J = 3–2 emission). Below, we compare the emission surfaces between this work and the previous work.

The previous DM Tau emission surfaces for ¹²CO and ¹³CO J = 2–1 were extracted at angular resolutions of ~0.13 and ~0.19, respectively. C. J. Law et al. (2023) finds a $\langle z/r \rangle$ of 0.41$_{-0.17}^{+0.7}$ in ¹²CO J = 2–1 and 0.22$_{-0.04}^{+0.8}$ in ¹³CO J = 2–1. The surface was fit out to 2.7, or 388.8 au. These results differ from ours; we find shallower ¹²CO and ¹³CO surfaces, with a $\langle z/r \rangle$

of 0.30 and 0.13, respectively. When compared to our parametric surfaces, that of C. J. Law et al. (2023) is consistent, albeit slightly higher and with a much shorter cutoff radius. The differences in emission surfaces can be partially explained by the difference in angular and spectral resolution between the observations, as well as the implementation of the iterative function in disksurf, which can lead to slight differences in the observed surface structure. In particular, the iterative surface process neglects the outer diffuse emission seen in DM Tau, and instead finds a tapered power-law surface. We have included the diffuse outer points in Figure 2 and Figure 9 in Appendix A; if we include these in our surface fit, we find results more consistent with that of C. J. Law et al. (2023).

The ¹²CO and ¹³CO emission of LkCa 15 was first mapped by M. Leemker et al. (2022), and C. J. Law et al. (2023) derived the emission height and 2D temperature profile using the J = 2–1 and J = 3–2 transitions of ¹²CO and ¹³CO. C. J. Law et al. (2023) finds $\langle z/r \rangle$ values of 0.26 and 0.13 for the respective J = 2–1 transitions, and $\langle z/r \rangle$ values of 0.23 and 0.15 for the J = 3–2 transitions. In comparison, we find elevated $\langle z/r \rangle$ values of 0.34 and 0.21 for ¹²CO and ¹³CO. When overplotted, our derived surface is consistent with that from C. J. Law et al. (2023), the main difference being their parametric surface is slightly below the one we find. These differences may be due to the resolution (~0.27 versus 0.15), or due to the iterative function in disksurf, which was not implemented in C. J. Law et al. (2023).

The final disk from C. J. Law et al. (2023) that overlaps with the exoALMA sample is HD 34282. They find $\langle z/r \rangle$ values of 0.46 and 0.38 for the ¹²CO J = 2–1 and J = 3–2 transitions. In comparison, we find a $\langle z/r \rangle$ value of 0.41, consistent with the 0.38 value Law et al. finds. When comparing emission surface morphology, our surface more closely matches that of their derived J = 3–2 surface, as expected, although we find our surface has a higher peak and a sharper taper.

The V4046 Sgr disk emission surface and properties were previously analyzed in C. J. Law et al. (2022). Our emission surface differs significantly from the previously derived one; C. J. Law et al. (2022) finds a maximum outer radius of ~160 au, while our surface extends out to ~400 au. Their reported outer radius of ~160 au matches closely with the large gap observed at 173 au in this work; it is likely that the difference in sensitivities allowed us to pick up the emission in the outer portion of this disk past the gap.

In general, our results agree with those previously reported. Some of the biggest differences in emission surface morphology and outer radius are likely due to differences in angular resolution, spectral resolution, and sensitivity.

### 6.2.2. Temperature Structure

Both DM Tau and LkCa 15 have previously derived 2D temperature profiles. For DM Tau, we find mostly consistent temperature profiles to that of C. J. Law et al. (2023). We find an atmospheric temperature of 37 K and a midplane temperature of 20 K; C. J. Law et al. (2023) finds $T_{\text{atm}} = 38$ and $T_{\text{mid}} = 26$ K. The 2D shape of the profile is also mostly consistent, with the edge of the ¹³CO surface having a temperature of ~20 K, and the rising portion of the ¹²CO surface at ~30 K.

For the temperature structure of LkCa 15, we find a hotter atmospheric temperature ($T_{\text{atm}}$) of 48 K, compared with 35 K from C. J. Law et al. (2023). The midplane temperatures are
consistent; we report 20 K, and they report 21 K. The discrepancy between atmospheric temperatures may be because C. J. Law et al. (2023) included all temperature data points, even those below 20 K, whereas we exclude these points; additionally, we use different equations to connect the temperature structure.

In C. J. Law et al. (2023), the temperature profile of HD 34282 in $^{12}$CO $J = 2-1$ and $J = 3-2$ is derived. We find generally consistent radial profile shapes, identifying the same bump as C. J. Law et al. (2023) at around 295 au. Our data extend out another ~300 au due to our achieved sensitivity. However, the overall temperature we measure for this peak intensity profile is considerably higher (by about 15 K), which is likely due to the difference in beam sizes between the data sets.

Any differences in the temperature profiles can be attributed to several factors: the CO transition used is different, the spectral and angular resolution of the cubes is not the same, and the functional form used to fit the 2D temperatures is not the same. Despite these differences, it seems that we have independently come to the same answers, implying that our temperature prescriptions are comparable.

### 6.3. Diffuse Emission

Some disks exhibit diffuse emission at large radii, visible in the channel maps (see C. Pinte et al. 2025), moment maps, and azimuthally averaged emission surfaces. We see a sharp drop-off in AA Tau and SY Cha in $^{12}$CO emission, and a similar morphology in LkCa 15 $^{13}$CO emission. In DM Tau, the outer surface does not taper, and instead extends outwardly as diffuse emission.

One explanation is that these disks may be experiencing UV photodesorption. When a strong enough external radiation field is present, gas may be irradiated in the outer disk, which would allow for CO to emit from regions at and near the midplane, extending out radially. The emission morphology of IM Lup is similar to that of the disks described here. L. I. Cleeves et al. (2016) modeled IM Lup and the surrounding radiation field and found that UV photodesorption can explain the extended emission in IM Lup. This was further confirmed by C. Pinte et al. (2018), who modeled IM Lup and found that a disk experiencing interstellar medium (ISM) radiation and UV chemistry most closely matched the observations. CO gas in the midplane indicative of photodesorption has also been observed in the edge-on Flying Saucer disk (A. Dutrey et al. 2017), which resides in the ρ Ophiuchus star-forming region. Photodesorption can be due to stellar irradiation, or from the surrounding ISM. AA Tau, LkCa 15, and DM Tau are all located in the Taurus star-forming region, which is lacking in massive stars compared to other star-forming regions (K. L. Luhman 2004). This means that Taurus has a lower far-ultraviolet (FUV) radiation field as compared to other star formation rates, such as Lupus; Taurus is generally thought to have no photoionizing sources.

However, the diffuse emission we observe suggest it may be possible that, even in low FUV areas like Taurus, disks still experience some photodesorption. The disks that we observe the diffuse emission of are some of the most massive in the exoALMA sample, with $R_{\text{edge}}$ ranging from 644 au (AA Tau) to 990 au (DM Tau). It may be the case that the outer regions of the disks are low density enough to be photodesorbed, even in an FUV field like Taurus's. Recent work in R. Anania et al.

(2025) indicates that regions with late-type B and early-type A stars such as Taurus can reach nonnegligible levels of FUV radiation, in agreement with the preliminary findings here. However, more work is needed to model these disks to determine if they are truly experiencing photodesorption, or if there are other potential causes (see Section 6.4.1).

### 6.4. Origins and Implications of Substructures in Emission Surfaces and Intensity Profiles

In Section 5, we detailed the substructures identified in the emission surfaces and peak intensities, and compared their radial locations to continuum rings/gaps and velocity kink signatures. Here, we discuss the potential origin of these substructures, as well as their implications.

Many of the emission surface dips are colocated between molecular tracers. This suggests that these features are likely physical and coherent over a range of heights. Several of the surface gaps also coincide within a beam size with peak intensity modulations in AA Tau, J1842, LkCa 15, and SY Cha. Planets forming within the disk are expected to open gaps, whose signatures could be imprinted on the azimuthally averaged emission surfaces. As shown in Figure 12 in Appendix A, surface dips in the disks of AA Tau, J1842, and LkCa 15 coincide with the radial locations of velocity kinks.

A few peak intensity bright spots coincide with proposed planet locations in AA Tau, J1615, LkCa 15, and SY Cha, which could be attributed to processes induced by planet formation. As a planet forms, there are several methods by which it is expected to alter the disk. Forming protoplanets will heat the surrounding gas (L. I. Cleeves et al. 2015; J. Bae et al. 2022). This can happen via accretion heating, as the planet accretes during its formation (J. Szulágyi et al. 2016), or shock heating as the disk material infalls onto the planet and/or circumplanetary disk (CPD; J. Szulágyi & C. Mordasini 2017). Despite this, the planet and subsequent CPD are expected to exist at the midplane, and their ability to cause localized heating in the molecular layers above this is unclear. The radial profiles are an azimuthal average, implying that the bright spots persist throughout the disk and do not just stem from a singular azimuthal location. Additionally, it is expected that a planet will carve a gap in the gas as it forms, which would likely manifest as a temperature deficit in the radial intensity profiles (S. Facchini et al. 2018). Whether or not the signatures identified here are the result of forming planets is yet to be seen. Further efforts to model the impact of planet formation on the gaseous components of thermally stratified disks is needed to characterize these peak intensity signatures.

Surface gaps could also be “artificially” induced when we make the assumption that the disk is rotating at the Keplerian velocity. If a disk is experiencing velocity perturbations, this could cause artifacts on the surface when the isovelocity curve is drawn (see also T. Paneque-Carreño et al. 2023, who discusses the same possibility). This could be another explanation as to why several of the surface gaps coincide with the proposed planet locations. However, even if the signatures are induced by non-Keplerian velocity deviations, their characterization is still helpful in pinpointing potential sites of planet formation.

LkCa 15. Besides the peak intensity bright spot coinciding with the proposed planet location, the LkCa 15 disk exhibits “wiggles” in the $^{12}$CO $J = 3-2$ peak intensity profile. The
intensity peaks show a striking wavelike pattern that extends radially, with the bright spots being placed every ~1751 au. Additionally, one of the peak intensity bright spots at 333 au coincides within a beam size with an emission surface gap at 342 au, and both of these signatures nearly coincide with the proposed planet location at 370 au (see Figures 12, 14 in Appendices A and C, respectively). A. F. Izquierdo et al. (2025) explores the 2D residual maps of several exoALMA sources; we find consistent radial peak intensity profiles, and the residual map of LkCa 15 clearly displays the peaks and bumps that persist azimuthally throughout the disk. In this work, we also see evenly spaced peak intensity “wiggles” to a lesser extent in AA Tau and SY Cha, but we focus on the clearest case of LkCa 15.

LkCa 15 is thought to have several planets embedded within its disk. First, there are several rings and a larger cavity in continuum emission (S. M. Andrews et al. 2011; S. Facchini et al. 2020; P. Curone et al. 2025), with dust trapping along two Lagrangian points (F. Long et al. 2022b), all indicative of potential ongoing planet formation processes impacting the disk structure. Several studies within the exoALMA collaboration find further evidence of ongoing planet formation in LkCa 15. The disk exhibits a high degree of $\delta v_\phi$ structure, derived in J. Stadler et al. (2025), implicative of pressure variations that could be linked to planet formation. C. H. Gardner et al. (2025) models the continuum and gaseous disk components, concluding that the observed signatures are consistent with massive planets orbiting between 10 and 30 au. Although there are no direct detections of any planet in the disk, it remains likely that LkCa 15 has at least one embedded protoplanet. This raises the possibility that the peak intensity wiggles are linked to the planet formation process.

One theory is that the signatures seen in the LkCa 15 peak intensity profiles are the result of tightly wound spirals. Several exoALMA sources display large-scale spirals (e.g., HD 135344B, MWC 758; A. F. Izquierdo et al. 2025; C. Pinte et al. 2025). While those sources have more obvious spirals, it is possible that the peak/bump structure in LkCa 15 actually stems from spirals, which could explain their even spacing. Simulations have established that planets can launch spiral arms throughout the disk (W. Kley 1999; R. Dong et al. 2015; J. Bae et al. 2016), through Linblad resonances (G. I. Ogilvie & S. H. Lubow 2002; J. Bae & Z. Zhu 2018a, 2018b) or buoyancy resonances (Z. Zhu et al. 2012; J. Bae et al. 2021). Observationally, temperature spirals have been characterized in TW Hya and MWC 480 (R. Teague et al. 2019b; R. Teague et al. 2021b; L. Wölfer et al. 2021; A. F. Izquierdo et al. 2023), where the authors postulate that a giant planet may have induced the signatures. However, if the temperature structure in LkCa 15 is driven by buoyancy spirals, we would expect the velocity residuals to be dominated by vertical gas motions, $v_z$. Based on the residuals in A. F. Izquierdo et al. (2025); see their Figure 10), this may not be the case in LkCa 15. Because the peak intensity structure aligns with some of the structure seen in $v_\phi$ (J. Stadler et al. 2025), it is more likely that these peaks are driven by pressure variations. More work is needed to determine if the disk structure is truly a spiral, or rather concentric peaks and troughs. If a planet is driving the temperature variations, it likely exists at large radii; modeling the impacts of giant planet formation on the outer gaseous disk remains a topic of interest.

Thermal wave instability, or irradiation instability, could also produce wavelike features. The central star is responsible for much of the disk heating; this is related to the disk thickness and orientation. If parts of the disk receive more thermal radiation, this can lead to a growing perturbation that will produce peaks and troughs in the surface density (P. D’Alessio et al. 1999; R. Siebenmorgen & F. Heymann 2012; T. Ueda et al. 2021; Y. Wu & Y. Lithwick 2021). However, these studies have generally focused on the inner disk; this instability is known to operate at radii closer to the star because of the thermal requirement. The wave structure in LkCa 15 extends out to 900 au, which may be too far for irradiation instability to effectively operate.

Of course, other physical mechanisms are expected to cause disk perturbations as well. Magnetorotational instability (MRI), which generates turbulence, can excite density perturbations, which could manifest as spirals (M. Flock et al. 2015). However, MRI is not expected to produce symmetric features, which makes it an unlikely candidate. Vertical shear instability (VSI; R. P. Nelson et al. 2013) can also induce surface density perturbations. Due to angular velocity, the surface layers of protoplanetary disks rotate more slowly than the midplane, which naturally induces VSI. Numerical simulations have shown that VSI can produce antisymmetric structures. M. Barraza-Alfaro et al. (2025) and M. Barraza-Alfaro et al. (2021) find that VSI can manifest as concentric rings in the velocity residuals (i.e., the “corrugation mode”), but these signatures are not seen in the modeled peak intensity maps. Concentric patterns in the surface density and velocity structure are also retrieved in simulations by S. Zhang et al. (2024), who investigates VSI in a thermally stratified disk. They find that perturbations in the gas surface density can develop, sometimes exhibiting wavelike structure. However, due to the high inclination of LkCa 15, it is difficult to definitively say that the signatures are driven by VSI. It is likely that multiple physical processes are impacting the disk, and further investigation and modeling work are needed.

### 6.4.1. *Projection Effects*

Throughout Section 5, we have alluded to the fact that some of the observed substructures can be explained by projection effects, and may not be physical. Here, we discuss the ways in which our viewing angle, optical depth effects, and temperature disparities could cause artificially induced dips.

*Diffuse emission.* First, the “diffuse” emission detected in AA Tau, LkCa 15, and SY Cha all have a characteristic shape, seen in Figures 9 and 10 in Appendix A. All of these disks are highly inclined (i > 50°), and all of the surfaces suffer an abrupt drop, at which point the more diffuse emission from closer to the midplane continues. In the channel maps, this emission does appear as a combination of emission from the front, back, and midplane, which implies that CO is not frozen at the midplane and indicative of photodesorption. However, at high inclinations, from our viewing angle, the back surface will appear to extend beyond the front surface, and in some cases, can appear just as bright. With real observations, this projection effect could mimic diffuse emission, particularly in the presence of noise or non-Keplerian rotation. One hint that this effect may be at play is that we only see a sharp surface drop in LkCa 15 $^{13}$CO emission, not $^{12}$CO. If the morphology is due to photodesorption, it is unlikely that it would only impact $^{13}$CO.
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Two channel maps of V4046 Sgr at 2.00 and 3.70 km s⁻¹. The left column shows two channels with no points overlaid, and the right column shows the channels with the raw emission surface points identified by *disksurf*. There is a notable gap in the emission surface points that is not visible in the channels.</figcaption></figure>

V4046 Sgr. The V4046 Sgr disk exhibits one of the most prominent emission surface substructures in both width and depth out of the exoALMA sample in both ¹²CO and ¹³CO at 174 au, where a large decrease in surface height can be seen in Figure 2. Curiously, the channel maps and peak intensity map of this disk appear relatively smooth, and no gap can be identified by eye. The intensity does appear to decrease more sharply past ~150 au (see Figure 4), but other disks that have similar intensity profile slope changes do not show large gaps. C. Pinte et al. (2025) presents the channel maps of V4046 Sgr, and notes that there are no structures visible by eye. Figure 8 presents two channels in ¹²CO; the left column shows the channels without anything overlaid, and the right column shows the channels with the *disksurf* emission surface points for the respective channels. There is no gap visible by eye in the left columns, but the emission surface points do have a gap at the same radial location in multiple channels. It is uncommon that such a large surface gap identified via *disksurf* is not visible in the channel maps.

It is evident from Figure 8 that the intensity of the disk drops off dramatically at the location of the gap. This intensity drop can also be seen in the peak intensity at around 150 au, shown in Figure 4. The gap could be related to the intensity drop; decreases in column density or optical depth could result in a decreased intensity. However, other disks that show the same level of intensity change do not exhibit large gaps as V4046 Sgr does. Additionally, although the intensity decreases, in ¹²CO, it remains above 20 K, which means the gas likely remains moderately optically thick.

The gap in the V4046 Sgr surface may be caused by projection effects. Without a visible gap in the channel or moment maps, it is unlikely that one actually exists, although there are other possibilities. *disksurf* works by identifying peak points along an isovelocity curve (C. Pinte et al. 2018); in general, this process can accurately identify real gaseous emission gaps (C. J. Law et al. 2021a, 2022, 2023; M. Galloway-Sprietsma et al. 2023; T. Paneque-Carreño et al. 2023). In

the case of V4046 Sgr, we hypothesize that this gap is actually caused by low optical depth (e.g., $\tau \simeq 1$) or low temperature of the front surface. If the optical depth and or temperature of the CO emission begins to fall off, this would occur at different radial locations from our viewing angle, due to the inclined angle of the disk. At a certain radial location where the front and back surface are best separated, the back surface would appear slightly brighter than the front. The fact that we only see the gap in the upper channels, or the nearside of disk, lends more credence to this theory; the opposite effect would be at play for the farside of the disk, and we would not see this effect. Inspecting the isovelocity contours from DISCMINER models reveals that the gap falls over the predicted back surface. Therefore, we conclude that this gap is likely not physically real, but rather an effect caused by the projection of the disk, although other explanations remain possible.

We believe that the gap in V4046 Sgr is the only obvious case of projection effects. The other surface modulations, discussed in Section 5, are found in disks that are not smoothly varying. Many of the substructures can be identified visually on the channel maps, and are of a true physical origin. For future high-resolution imaging campaigns, it will be critical to identify and understand the role that projection plays on our observations.

## 7. Conclusions

This work presents the emission surfaces, 2D temperature profiles, and radial intensity profiles of the exoALMA sources. Our analysis leads to the following conclusions:

1. Nearly every disk in the sample exhibits stratified molecular layers, with ¹²CO $J=3-2$ tracing elevated regions of the disk, with an average $\langle z/r \rangle$ of ≈ 0.28. ¹³CO $J=3-2$ traces layers below that of ¹²CO, with an average $\langle z/r \rangle$ of ≈ 0.16. CS $J=7-6$ lies closest to the midplane, but still exhibits elevated layers in several disks, with an average $\langle z/r \rangle$ of ≈0.18.

2. We search for “substructure,” i.e., gaps and dips, in the emission surfaces and radial profiles. Most disks exhibit substructure in both the emission surfaces and temperature profiles, often coincident between molecules and in some cases coinciding with proposed planet locations.

a. In AA Tau, DM Tau, LkCa 15, and SY Cha, we see evidence of diffuse emission at large radii, which generally manifest in the emission surfaces as a sharp drop with subsequent midplane emission. This structure may be due to photodesorption, which implies that this effect may be common even in low FUV star-forming regions.

b. In CQ Tau, DM Tau, J1842, J1852, LkCa 15, MWC 758, and V4046 Sgr, we find a correlation of emission surface substructures between separate molecular layers, potentially indicative of physical processes, which impact the vertical height of the disk at these locations.

c. In AA Tau, DM Tau, J1615, J1842, LkCa 15, MWC 758, and SY Cha, we find a correlation between peak intensity substructures and emission surface substructures.

d. We do not find a statistically significant correlation between the continuum rings and gaps and any surface
substructure. The notable exceptions are AA Tau, J1615, and MWC 758. In the CS surfaces of AA Tau and J1615, we see hints of surface height changes at the location of gaps in the continuum. In MWC 758, the continuum rings and gaps overlap with those identified in the ¹²CO and ¹³CO surfaces.

e. For some disks, structures in the emission surfaces and peak intensities overlap with velocity perturbations seen in the channel maps, discussed in C. Pinte et al. (2025). Whether the surface substructures are truly physical decreases of column density or intensity within the disk, or are instead an “artificial” detection driven by the velocity perturbations themselves, is still unclear. Future work to characterize the origin of these substructures is needed.

f. Projection effects due to inclination and viewing angle may cause artificial structures in our observations. More work is needed to characterize the interplay between disk position and identified substructures.

3. The 2D temperature profiles derived suggest the average $T_{\text{atm}}$ temperature is ~43 K, whereas the average $T_{\text{mid}}$ temperature is found to be ~24 K. This midplane temperature is hotter than that predicted by the scale height equation and simulations, and suggests a need for an improved temperature profile prescription and methods to better constrain these temperatures.

4. We compare our work to previous studies and find consistent results. The stellar mass and average ¹²CO temperature is loosely correlated with ¹²CO surface height, but no definite relations can be identified. Surface height and the outer disk radius $R_{\text{CO}}$ show a strong linear relation as seen in C. J. Law et al. (2021a, 2022, 2023) and T. Paneque-Carreño et al. (2023). Additional linear relations are found between $\langle z/r \rangle$, disk mass, and surface density, suggesting that the emission height is inherently related to these quantities. We propose that $\langle z/r \rangle$ could be used to infer the disk mass.

5. The data discussed in this Letter are provided as a publicly available value-added data product (VADP). The nonparametric emission surfaces (raw, binned, and rolling points, along with their temperatures), parametric surfaces, peak intensity radial profiles, integrated intensity radial profiles, and the radial locations of the surface and intensity modulations are all included in the VADP release.

## Acknowledgments

This Letter 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 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. Support for AFI 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. J.S. has 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 have been done on the "Mesocentre SIGAMM" machine, hosted by Observatoire de la Cote d'Azur. 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 (DUSTBUS-TERS) and by the UK Science and Technology research Council (STFC) via the consolidated grant ST/W000997/1. 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. M.B. has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (PROTOPLANETS, grant agreement No. 101002188). S.F. is funded by the European Union (ERC, UNVEIL, 101076613). Views and opinions expressed 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. S.F. acknowledges financial contribution from PRIN-MUR 2022YP5ACE. 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). Views and opinions expressed are however those of the author(s) only, and do not necessarily reflect those of the European Union or the European Research Council Executive Agency. Neither the European Union nor the granting authority can be held responsible for them. Support for B.Z. was provided by The Brinson Foundation. C.P. acknowledges Australian Research Council funding via FT170100040, DP18010423, DP220103767, and DP240103290. D.F. has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (PROTOPLANETS, grant agreement No. 101002188).

N.C. has received funding from the European Research Council (ERC) under the European Union Horizon Europe research and innovation program (grant agreement No. 101042275, project Stellar-MADE). 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. M.F. is supported by a Grant-in-Aid from the Japan Society for the Promotion of Science (KAKENHI: No. JP22H01274). M.F. is supported by a Grant-in-Aid from the Japan Society for the Promotion of Science (KAKENHI: No. JP22H01274). C.H.G. acknowledges support from the National Aeronautics and Space Administration under grant No. 80NSSC18K0828. C.H.G. gratefully acknowledges the support by an LANL/CSES Student Fellow project. C.H. acknowledges support from NSF AAG grant No. 2407679. J.D.I. acknowledges support from an STFC Ernest Rutherford
Fellowship (ST/W004119/1) and a University Academic Fellowship from the University of Leeds. 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). F.M. acknowledges funding from the European Research Council (ERC) under the European Union's Horizon Europe research and innovation program (grant agreement No. 101053020, project Dust2Pla-nets). D.P. acknowledges Australian Research Council funding via DP18010423, DP220103767, and DP240103290. 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. 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). T.C.Y. acknowledges support by Grant-in-Aid for JSPS Fellows JP23KJ1008.

## Appendix A
Emission Surfaces

Figures 9, 10, and 11 depict the emission surfaces for ¹²CO *J* = 3–2, ¹³CO *J* = 3–2, and CS *J* = 7–6, respectively. In all three figures, the “raw” emission surface *r* — *z* points are shown in gray. The larger black points show these raw surfaces binned by a quarter of the beam size. The dotted lines show the tapered power-law fit from *disksurf*, fit to the background gray points. The solid lines show the tapered power-law fit from DISCMINER, derived via line profile modeling. Table 1 presents the tapered power-law fits from both codes for all molecules.

Figure 12 shows the rolling surfaces derived from the nonparametric surfaces. These rolling surfaces were used to identify substructures in the emission surfaces as discussed in Section 5. The surface substructures are plotted as colored lines and denoted with a Z followed by the radial location in astronomical units, following the nomenclature from C. J. Law et al. (2021a). The continuum rings and gaps from P. Curone et al. (2025) are also plotted as gray solid and dashed lines.
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Nonparametric (data points) and parametric surfaces (solid and dashed curves) for $^{12}$CO J = 3–2. The larger black points represent radially binned data. We exclude two disks (HD 135344B and J1604) whose surfaces could not be extracted due to a face-on geometry. Background gray lines show scale heights (z/r) from 0.1 to 0.6 with 0.1 increments. The darker blue solid line shows the tapered power-law surface fit from DISCMINER, and the light blue and gray dashed line shows the tapered power-law fit from *disksurf*, which uses the background data points. The gray cross markers for AA Tau and SY Cha show the diffuse outer surface points that we do not include in the parametric fits.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>**Figure 10.** The same as Figure 9, but for $^{13}$CO J = 3–2. The orange solid line shows the tapered power-law surface fit from DISCMINER, and the gray and peach dashed line shows the tapered power-law fit from *disksurf*.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 11.** The same as Figure 9, but for CS J = 7–6. The purple solid line shows the tapered power-law surface fit from DISCMINER, and the gray and lavender dashed line shows the tapered power-law fit from *disksurf*.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>**Figure 12.** Rolling surfaces made from the raw *r-z* *disksurf* emission surface points for ¹²CO *J* = 3–2 (blue), ¹³CO *J* = 3–2 (orange), and CS *J* = 7–6 (pink). Vertical dashed lines denote emission surface substructure, whose identification is discussed in Section 5. We label dips in the emission surfaces with the nomenclature first used in C. J. Law et al. (2021a), where “Z” denotes a dip in the emission surface, followed by the radial location in astronomical units. We also plot the continuum rings (solid gray lines) and gaps (dashed gray lines) identified in P. Curone et al. (2025). Thin background gray lines show scale heights (*z/r*) from 0.1 to 0.6 with 0.1 increments. The beam size of 0.15 is shown in the upper right corner for each disk. AA Tau is zoomed in to better show structure.</figcaption></figure>
## Appendix B
Temperature Profiles

Figure 13 shows the nonparametric surface points used to calculate the 2D temperature profiles. The corresponding color of each point is calculated using $\frac{T_{\text{obs}} - T_{\text{fit}}}{T_{\text{obs}}}$, where $T_{\text{obs}}$ is the observed brightness temperature of the point, and $T_{\text{fit}}$ is the corresponding value that the 2D temperature fit finds. Values closer to zero indicate the 2D fit more closely matches the observed temperature structures. Colors closer to blue imply that the temperature fit is underestimating the observed temperature, and colors closer to red imply that the temperature fit is overestimating the observed temperature. It is evident that the 2D temperature profiles cannot account for areas where molecular layers overlap; these areas are visible particularly

where the top of the ¹³CO emission layer overlaps with the bottom of the ¹²CO layer. Additionally, many of the points close to the midplane are overestimated in the 2D temperature fits.

These results highlight a need for a better analytical temperature profile. Neither of most commonly used temperature profile equations, i.e., the sinusoidal form from E. Dartois et al. (2003), and the tangential form from C. P. Dullemond et al. (2020), cannot adequately describe the temperature change between overlapping molecular layers. Additionally, both functional forms tend to be highly sensitive to the imposed priors, which means much of the results are imparted from our initial assumptions of the atmospheric and midplane temperature functions.

<figure><img src="image_15.png" /><figcaption>**Figure 13.** Nonparametric surface points used to calculate the 2D temperature in *r* — *z* space. The color of the points corresponds to the difference between the observed brightness temperature (*T*<sub>obs</sub>) and the temperature found by the 2D fit (*T*<sub>fit</sub>). The boxed numbers in the upper right corners represent the absolute value of the average $\frac{T_{\text{obs}} - T_{\text{fit}}}{T_{\text{obs}}}$; values closer to zero indicate the 2D fit more closely matches the observed temperature structures.</figcaption></figure>
## Appendix C
### Radial Profiles

Figure 14 shows the sources with identified peak intensity substructures. Table 4 lists the substructures identified in the emission surfaces and peak intensity profiles. Figures 15, 16, and 17 show the integrated intensity in Jy beam⁻¹ km s⁻¹, calculated using the fiducial image cubes with GoFish. We leave an analysis of the integrated intensity profiles for future work. Figure 18 compares the peak intensities retrieved from GoFish (black lines), disksurf (black points), and DISCMINER (red lines). All peak intensities in this figure were derived using the fiducial images (see R. Teague et al. 2025 for image details), but the DISCMINER intensities were calculated using the continuum-subtracted cubes. The three codes used to derive the peak intensities are in general agreement, and nearly all are within error bars. Most of the intensities diverge in the inner disk, where beam-dilution causes the observed intensity

to drop. C. J. Law et al. (2021a) hypothesized that the decrease in intensity could be due to the lines becoming optically thin, unresolved substructures, or emission being absorbed by dust. A combination of these effects is a possibility for each disk, but we save an in-depth diagnosis of these individual disks for future work. A few of the disks presented in Figure 14 exhibit non-Kelperian motion, as discussed in C. Pinte et al. (2025), A. F. Izquierdo et al. (2025), and J. Stadler et al. (2025). Nearly all non-Keplerian motion can impact the disk temperature structure; likewise, any temperature variation can induce some non-Keplerian motion. Because of this, we do not definitively associate any of the peak intensity substructures with the observed non-Keplerian motion. Additionally, because these radial profiles are an azimuthal average, only large-scale peak intensity perturbations will be detected. We refer the reader to future exoALMA Letters for a discussion of the signatures of velocity perturbations on peak intensity maps.

<figure><img src="image_16.png" /><figcaption>**Figure 14.** Peak intensities derived from GoFish for disks with identified intensity structure. Intensity bright spots are marked with solid lines and labeled with a “B,” and intensity dips are marked with dashed lines and labeled with a “D.” Several bright spots and gaps are unlabeled due to overlap; we refer the reader to Table 4 for the full list of radial locations. Blue corresponds to ¹²CO, orange corresponds to ¹³CO, and pink corresponds to CS. Proposed planet locations, discussed in C. Pinte et al. (2025), are shown at their suggested radial locations in green. Continuum rings and gaps, identified and discussed in P. Curone et al. (2025), are shown as gray solid lines (rings) and dashed lines (gaps).</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>**Figure 15.** Integrated intensity profiles for $^{12}$CO J = 3–2. The shaded regions represent the upper and lower uncertainties. The beam size of 0.15 is shown in the top right corner.</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>Figure 16. The same as Figure 15, but for $^{13}$CO J = 3–2.</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>Figure 17. The same as Figure 15, but for CS J = 7–6.</figcaption></figure>
<figure><img src="image_20.png" /><figcaption>**Figure 18.** Peak intensities derived from GoFish (black lines), disksurf (white circles), and DISCMINER (red lines). The beam size of 0.15" is shown in the upper right of each plot.</figcaption></figure>
## Appendix D
### Comparing Disksurf to DISCMinER

The exoALMA collaboration utilized two open-source Python codes for the analysis of the emission surfaces: `disksurf` (R. Teague et al. 2021a) and `DISCMinER` (A. F. Izquierdo et al. 2021, 2023). As discussed in Section 3, the emission surfaces retrieved between the two are generally consistent, but there are notable differences (see Figures 9, 10, and 11). Below, we discuss and demonstrate the differences between `disksurf` and `DISCMinER` on simulated and real data.

Figure 19 presents the results of this test on two test cubes with varying inclination and molecular line emission. The left panels shows the results from the 30° inclined simulated disk, and the right panel shows the same results from the same procedure, but on a 60° inclined simulated disk. The top two plots are for ¹²CO J = 3–2 emission, and the bottom two plots show the results for ¹³CO J = 3–2 emission. The red line in both of the plots depicts the τ = 2/3 surface at an inclination of zero; therefore, the τ = 2/3 surface in these plots represents more of an upper limit. The simulated cubes were made with MCFOST (C. Pinte et al. 2006, 2009) and have identical parameters apart from the inclination and molecular composition. It is immediately evident that the inclination has an impact on the derived surface for both codes. For the 30° inclined disk, the `disksurf` surface aligns more closely with the τ = 2/3 surface; for the 60° inclined disk, the opposite is true, with the `DISCMinER` surface falling closely in line with the τ = 2/3 surface. The `disksurf` surface scatter in the z-direction also changes between the two cubes, as discussed in Section 3.3. More highly inclined disks will have a larger surface scatter when using the method implemented in `disksurf`, due to the

fact that the front and back emission layers are less separated, causing some points that belong below the front surface to be picked up.

It is critical to interpret these results from simulated cubes carefully. For the 30° inclined disk, it appears that `disksurf` does a better job at finding the τ = 2/3 surface, and vice versa for the 60° inclined disk. However, the parametric surfaces for the real data, shown in Figures 9, 10, and 11, do not show large differences between `disksurf` and `DISCMinER` derived tapered power laws. Some of the largest differences in parametric surfaces are seen in J1615 ¹²CO emission, where the radial location of the tapering between `disksurf` and `DISCMinER` differs by ~130 au, and the vertical height of the emission surface found by `DISCMinER` is systematically lower. The other disks that show the largest differences between the parametric surfaces are those that suffer from midplane, backside, or diffuse contamination, namely, AA Tau, DM Tau, and SY Cha. Finally, disks with mostly flat surfaces, like those of MWC 758 (¹³CO), PDS 66 (¹³CO, CS), and V4046 Sgr (CS), also suffer from some differences between the two codes. However, unlike the results from the MCFOST cubes, nearly all of the `DISCMinER` parametric surfaces align with the `disksurf` surfaces and are within the area shown by the raw r-z points for the real data. This may indicate that the two codes function differently with simulated cubes, but have a similar performance on real data.

The root cause of the underlying differences between the surfaces, both for the real data and simulated cubes, is still unknown. It is likely that there are several factors at play, and disentangling their individual impact on the retrieved surfaces is beyond the scope of this Letter. As discussed in Section 3.2, there are several factors that may contribute to the differences. First is the innate differences between how `disksurf` and

<figure><img src="image_21.png" /><figcaption>Figure 19. Emission surface comparison between `disksurf` (dashed line and background points) and `DISCMinER` (solid blue line) for two model cubes at 30° and 60°, respectively. The model cubes are identical besides their inclination. The red line denotes the known τ = 2/3 surface. The top two panels show the results for ¹²CO J = 3–2, and the bottom two panels show the results for ¹³CO J = 3–2.</figcaption></figure>
<figure><img src="image_22.png" /><figcaption>**Figure 20.** Velocity residuals calculated with DISCMINER for J1615 (top panels) and LkCa 15 (bottom panels), derived with the *disksurf* parametric emission surface and the DISCMINER parametric emission surface. We subtract the velocity residuals (labeled “Residual”), and subtract the absolute value of the velocity residuals (labeled “Absolute Residual”).</figcaption></figure>

DISCMINER work; *disksurf* utilizes the method outlined in C. Pinte et al. (2018) to derive *r* — *z* points, which relies on finding peaks along isovelocity contours. If the front and back surfaces are not well separated, this procedure can easily confuse points belonging to regions below the emission surface as part of the upper surface, which can lead to artificially lower surfaces. Additionally, since *disksurf* is reliant on intensity peaks, the shape of the line profile can induce biases. On the other hand, DISCMINER uses a line-profile analysis to distinguish whether emission belongs to the front or back surface (see Figure 3 from A. F. Izquierdo et al. 2021), so the surfaces derived by DISCMINER should not be as contaminated with midplane and backside emission. Another potential cause for differences is the velocity; *disksurf* assumes purely Keplerian rotation, whereas DISCMINER takes into account any velocity deviations during modeling. This could potentially explain some of the differences seen in the surfaces, but even so, J1615, which has one of the biggest differences in parametric surfaces, does not have velocity deviations that are significantly larger than those of other disks (see J. Stadler et al. 2025).

Figure 20 shows velocity residuals for J1615 (top panels) and LkCa 15 (bottom panels) calculated with DISCMINER, one using the *disksurf* parametric emission surface and the other using the DISCMINER surface. We take two residuals to compare the resulting differences; the first residual is simply a subtraction between residuals, and the second residual is a subtraction of the absolute values of the residuals. The differences between the results can be seen in several areas of the disks. The J1615 emission surfaces differ significantly (as seen in Figure 9), and this effect can be seen in the velocity residuals. The *disksurf* surface tapers more sharply, resulting in stronger velocity residuals in the outer portions of the disk. Additionally, the *disksurf* surface is shifted compared to that of DISCMINER, resulting in overall stronger residuals throughout the disk. The strongest residuals can be seen in the farside of the disk. The J1615 surfaces are some of the most discrepant between *disksurf* and DISCMINER; clearly, the velocity residuals are sensitive to the chosen emission surfaces. For LkCa 15, there are slight differences between the *disksurf* and DISCMINER surfaces, but not as large as those seen in J1615. Subsequently, the differences in
the residuals of LkCa 15 are not as severe. The sharper tapering of the *disksurf* surface again leads to larger residuals toward the farthest minor axis of the disk.

The question of which surface to utilize will depend on the individual applications. The *disksurf* nonparametric surfaces allow for the characterization of substructures within the surfaces such as bumps and dips; these can be used in chemical models to account for areas of lowered column density, and are helpful visual indicators for gaseous emission structures that may not have been seen on a visual inspection of the data cube. Additionally, the *disksurf* surfaces can be more readily used for the calculation of surface errors because they provide scatter in the z-direction, which is more representative of physical reality. The DISCMINER surfaces are useful when parametric surfaces are needed, particularly for a velocity residual analysis. Because DISCMINER models multiple components of the disk at once, including the emission surfaces, using these parametric surfaces derived during the modeling is the most consistent method when utilizing other DISCMINER outputs. If the *disksurf* surface is artificially imposed in DISCMINER, this can lead to nonphysical residuals, because DISCMINER relies on the surface modeling to minimize these residuals.

## ORCID iDs

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

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

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

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

Cristiano Longarini https://orcid.org/0000-0003-4663-0318

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

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

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

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

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

Giovanni Rosotti https://orcid.org/0000-0003-4853-5736

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

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

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

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

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

Nicolás Cuello https://orcid.org/0000-0003-3713-8073

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

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

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

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

Charles H. Gardner https://orcid.org/0009-0003-8984-2094

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

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

François Menard https://orcid.org/0000-0002-1637-7393

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

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

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

David J. Wilner https://orcid.org/0000-0003-1526-7587

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

## References

Aikawa, Y., & Herbst, E. 1999, A&A, 351, 233

Anania, R., Winter, A. J., Rosotti, G., et al. 2025, arXiv:2501.18752

Andrews, S. M., Huang, J., Pérez, L. M., et al. 2018, ApJL, 869, L41

Andrews, S. M., Rosenfeld, K. A., Wilner, D. J., & Bremer, M. 2011, ApJL, 742, L5

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

Ansdell, M., Williams, J. P., Trapman, L., et al. 2018, ApJ, 859, 21

Bae, J., Nelson, R. P., Hartmann, L., & Richard, S. 2016, ApJ, 829, 13

Bae, J., Teague, R., Andrews, S. M., et al. 2022, ApJL, 934, L20

Bae, J., Teague, R., & Zhu, Z. 2021, ApJ, 912, 56

Bae, J., & Zhu, Z. 2018a, ApJ, 859, 118

Bae, J., & Zhu, Z. 2018b, ApJ, 859, 119

Barraza-Alfaro, M., Flock, M., Béthune, W., et al. 2025, ApJL, 984, L21

Barraza-Alfaro, M., Flock, M., Marino, S., & Pérez, S. 2021, A&A, 653, A113

Barrière-Fouchet, L., Gonzalez, J. F., Murray, J. R., Humble, R. J., & Maddison, S. T. 2005, A&A, 443, 185

Calahan, J. K., Bergin, E. A., Zhang, K., et al. 2021, ApJS, 257, 17

Cleeves, L. I., Bergin, E. A., & Harries, T. J. 2015, ApJ, 807, 2

Cleeves, L. I., Öberg, K. I., Wilner, D. J., et al. 2016, ApJ, 832, 110

Curone, P., Facchini, S., Andrews, S. M., et al. 2025, ApJL, 984, L9

D'Alessio, P., Calvet, N., & Hartmann, L. 2001, ApJ, 553, 321

D'Alessio, P., Cantó, J., Hartmann, L., Calvet, N., & Lizano, S. 1999, ApJ, 511, 896

Dartois, E., Dutrey, A., & Guilloteau, S. 2003, A&A, 399, 773

Dong, R., Zhu, Z., & Whitney, B. 2015, ApJ, 809, 93

Duchêne, G., Ménard, F., Stapelfeldt, K. R., et al. 2024, AJ, 167, 77

Dullemond, C. P., Isella, A., Andrews, S. M., Skobleva, I., & Dzyurkevich, N. 2020, A&A, 633, A137

Dutrey, A., Guilloteau, S., Piétu, V., et al. 2017, A&A, 607, A130

Facchini, S., Benisty, M., Bae, J., et al. 2020, A&A, 639, A121

Facchini, S., Pinilla, P., van Dishoeck, E. F., & de Juan Ovelar, M. 2018, A&A, 612, A104

Flock, M., Ruge, J. P., Dzyurkevich, N., et al. 2015, A&A, 574, A68

Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306

Galloway-Sprietsma, M., Bae, J., Teague, R., et al. 2023, ApJ, 950, 147

Gardner, C. H., Isella, A., Li, H., et al. 2025, ApJL, 984, L16

Garrod, R. T., & Herbst, E. 2006, A&A, 457, 927

Hernández-Vera, C., Guzmán, V. V., Artur de la Villarmois, E., et al. 2024, ApJ, 967, 68

Izquierdo, A. F., Stadler, J., Bae, J., 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

Kley, W. 1999, arXiv:astro-ph/9909394

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., Loomis, R. A., Teague, R., et al. 2021b, ApJS, 257, 3

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

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

Leemker, M., Booth, A. S., van Dishoeck, E. F., et al. 2022, A&A, 663, A23

Lodato, G., Rampinelli, L., Viscardi, E., et al. 2023, MNRAS, 518, 4481

Long, F., Andrews, S. M., Rosotti, G., et al. 2022a, ApJ, 931, 6

Long, F., Andrews, S. M., Zhang, S., et al. 2022b, ApJL, 937, L1

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

Loomis, R. A., Facchini, S., Benisty, M., et al. 2025, ApJL, 984, L7

Luhman, K. L. 2004, ApJ, 617, 1216

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

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
Nelson, R. P., Gressel, O., & Umurhan, O. M. 2013, MNRAS, 435, 2610

Öberg, K. I., Facchini, S., & Anderson, D. E. 2023, ARA&A, 61, 287

Ogilvie, G. I., & Lubow, S. H. 2002, MNRAS, 330, 950

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. 2022, A&A,
666, A168

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

Pinte, C., Harries, T. J., Min, M., et al. 2009, A&A, 498, 967

Pinte, C., Ilee, J., Huang, J., et al. 2025, ApJL, 984, L15

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

Pinte, C., Ménard, F., Duchêne, G., & Bastien, P. 2006, A&A, 459, 797

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

Rich, E. A., Teague, R., Monnier, J. D., et al. 2021, ApJ, 913, 138

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

Rosotti, G., Longarini, C., Paneque-Carreño, T., et al. 2025, ApJL, 984, L20

Sanchis, E., Testi, L., Natta, A., et al. 2020, A&A, 633, A114

Siebenmorgen, R., & Heymann, F. 2012, A&A, 539, A20

Stadler, J., Benisty, M., Winter, A. J., et al. 2025, ApJL, 984, L11

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

Szulágyi, J., Masset, F., Lega, E., et al. 2016, MNRAS, 460, 2853

Szulágyi, J., & Mordasini, C. 2017, MNRAS, 465, L64

Teague, R. 2019, JOSS, 4, 1632

Teague, R., Bae, J., Aikawa, Y., et al. 2021b, ApJS, 257, 18

Teague, R., Bae, J., & Bergin, E. A. 2019a, Natur, 574, 378

Teague, R., Bae, J., Huang, J., & Bergin, E. A. 2019b, ApJL, 884, L56

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

Teague, R., & Foreman-Mackey, D. 2018, bettermoments: A Robust Method
to Measure Line Centroids, v1.0, Zenodo, doi:10.5281/zenodo.1419754

Teague, R., Guilloteau, S., Semenov, D., et al. 2016, A&A, 592, A49

Teague, R., Law, C., Huang, J., & Meng, F. 2021a, JOSS, 6, 3827

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

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

Ueda, T., Flock, M., & Birnstiel, T. 2021, ApJL, 914, L38

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

van der Plas, G., Casassus, S., Ménard, F., et al. 2014, ApJL, 792, L25

van Zadelhoff, G. J., van Dishoeck, E. F., Thi, W. F., & Blake, G. A. 2001, A&A, 377, 566

Villenave, M., Ménard, F., Dent, W. R. F., et al. 2020, A&A, 642, A164

Villenave, M., Stapelfeldt, K. R., Duchêne, G., et al. 2024, ApJ, 961, 95

Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, NatMe, 17, 261

Walsh, C., Nomura, H., & van Dishoeck, E. 2015, A&A, 582, A88

Weaver, E., Isella, A., & Boehler, Y. 2018, ApJ, 853, 113

Wölfer, L., Facchini, S., Kurtovic, N. T., et al. 2021, A&A, 648, A19

Wu, Y., & Lithwick, Y. 2021, ApJ, 923, 123

Yen, H.-W., Koch, P. M., Liu, H. B., et al. 2016, ApJ, 832, 204

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

Zhang, S., Zhu, Z., & Jiang, Y.-F. 2024, ApJ, 968, 29

Zhu, Z., Stone, J. M., & Rafikov, R. R. 2012, ApJL, 758, L42