exoALMA. XIV. Gas Surface Densities in the RX
J1604.3-2130 A Disk from Pressure-broadened CO Line
Wings

Tomohiro C. Yoshida, Pietro Curone, Jochen Stadler, Stefano Facchini,
Richard Teague, Munetake Momose, Sean M. Andrews, Jaehan Bae, Marcelo
Barraza-Alfaro, Myriam Benisty, et al.

► To cite this version:

Tomohiro C. Yoshida, Pietro Curone, Jochen Stadler, Stefano Facchini, Richard Teague, et al. exoALMA. XIV. Gas Surface Densities in the RX J1604.3-2130 A Disk from Pressure-broadened CO Line Wings. The Astrophysical Journal Letters, 2025, 984, 10.3847/2041-8213/adc42f. insu-05099873

HAL Id: insu-05099873

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

Submitted on 6 Jun 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

exoALMA. XIV. Gas Surface Densities in the RX J1604.3–2130 A Disk from Pressure-broadened CO Line Wings

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

<sup>1</sup> National Astronomical Observatory of Japan, 2-21-1 Osawa, Mitaka, Tokyo 181-8588, Japan; tomohiroyoshida.astro@gmail.com

<sup>2</sup> Department of Astronomical Science, The Graduate University for Advanced Studies, SOKENDAI, 2-21-1 Osawa, Mitaka, Tokyo 181-8588, Japan
<sup>3</sup> Dipartimento di Fisica, Università degli Studi di Milano, Via Celoria 16, 20133 Milano, Italy

<sup>4</sup> Departamento de Astronomía, Universidad de Chile, Camino El Observatorio 1515, Las Condes, Santiago, Chile

<sup>5</sup> Laboratoire Lagrange, Université Côte d'Azur, CNRS, Observatoire de la Côte d'Azur, 06304 Nice, France

<sup>6</sup> Department of Earth, Atmospheric, and Planetary Sciences, Massachusetts Institute of Technology, Cambridge, MA 02139, USA

<sup>7</sup> College of Science, Ibaraki University, 2-1-1 Bunkyo, Mito, Ibaraki 310-8512, Japan

<sup>8</sup> Center for Astrophysics | Harvard & Smithsonian, Cambridge, MA 02138, USA

<sup>9</sup> Department of Astronomy, University of Florida, Gainesville, FL 32611, USA

<sup>10</sup> Max-Planck Institute for Astronomy (MPIA), Königstuhl 17, 69117 Heidelberg, Germany

<sup>11</sup> School of Physics and Astronomy, Monash University, Clayton VIC 3800, Australia

<sup>12</sup> Department of Physics and Astronomy, The University of Georgia, Athens, GA 30602, USA

<sup>13</sup> Center for Simulational Physics, The University of Georgia, Athens, GA 30602, USA

<sup>14</sup> Institute for Artificial Intelligence, The University of Georgia, Athens, GA, 30602, USA

<sup>15</sup> Department of Astronomy, Columbia University, 538 W. 120th Street, Pupin Hall, New York, NY, USA

<sup>16</sup> School of Physics and Astronomy, University of Leeds, Leeds, LS2 9JT, UK

<sup>17</sup> Leiden Observatory, Leiden University, P.O. Box 9513, NL-2300 RA Leiden, The Netherlands

<sup>18</sup> European Southern Observatory, Karl-Schwarzschild-Str. 2, D-85748 Garching bei München, Germany

<sup>19</sup> Univ. Grenoble Alpes, CNRS, IPAG, 38000 Grenoble, France

<sup>20</sup> Institute of Astronomy, University of Cambridge, Madingley Road, CB3 0HA, Cambridge, UK

<sup>21</sup> National Radio Astronomy Observatory, 520 Edgemont Rd., Charlottesville, VA 22903, USA

<sup>22</sup> Academia Sinica Institute of Astronomy & Astrophysics, 11F of Astronomy-Mathematics Building, AS/NTU, No.1, Sec. 4, Roosevelt Rd, Taipei 10617, Taiwan
<sup>23</sup> Department of Astronomy, Van Vleck Observatory, Wesleyan University, 96 Foss Hill Drive, Middletown, CT 06459, USA

<sup>24</sup> Department of Astronomy & Astrophysics, 525 Davey Laboratory, The Pennsylvania State University, University Park, PA 16802, USA
Received 2024 November 26; revised 2025 February 26; accepted 2025 March 14; published 2025 April 28

Abstract

Gas surface density is one of the most relevant physical quantities in protoplanetary disks. However, its precise measurement remains highly challenging due to the lack of a direct tracer. In this study, we report the spatially resolved detection of pressure-broadened line wings in the CO $J = 3-2$ line in the RX J1604.3–2130 A transition disk as part of the exoALMA large program. Since pressure-broadened line wings are sensitive to the total gas volume density, we robustly constrain the radial dependence of the gas surface density and midplane pressure in the region located 50–110 au from the central star, which encompasses the dust ring of the system. The peak radius of the midplane pressure profile matches the dust ring radial location, directly proving radial dust trapping at a gas pressure maximum. The peak gas surface density is 18–44 g cm<sup>-2</sup> and decreases at radii interior to and exterior of the dust ring. A comparison of the gas and dust surface densities suggests that the disk turbulence is as low as $\alpha_{turb} \sim 2 \times 10^{-4}$. Despite dust trapping, the gas-to-dust surface density ratio at the ring peak is 70–400, which implies already-formed protoplanets and/or less efficient dust trapping. The gas surface density drop at radii interior to the ring is consistent with a gas gap induced by a Jupiter-mass planet. The total gas mass within $50 < r < 110$ au is estimated to be $\sim 0.05-0.1 M_{\odot}(50-100 M_{Jup})$, suggesting that planetary system formation is possible.

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

1. Introduction

Planetary systems are born in protoplanetary disks, where gas and dust grains orbit around a young star. To understand the mechanisms of planet formation, it is crucial to observationally reveal the physical structure of these disks. Since the Atacama Large Millimeter/submillimeter Array (ALMA) began operating, it has been revealed that substructures in dust continuum emission in extended systems are ubiquitous

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

<footnote n="25"><sup>25</sup> NASA Hubble Fellowship Program Sagan Fellow.</footnote>
(e.g., S. M. Andrews et al. 2018). Among the different observed substructures, some of the most prominent ones are seen in transition disks, which often exhibit a millimeter dust cavity of about 1–100 au and a radially confined dust ring. It is believed that transition disks appear after the evolution of monolithic full disks via several potential mechanisms, such as grain growth in the inner region, disk clearing by planets, photoevaporation, and/or dust trapping at the outer edge of the magnetorotational instability dead zone (e.g., N. van der Marel 2023). Furthermore, in some cases, the radial accumulation of (sub) millimeter-sized dust grains may cause dust coagulation and could trigger streaming instability, which is one of the most promising formation routes of terrestrial planets and cores of gas giant planets (e.g., A. N. Youdin & J. Goodman 2005; P. Pinilla et al. 2012).

One of the most essential physical quantities to reveal the planet-forming potential in such rings is the gas surface density profile. The gas in such regions is the main ingredient and, as such, it sets the mass budget available for giant planet formation. In addition, the gas affects the dynamics of dust grains. Furthermore, gas dissipation after planet formation may affect the planets' orbits (e.g., P. J. Armitage & W. Kley 2019). Therefore, it is crucial to measure the gas surface density profile. However, despite its importance, measuring it is generally very challenging. This is because the main component of the gas, H₂, does not efficiently emit radiation under typical conditions of protoplanetary disks due to the lack of a permanent dipole moment and high upper-level energies for the rotational ladder. To trace the gas surface density profile, various indirect tracers such as dust thermal emission, optically thin CO isotopologue emission (e.g., A. Miotello et al. 2016), and HD emission (e.g., E. A. Bergin et al. 2013) have conventionally been used. However, even for the total gas mass of one disk, the measured values are significantly scattered (a factor of ~10³) since the conversion factors such as CO/H₂ are highly uncertain (A. Miotello et al. 2023). We note that the relative radial and vertical dependence of the gas pressure (and therefore density) structure can be constrained by analyzing the rotation curve of molecular lines (e.g., R. Teague et al. 2018a, 2018b; G. P. Rosotti et al. 2020; H. Yu et al. 2021; J. Stadler et al. 2025; V. Pezzotta et al. 2025).

Recently, T. C. Yoshida et al. (2022) found very broad wings in the CO *J* = 3–2 line in the inner region of the TW Hya disk and interpreted it as a result of pressure broadening. Since pressure broadening occurs mainly due to collisions between CO and H₂ molecules, the pressure-broadened wings are sensitive to the total gas volume and surface densities. They indeed constrain it without an assumption of the CO/H₂ ratio.

In this Letter, we focus on the transition disk around RX J1604.3–2130 A (hereafter J1604). The central source is a K2 type pre-main-sequence star (K. L. Luhman & E. E. Mamajek 2012) with a mass of $M_* = 1.24 M_\odot$ and an age of ~5–10 Myr (C. F. Manara et al. 2020), which is located at 144.6 pc from the Earth (Gaia Collaboration et al. 2023) in the Upper Scorpius star-forming region (P. De Zeeuw et al. 1999; T. Preibisch & H. Zinnecker 1999). The star is surrounded by a transition disk with a large central cavity (e.g., K. Zhang et al. 2014; S. Mayama et al. 2018; N. van der Marel et al. 2021). The dust continuum emission consists of a single ring at a radius of 81 au from the central star (Figure 1). The intensity of the ring sharply peaks just outside the cavity and steeply decreases toward the outside, suggesting dust trapping at the

<figure><img src="image_4.png" /><figcaption>**Figure 1.** (Top) 331 GHz continuum image. (Bottom) The integrated intensity map of the $^{12}$CO *J* = 3–2 line. Both images are obtained as part of the exoALMA project (P. Curone et al. 2025; R. Teague et al. 2025). The J1604 disk is close to face on and consists of a single well-defined continuum ring.</figcaption></figure>

cavity edge (J. Stadler et al. 2023). This disk is one of the brightest disks in the Upper Scorpius star-forming region in terms of dust and gas emission (G. S. Mathews et al. 2013; J. M. Carpenter et al. 2014; S. A. Barenfeld et al. 2016). R. Dong et al. (2017) observed CO isotopologue lines as well as infrared and millimeter continuum emission and constrained the surface density structure via physical and chemical modeling.

They found that the gas surface density peaks at the dust ring and drops by 3 orders of magnitude in the inner cavity. The peak gas surface density at the ring and the total gas mass were estimated to be ~1 g cm⁻² and 2.5 M<sub>jup</sub>, respectively, although it relies on the potentially optically thick C¹⁸O line. Inside and outside the dust cavity, J. Stadler et al. (2023) detected a velocity deviation of CO gas from the Keplerian rotation, which can be caused by an embedded Jupiter-mass planet. J. B. Bergner et al. (2019) conducted observations of C₂H, HCN, and C¹⁸O lines and detected all molecules, suggesting that the disk is generally molecular rich. As part of the ALMA large program exoALMA (R. Teague et al. 2025), we observed the J1604 disk in the CO *J* = 3–2 line with high angular and spectral resolution and sensitivity, which enables us to search for faint pressure-broadened line wings. We note that this disk is the lowest-inclination disk among the exoALMA sample
(P. Curone et al. 2025), which is preferable for the analysis as we can avoid the beam-smearing effect on the spectral shape. Figure 1 shows the 331 GHz continuum and the integrated CO J = 3–2 line images obtained as part of the exoALMA project. The structure of this Letter is as follows. In the next section, we describe our analysis of the dust continuum observations, which is necessary to interpret the gas observations. In Section 3, we report the wings detected in the ¹²CO line and explain the model-fitting procedure. We present the results and discuss our interpretation in Section 4. We summarize our results in Section 5.

## 2. Analysis of the Dust Emission

Before diving into the analysis of the CO emission line, which is the primary focus of this Letter, we first present our modeling of the dust continuum disk and its corresponding results. The continuum modeling is important for the spectral line fitting (Section 3). Additionally, the derived dust and gas structures are compared in Section 4.

### 2.1. Spectral Energy Distribution Modeling

We analyze the spectral energy distribution (SED) of multiband dust continuum observations to retrieve the dust disk properties, such as the dust surface density profile, maximum dust size, and temperature profile. We use continuum images previously observed with ALMA at 1.3 mm (Band 6) and 3.1 mm (Band 3) in addition to the new exoALMA continuum image at 0.9 mm (Band 7, a beam size of 0.′′14 × 0.′′12; P. Curone et al. 2025). For Band 6, we use the image published by J. Stadler et al. (2023), which has a beam size of 0.′′060 × 0.′′039. We newly analyze archival observations (#2015.1.00819.S; PI L. Ricci) for the Band 3 image, which has a beam size of 2.′′7 × 1.′′7. The reduction processes of the Band 3 data are described in Appendix A. The CLEAN images are reconvolved to have circularized beams of 2.′′7, 0.′′061, and 0.′′14, for Bands 3, 6, and 7, respectively.

The beam sizes are different in the three wavelengths; the Band 6 image clearly resolves the narrow ring structure, but the Band 3 image does not even resolve the whole disk. Therefore, we perform simple parametric modeling to retrieve the dust disk properties. Since the J1604 disk consists of a single ring with a shoulder on the outer side of the ring (S. Mayama et al. 2018; J. Stadler et al. 2023), we assume that the dust surface density $\Sigma_d$ can be expressed by superposition of two Gaussian rings, that is,

$$ \Sigma_d(r) = \Sigma_{d,1} \exp \left[ -\frac{1}{2} \left( \frac{r - r_{d,1}}{\sigma_{d,1}} \right)^2 \right] \quad (1) $$

$$ + \Sigma_{d,2} \exp \left[ -\frac{1}{2} \left( \frac{r - r_{d,2}}{\sigma_{d,2}} \right)^2 \right], \quad (2) $$

where $\Sigma_{d,i}$, $r_{d,i}$, and $\sigma_{d,i}$ represent the peak dust surface density, ring radius, and ring width for the ring index $i = (1, 2)$, respectively. The midplane temperature profile is assumed to be a power law,

$$ T_{\text{mid}}(r) = T_{100}(r/100\text{au})^{-0.5}. \quad (3) $$

In addition, we assume that the dust size distribution is homogeneous in the ring. We set a grain size distribution as a power law with a fixed index of -3.5 and a flexible maximum grain size of $a_{\text{max}}$. Note that the minimum grain size is fixed to

$10^{-5}$ cm, but this assumption does not affect the results. The dust opacities at each wavelength are calculated using dsharp_opac, assuming the DSHARP dust properties (T. Birnstiel et al. 2018).

The dust emission for each location of the disk is then calculated in the same manner as in recent similar studies (e.g., E. Macías et al. 2021; A. Sierra et al. 2021; G. Guidi et al. 2022) based on a solution of the radiative transfer equation when considering self-scattering by dust grains (K. Miyake & Y. Nakagawa 1993; H. B. Liu 2019; Z. Zhu et al. 2019). First, we calculate the intensity of a dust slab at each wavelength as a function of radius according to Equation (11) of Z. Zhu et al. (2019). Then, we assume the geometry of the J1604 disk ($i = 9.°5$, PA = 124.°9, consistent with P. Curone et al. 2025) and distribute the intensity on the 2D image plane. Finally, we convolve the images with the beams and resample the image with intervals of the beam size to properly treat the noise correlation in a beam. Our dust emission model has eight free parameters: six parameters for the surface density profile in Equation (1), $T_{100}$, and $\log_{10}a_{\text{max}}$.

The likelihood function of the model images is calculated by

$$ \chi^2 = \sum_{\nu} \mathbf{R}_{\nu}^{T} C_{\nu}^{-1} \mathbf{R}_{\nu}, \quad (4) $$

where $\mathbf{R}_{\nu}$ is the residual vector between the modeled and observed image at a frequency $\nu$. Note that the observed images are also resampled in the same manner as the model images. $C_{\nu}$ is the variance-covariance matrix at $\nu$. The summation is taken for three frequencies. In general, the absolute flux error dominates the uncertainty in the case of SED analysis. This uncertainty can be treated as a spatially correlated noise over the image assuming that the flux uncertainty follows a known probability distribution. We adopt the following formulation for an $(i, j)$ element of $C$:

$$ C_{i,j} = \sigma^2 \delta_{i,j} + s_f^2 M_i M_j, \quad (5) $$

where $\delta_{i,j}$ is the Kronecker delta, $\sigma$ is the image noise level, $s_f$ is the flux scaling uncertainty (i.e., $s_f = 0.1$ for a flux uncertainty of 10%), and $M_i$ is the $i$th element of the model image.

In practice, we set $s_f = 0$ for all bands first and sampled the posterior probability distribution by the Markov Chain Monte Carlo (MCMC) method using emcee (D. Foreman-Mackey et al. 2013). We adopt 5000 steps with 24 walkers for our MCMC chain. After confirming convergence, we calculate the median values from the last 4200 steps. The model images created by the median values are used for $M_i$ and $M_j$. Then, we allow $s_f$ to be nonzero to estimate the uncertainty of each parameter. While the ALMA technical handbook states that the flux accuracy for flux calibrators is 2.5% for Band 3 and 5% for Bands 6 and 7, the resultant uncertainty on the image plane could be larger due to, for example, phase error and imaging processes. According to L. Francis et al. (2020) and T. C. Yoshida et al. (2025), we conservatively assume that the resultant uncertainty is 3 times larger than the nominal value, substitute $s_f = [0.075, 0.15, 0.15]$ for Bands 3, 6, and 7, respectively, and perform the second round of MCMC sampling. In the second round, the uncertainties of some parameters became larger, as expected. We compare the best-fit model and observations without pixel resampling in Appendix B. The relative difference between observations and models is almost within the flux uncertainties.
<table><caption>**Table 1**
Best-fit Parameters of the Multiband Continuum Modeling</caption><thead><tr><th>Parameter</th><th>Best-fit Value</th><th>Unit</th></tr></thead><tbody><tr><td>r<sub>d,1</sub></td><td>82.81<sup>+0.01</sup><sub>-0.01</sub></td><td>au</td></tr><tr><td>log<sub>10</sub>Σ<sub>d,1</sub></td><td>-0.50<sup>+0.04</sup><sub>-0.03</sub></td><td>g cm<sup>-2</sup></td></tr><tr><td>σ<sub>d,1</sub></td><td>3.09<sup>+0.03</sup><sub>-0.04</sub></td><td>au</td></tr><tr><td>r<sub>d,2</sub></td><td>96.65<sup>+0.09</sup><sub>-0.08</sub></td><td>au</td></tr><tr><td>log<sub>10</sub>Σ<sub>d,2</sub></td><td>-1.31<sup>+0.02</sup><sub>-0.02</sub></td><td>g cm<sup>-2</sup></td></tr><tr><td>σ<sub>d,2</sub></td><td>15.29<sup>+0.05</sup><sub>-0.05</sub></td><td>au</td></tr><tr><td>T<sub>100</sub></td><td>14<sup>+1</sup><sub>-1</sub></td><td>K</td></tr><tr><td>log<sub>10</sub>a<sub>max</sub></td><td>-0.85<sup>+0.05</sup><sub>-0.03</sub></td><td>cm</td></tr></tbody></table>

**Note.** The uncertainties are 16th and 84th percentiles of the posterior probability distributions.

## 2.2. Results of the SED Modeling

The simple model described in Section 2.1 can reproduce the observed images, and the parameters are well constrained. The best-fit parameters of the multiband continuum modeling are presented in Table 1. The resulting dust surface density profile is plotted in Figure 5(a). The peak dust surface density is estimated to be $0.35_{-0.02}^{+0.03}$ g cm$^{-2}$ at a radius of ≈82 au, which is in a typical range found in other (bright) protoplanetary disks (e.g., E. Macías et al. 2021; A. Sierra et al. 2021; G. Guidi et al. 2022). The midplane temperature at $r = 82$ and 100 au is constrained to be $15 \pm 1$ K and $14 \pm 1$ K, respectively. The maximum dust size $a_{max}$ is estimated to be $1.4_{-0.1}^{+0.2}$ mm. Although we assume that $a_{max}$ is uniform across the ring, there is no strong residual in all images (Figure 8), suggesting the assumption is reasonable. From the dust surface density profile, we estimate a total dust mass of $2.7_{-0.1}^{+0.2} \times 10^{-4} M_{\odot}$.

We also test a model with the temperature radial slope being $-0.4$ instead of $-0.5$. While the best-fit temperature increases by ~10%, the peak dust surface density decreases by ~10%. We note that the assumption of the power-law temperature profile might be not very accurate, given that the temperature structure is no longer radially smooth in transition disks (e.g., L. I. Cleeves et al. 2011; S. Bruderer 2013). Higher spatial resolution observations (especially in Band 3) are needed to more robustly determine the dust disk properties.

# 3. Analysis of the CO Emission

In this section, we explain how we analyze the CO line data. The results are presented in Section 4.

## 3.1. Stacked Spectra and Visual Inspection

In the following analysis of the CO line, we use one of the fiducial exoALMA image cubes of the CO $J = 3-2$ line (R. Loomis et al. 2025). The synthesized beam is circularized and has a full width at half-maximum (FWHM) of 0.15. The channel width is set to 0.1 km s$^{-1}$.

To search for the faint wing feature due to pressure broadening, we need to obtain a high signal-to-noise ratio, and therefore we azimuthally stack the spectra. First, we create the line-of-sight velocity map by fitting a tapered Gaussian with an optically thick core to the observed spectrum at each position using better-moments (R. Teague & D. Foreman-Mackey 2018). Then, the velocity map is fitted by the geometrically thin Keplerian rotation

model using eddy (R. Teague 2019) to determine the position angle, stellar mass, and central position of the disk. The inclination angle is fixed to $i = 6.0^{\circ}$ (J. Stadler et al. 2023). We note that a possible error in the inclination angle does not affect the results (Appendix C). Adopting the best-fit parameters, which are consistent with results based on discminer by A. F. Izquierdo et al. (2025), we pick up an annulus with a radial width of 0.05 and fit the velocity variation as a function of the azimuthal angle by

$$v_{\text{LOS}} = v_{\phi} \cos(\phi) \sin(i) - v_r \sin(\phi) \sin(i) + v_{\text{LSR}}, \quad (6)$$

where $v_{LOS}$ and $\phi$ are the line-of-sight velocity and azimuthal angle on the disk coordinate. $v_{LSR}$ denotes the systemic velocity but contains any vertical motion within the disk as well. Fitting parameters $v_{\phi}$ and $v_r$ represent the azimuthal and radial velocity, respectively. We use the module eddy linecube to do this fitting and repeat the procedure for all radii. The original spectra are then shifted by the derived line-of-sight velocities. In Appendix C (Figure 9), we present the shifted spectra as a function of radius, azimuthal angle, and velocity. It is clear that the shifting works well. Finally, we average the spectra along the azimuthal axis and obtain the azimuthally stacked CO spectra, as shown in Figure 2. When stacking the spectra, we make the assumption that the line wings are identical along the azimuthal direction. We have checked that there is no significant azimuthal dependence (Appendix C).

In the stacked spectra, the peak intensity presents its maximum at ~0.4, which is consistent with J. Stadler et al. (2023) and M. Galloway-Sprietsma et al. (2025). Interestingly, the spectra have relatively broad line wings at $r = 0.4 - 0.8$.

The azimuthally stacked spectra are radially averaged at radii between $r = 0.4$ and $r = 0.8$, and are plotted in Figure 3, which shows the wing feature.

In Figure 4, the averaged spectrum is plotted against the absolute velocity on a log-log scale. The wing emission at $|v| \gtrsim 1$ km s$^{-1}$ is almost proportional to $v^{-2}$. If an optical depth profile follows the Gaussian function, the line wings should damp with $\propto \exp(-v^2/2\sigma)$, where $\sigma$ is the line width. Therefore, the observed line profile implies a non-Gaussian line-broadening process.

To quantitatively test if the observed line wing can be explained only by Gaussian broadening, we model the wing emission by an isothermal slab,

$$I = B(T)(1 - e^{-\tau}), \quad (7)$$

where

$$\tau = \tau_0 \exp\left(-\frac{m_{\text{CO}}v^2}{2k_B T}\right), \quad (8)$$

$B$ is the Planck function, $T$ is the temperature, $\tau_0$ denotes the optical depth at the line center, and $m_{\text{CO}}$ and $k_B$ are the molecular mass of CO and the Boltzmann constant, respectively. We fitted the high-velocity range ($|v| > 0.8$ km s$^{-1}$) of the averaged spectrum by Equation (7) using scipy optimize . curve_fit , after converting the units to brightness temperature. The best-fit parameters are $T \simeq 4.8 \times 10^3$ K and $\tau_0 \simeq 1.9 \times 10^{-4}$. The best-fit spectrum is also plotted with the orange dotted line in Figure 4, which matches the observations at small $v$ but deviates at larger $v$. The best-fit
<figure><img src="image_5.png" /><figcaption>**Figure 2.** (Left) Azimuthally stacked CO spectra. The upper and lower panels show the same data, but in the lower panel, we convolved the spectrum at each radius with a flat kernel of three channels and use the symmetric logarithm color map in order to enhance the faint emission in the line wings. The horizontal dotted lines indicate the radial range of the dust ring. The dashed ellipse marks the faint line wings we analyze in this Letter. (Middle) The best-fit model of the CO spectra. (Right) The residual between the observed and model spectra.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 3.** Radially averaged spectrum at $r = 0.4'' - 0.8''$ (see Figure 2). The lower panel shows a zoomed-in view of the upper panel. Colored lines indicate models with different gas surface densities, as discussed in Section 4.1.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 4.** Log-log plot of the radially averaged spectrum at $r = 0.4'' - 0.8''$ as a function of $|v|$. The wing ($|v| > 1 \text{ km s}^{-1}$) emission follows $v^{-2}$, which is predicted for pressure-broadened line wings. The orange dotted line indicates the best-fit model if we interpret the wing as a high-temperature component. As discussed in Section 3.1, the derived temperature from the line width is $4.8 \times 10^3 \text{ K}$, which is unrealistically high for CO molecules.</figcaption></figure>

temperature of $T \approx 4.8 \times 10^3 \text{ K}$ might suggest that the emission originates from the uppermost layer of the disk. However, since CO molecules are thermally dissociated at ~2000 K, such a high temperature is unphysical. Therefore, the observed broad line wings cannot be explained by a high-temperature component. We note that subsonic turbulent
broadening (e.g., K. M. Flaherty et al. 2018) is also unlikely as a cause of the observed wings because such broadening effects should be substantially smaller than thermal broadening, which means that the very high temperature is still needed even if turbulence exists. In addition to turbulence, we note that any microscopic non-Keplerian motion can be excluded as long as it is subsonic.

Larger-scale non-Keplerian motion of the CO gas could also enhance the line width. However, the velocity deviation from the Keplerian rotation does not exceed ~0.1 km s⁻¹ at these radii (J. Stadler et al. 2023; R. Teague et al. 2025), which is much smaller than the observed width.

Another potential scenario is a disk wind that is vertically mirror symmetric against the disk midplane. However, the continuum modeling (Section 2.1) suggests that the optical depth of the ~80 au continuum ring is ~5 in Band 7. Therefore, the emission from the redshifted side of the disk wind should be absorbed by the continuum disk. The observed symmetric spectrum shows that the disk winds cannot explain the wings.

A candidate of the function that has a damping wing ∝ v⁻² is the Lorentzian function, or its convolution with a Gaussian, the Voigt function. The most promising mechanism that makes the line profile Voigt-like is pressure broadening (e.g., G. B. Rybicki & A. P. Lightman 1991; T. C. Yoshida et al. 2022). Pressure broadening is a general process of line broadening in astrophysics. When a CO molecule collides with another molecule, mainly H₂ or He, the radiation field is disrupted. If the collision randomly occurs by following the Poisson probability distribution as a function of time, the line profile becomes the Lorentzian function. Since the line profile without collision is the Gaussian function, the resultant line profile is the convolution of those two functions, the Voigt function. Indeed, T. C. Yoshida et al. (2022) found broad line wings in the TW Hya disk and reproduced the observed spectrum by taking pressure broadening into account. Importantly, the intensity of the pressure-broadened line wings depends on the gas volume density near the midplane, which is challenging to trace with other methods.

### 3.2. Modeling the Spectra

To derive the gas surface density profile, we model the azimuthally stacked spectra after considering pressure broadening. First, we assume the vertical temperature structure. We used a parametric function that describes a layered temperature structure (e.g., E. Dartois et al. 2003; S. M. Andrews et al. 2012; C. J. Law et al. 2022):

$$T(z) = T_0 \qquad (9)$$

$$+ \left\{ \begin{array}{ll} T_{\text{mid}} & (|z| < z_1) \\ \frac{T_{\text{atm}} - T_{\text{mid}}}{2} \cos\left(\frac{\pi |z| - z_1}{2 z_2 - z_1}\right) & (z_1 \le |z| < z_2), \\ T_{\text{atm}} & (z_2 \le |z|) \end{array} \right. \qquad (10)$$

Here, z is the height from the disk midplane. T_atm represents the temperature in the atmosphere, and z₁ and z₂ are the characteristic heights of transition from the midplane to the transition layer and the transition layer to the atmosphere, respectively. T mid is given in Equation (3). We note that z₁ and z₂ are practically given, as follows:

$$z_1 = f_{z_1} H_g \qquad (11)$$

$$z_2 = f_{z_2} H_g \qquad (12)$$

Here, f_z1 and f_zw ≡ f_z₂ − f_z₁ are free parameters in the fitting. H_g is the gas scale height and is given by

$$H_g = \sqrt{\frac{k_B T_{\text{mid}} r^3}{\mu m_p G M_*}}, \qquad (13)$$

where G is the gravitational constant, and m_p denotes the proton mass. μ is the mean molecular mass and is set to 2.37. Additionally, we introduced

$$T_0 = \begin{cases} 0 & (|z| < z_1) \\ p \frac{|z| - z_1}{z_2} & (|z| \ge z_1) \end{cases} \qquad (14)$$

The background temperature gradient in addition to the two-layer structure with arctangent-like transition has been seen in radiative transfer calculations (C. P. Dullemond et al. 2002). Indeed, this parameter is needed to fit the core component of the spectra based on our tests.

In our model-fitting procedure, we fix T₁₀₀ = 14 K as a fiducial case. This is not only based on our continuum modeling in Section 2.2, but also a typical value around r = 100 au in a T Tauri disk without a cavity (e.g., C. P. Dullemond et al. 2002). However, considering that the dust disk modeling presents its own uncertainties (Section 2.2), this value may be underestimated for three reasons. First, the detection of the pressure-broadened line wings suggests that the gaseous CO molecules survive around the disk midplane below one gas scale height, where the gas density is high. However, since the CO sublimation temperature is generally thought to be ~20 K (e.g., M. Minissale et al. 2022), the CO molecules may be frozen near the midplane, although we note that observations imply some diversity depending on the disks (~13–28 K; C. Qi et al. 2013; K. R. Schwarz et al. 2016; K. Zhang et al. 2017; C. Pinte et al. 2018; C. Qi et al. 2019). Second, it is theoretically predicted that the midplane temperature rises to ~30 K around a cavity edge of a transition disk (L. I. Cleeves et al. 2011; S. Bruderer 2013), which should be in similar conditions to the J1604 disk. Third, a disk survey by J. B. Bergner et al. (2019) suggests that the J1604 disk has the highest excitation temperatures of the HCN and C₂H lines among their samples that include Herbig disks (see also Appendix H). Although these lines may trace relatively elevated layers (e.g., T. Paneque-Carreño et al. 2023), it is likely that the midplane temperature is also relatively higher than typical T Tauri disks. Therefore, we also consider another case with T₁₀₀ = 30 K as the higher end of the midplane temperature.

The H₂ gas density structure as a function of height, nH₂(z), is determined by assuming the vertical hydrostatic equilibrium self-consistent with the temperature structure (C. Longarini et al. 2025). The CO density structure is written as

$$n_{\text{CO}}(z) = X_{\text{CO}} n_{\text{H}_2}(z), \qquad (15)$$

where XCO is the CO/H₂ abundance ratio. We tentatively assume XCO = 10⁻⁴ for fitting and discuss possible uncertainties in the following sections.

We integrate the radiative transfer equation using the above CO density and temperature structures, and then create a synthetic spectrum. Here, we also take the continuum absorption into account. In practice, we first calculate the line emission from the back side of the disk. Then, we multiply it by
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Fitting results for the fiducial temperature assumption. (a) Gas (blue) and dust (orange) surface density profiles. The dust surface density is multiplied by 100 for visibility purposes. The black bar at the upper-right corner indicates the spatial resolution of the CO image. (b) Midplane gas pressure profile (blue) with an arbitrarily scaled dust surface density profile (orange). The dust surface density peaks at the gas pressure maximum. Note that there are uncertainties from the undetermined $X_{\text{CO}}$ in panels (a) and (b). (c) Gas surface density profile for various $X_{\text{CO}}$ with criteria for gravitational instability and streaming instability.</figcaption></figure>

a factor of $\exp(-\tau_d)$, where $\tau_d$ is the corresponding optical depth of the continuum emission at 346 GHz derived by the multiband continuum modeling. Finally, we add the continuum intensity and restart the integration of the radiative transfer equation for the front side to take the potential absorption by the continuum disk into account. The continuum emission is then subtracted from the spectra for a fair comparison.

For the line opacity profile, we adopt the Voigt function and follow T. C. Yoshida et al. (2022). For the spectroscopic constants and broadening parameters, we use the LAMDA database (F. L. Schöier et al. 2005; B. Yang et al. 2010) and HITRAN database (I. E. Gordon et al. 2022; Y. Tan et al. 2022), respectively. We approximate that the (bulk) velocity of CO gas is constant along the line of sight at each position since the disk is nearly face on (B. Lankhaar & R. Teague 2023). In addition, we convolve the spectrum with a Gaussian that represents the effective velocity resolution. Note that the velocity resolution is limited not only by the instrumental resolution but also by the intensity gradient in a beam (M. Keppler et al. 2019; V. Pezzotta et al. 2025). Finally, we shift the spectrum along the velocity axis by $\Delta v$ and interpolate at the observed velocity channel.

The convolved synthetic spectrum at each radius $r$ is then compared with the azimuthally stacked spectrum at the same radius. Here, there are seven free parameters, which are the gas surface density $\log_{10}\Sigma_g$, the difference between atmospheric temperature and the midplane temperature $\Delta T = T_{\text{atm}} - T_{\text{mid}}$, parameters on the temperature structure, $f_{z_t}$, $f_{z_w}$, and $p$, velocity shift $\Delta v$, and the effective velocity resolution $\sigma_v$. We sample the posterior probability distribution of the model spectrum at each radius using the MCMC method with emcee (D. Foreman-Mackey et al. 2013). We set uniform priors for all parameters with boundaries of $\log_{10}\Sigma_g = [-3, 3] \text{ g cm}^{-2}$, $\Delta T = [0, 2000] \text{ K}$, $f_{z_t} = [0, 50]$, $f_{z_w} = [0, 50]$, $p = [0, 50]$, $\Delta v = [-0.1, 0.1] \text{ km s}^{-1}$, and $\sigma_v = [0.02, 1.0] \text{ km s}^{-1}$. For each radius, our MCMC chain consists of 2000 steps with 32 walkers. We run the sampling process, discard the first 1000 steps, and calculate the median values. Then, we initialize the chains again using the median values. We rerun MCMC sampling for 6000 steps and check that the chains are well converged. Finally, after discarding the first 3000 steps, the median values of the chain are adopted as fiducial parameters. We present the typical marginalized posterior distribution in Appendix D (Figure 12).

Furthermore, to check if pressure broadening is indeed required, we turn off the Voigt function but use the Gaussian

function as the opacity profile and perform the same fitting procedure as in the above case. The resultant best-fit model and its residual are plotted in Appendix E. We confirm that the spectra cannot be well fitted without pressure broadening.

# 4. Results and Discussion

## 4.1. Gas Surface Density Profile

We compare the best-fit CO spectrum at $r = 0.4 - 0.8$ for the fiducial temperature ($T_{100} = 14$ K) with observations in Figure 3. Note that the best-fit spectrum for the high-temperature case with $T_{100} = 30$ K is identical to the above case. If we reduce the gas surface density by a factor of 2, the observed wings are not reproduced. On the other hand, if we scale the gas surface density by a factor of 2, the wings are much brighter than the observations (colored lines in Figure 3). Clearly, the broad line wings are extremely sensitive to the gas surface density. This is because the optical depth in the wing is proportional to the squared gas volume density, and the wings are optically thin (T. C. Yoshida et al. 2022). We discuss the CO-emitting layer as a function of $v$ in Appendix F. The two-dimensional temperature structure simultaneously constrained by the model fitting is discussed in Appendix G.

We plot the best-fit gas and dust surface density profiles in Figures 5(a) and 6(a) for the fiducial and higher-temperature cases, respectively. In both cases, the gas surface density decreases at the radii within and outside the dust ring. The maximum value at $r \sim 87$ au is $\sim 80(X_{\text{CO}}/10^{-4})^{-0.5} \text{ g cm}^{-2}$ and $\sim 56(X_{\text{CO}}/10^{-4})^{-0.5} \text{ g cm}^{-2}$ for the fiducial and high-temperature cases, respectively. A factor of $(X_{\text{CO}}/10^{-4})^{-0.5}$ is included here because we assume $X_{\text{CO}} = 10^{-4}$ in the model fitting, and the intensity in the pressure-broadened line wings depends on the product of the squared gas surface density and the CO/H$_2$ ratio (T. C. Yoshida et al. 2022). This factor can be reduced by measuring the column density of CO via observations of optically thin CO isotopologue lines such as $^{13}\text{C}^{18}\text{O}$ or $^{13}\text{C}^{17}\text{O}$. As of writing this Letter, there is no observation of these isotopologues in the J1604 disk. The $^{13}\text{CO}$ line observed by the exoALMA project is significantly optically thick (M. Galloway-Sprietsma et al. 2025). We also check the $^{18}\text{O}$ line $J=2-1$ in the ALMA Science Archive, and the peak brightness temperature around the ring exceeds $\sim 30$ K in the full Planck function. This suggests that the $^{18}\text{O}$ line is also highly optically thick (Appendix H). It is also notable that we can eliminate the uncertainty from the temperature assumption by observing an optically thin isotopologue line because such lines should have the same source function as the line
<figure><img src="image_9.png" /><figcaption>Figure 6. The same as Figure 5, but for the high midplane temperature case with T$_{100}$ = 30 K.</figcaption></figure>

wing of $^{12}$CO and the ratio of them determines the gas density (T. C. Yoshida et al. 2022).

Future observations will allow us to robustly determine the absolute value of the gas surface density profile.

However, there are some constraints on XCO (K. I. Öberg et al. 2023). According to detailed chemical network calculations with dust evolution by S. Krijt et al. (2018, 2020), the gas phase XCO can vary from 10⁻³ to 10⁻⁶ in protoplanetary disks. The lower value can result from the locking of gas phase CO to large dust grains and their radial drift. On the other hand, the higher value may arise interior to the CO snowline due to transportation and sublimation of the CO ice from the outer disk. Observations by K. Zhang et al. (2020b) also suggest that XCO decreases from 10⁻⁴ to 10⁻⁶ as a function of age. If we take a conservative range of 10⁻⁶ < XCO < 10⁻³, the remaining factor in the gas surface density ranges from 0.3 to 10. Therefore, we can reduce the uncertainty by a factor of ~30 compared to conventional gas surface density measurements using CO isotopologue lines that are proportional to XCO.

We calculate Toomre's Q (A. Toomre 1964) from the inferred gas surface density profiles. In the case of the fiducial temperature assumption, Toomre's Q is estimated to be lower than unity at r > 70 au if XCO = 10⁻⁴, which might indicate that the disk is gravitationally unstable. However, there is no significant asymmetric structure such as a spiral arm in the J1604 disk, suggesting the disk is gravitationally stable and Q > 1 in reality. If we assume Q > 1 at the peak of the gas surface density, we can obtain XCO > 7 × 10⁻⁴. By combining the conservative range of XCO, it is suggested that the CO is not depleted around the ring, with 7 × 10⁻⁴ < XCO < 10⁻³. Even the lower limit is higher than the elemental abundance of carbon in the interstellar medium (ISM), C/H ≈ 3.2 × 10⁻⁴ (B. S. Hensley & B. T. Draine 2021), implying an additional supply of carbon by the radial drift and trapping of dust grains (S. Krijt et al. 2018; K. Zhang et al. 2020a). Adopting the range of XCO, the peak gas surface density is constrained to be 25–30 g cm⁻². On the other hand, assuming the high midplane temperature, the peak gas surface density is constrained to be 18–44 g cm⁻² by similarly assuming that the disk is gravitationally stable. In addition, XCO can be estimated to be 1.6 × 10⁻⁴ < XCO < 10⁻³ in this case. The lower side becomes consistent with the interstellar abundance. Considering both temperature assumptions, we conclude the peak gas surface density of 18–44 g cm⁻² with 1.6 × 10⁻⁴ < XCO < 10⁻³.

We can also estimate the gas mass within 50 < r < 110 au of 0.05–0.1 M⊙ (or 50–100 M$_{Jup}$) by integrating the derived gas surface density profile and assuming XCO is constant over this radial range. The gas mass within the continuum ring

(75 < r < 100 au) is also estimated to be ~0.02–0.06 M⊙ (or ~20–60 M$_{Jup}$). As the mass is much larger than several Jupiter masses, we conclude that there is enough gas for planetary system formation in the J1604 disk. The gas surface density and mass is at least ~20 times larger than the estimate by R. Dong et al. (2017). This discrepancy may be explained by the optically thick C¹⁸O line (Appendix H), which was effectively used as a mass tracer in R. Dong et al. (2017).

We also show the midplane gas pressure profile,

$$P_{\text{mid}} = n_{\text{H}_2}(0)k_{\text{B}}T, \qquad (16)$$

in Figures 5(b) and 6(b) for the different temperature cases. The P$_{mid}$ profile peaks at r ~ 80 au with ~21(X$_{CO}$/10⁻⁴)⁻⁰.⁵ µPa and ~23(X$_{CO}$/10⁻⁴)⁻⁰.⁵ µPa for the fiducial and higher-temperature cases, respectively, which matches the peak of the dust surface density. This is direct evidence that the dust grains are radially trapped at the pressure maximum. In addition, there is no detection of a significant substructure in the gas and dust surface density beyond the peak. Therefore, it is suggested that dust grains radially drifted from the outer region toward the ring. However, we estimate the gas-to-dust surface density ratio at the dust surface density peak to be 70–400. Here, the dust surface density for the higher-temperature case is estimated by multiplying B(T = 14 K)/B(T = 30 K), where B(T) is the Planck function, to the fiducial case to compensate for the temperature difference. This gas-to-dust surface density ratio is broadly similar to or larger than that in the interstellar medium. In the case of dust trapping in gas pressure maxima, however, significant dust enrichment by a factor of 10–100 is predicted (e.g., T. Birnstiel et al. 2010). The moderate or no dust enrichment in dust rings could be explained by two scenarios. First, it is possible that small dust grains can still be advected across the ring, suppressing the dust enrichment. Second, there might be significant hidden masses, such as planetesimals, that cannot be observed at millimeter wavelengths (S. M. Stammler et al. 2019).

The gas surface density drops by a factor of ~2–2.5 at r ~ 60 au compared to the ring peak. This is consistent with J. Stadler et al. (2025), who found evidence of a gas density cavity in the inner region by analyzing the rotation curves of $^{12}$CO and $^{13}$CO. This gas drop could be interpreted as a gap induced by a giant planet. Indeed, J. Stadler et al. (2023) suggest that the velocity deviation in this disk can be explained by a giant planet with a mass of 1.6–2.9 M$_{Jup}$ at r = 41 ± 10 au. Assuming that the background gas surface density at r ~ 60 au (before being carved by the planet) is ~110(X$_{CO}$/10⁻⁴)⁻⁰.⁵ g cm⁻² by a rough
<figure><img src="image_10.png" /><figcaption>Figure 7. Observed dust surface density profile (black) and calculated profiles from the observed gas surface density for different α<sub>turb</sub>.</figcaption></figure>

extrapolation from the outer region, the depth of the gas gap is estimated to be ~0.3. According to K. D. Kanagawa et al. (2015), we can estimate the mass of the planet of 0.1–2 M<sub>jup</sub> for a range of turbulence parameters α = 10<sup>-4</sup>–10<sup>-2</sup>, which is also roughly consistent with the results of J. Stadler et al. (2023). Note that J. Stadler et al. (2023) fixed the disk aspect ratio at the planet location to 0.1 while we used self-consistent values. If we adopt the disk aspect ratio of 0.1 instead, the mass of the planet is estimated to be 0.3–3 M<sub>jup</sub>.

## 4.2. Turbulence and Dust Fragmentation

The degree of accumulation of the dust grains at the gas pressure maximum depends on the turbulence strength. Therefore, with the derived gas and dust surface density profiles, we can estimate the disk turbulence strength (C. P. Dullemond et al. 2018; G. P. Rosotti et al. 2020). The gas surface density profile can be related to the dust surface density profile by the advection-diffusion equation. We can write its steady-state solution as

$$ \Sigma_d = \epsilon_0 \Sigma_g \exp \left[ \int_{r_0}^{r} \frac{St}{\alpha_{\text{turb}}} \frac{d \ln P_{\text{mid}}}{dr} dr \right], \quad (17) $$

where $\epsilon_0$ is the dust-to-gas surface density ratio at an arbitrary radius $r_0$, and $\alpha_{\text{turb}}$ is the turbulence strength (A. Sierra et al. 2019). St is the Stokes parameter and is given by

$$ St = \frac{\pi \rho_m a_{\text{rep}}}{2 \Sigma_g}. \quad (18) $$

Here, $\rho_m$ is the dust material density. We introduce a representative dust radius $a_{\text{rep}}$, to which the observations are most sensitive, rather than the maximum dust size. We adopt $a_{\text{rep}} = 0.6$ mm by converting $a_{\text{max}} = 1.4$ mm, which is derived from the multiband continuum modeling, using Figure 10 of K. Doi & A. Kataoka (2023).

Using the observed gas surface density and midplane pressure profiles for the fiducial temperature case with $X_{\text{CO}} = 7 \times 10^{-4}$, we calculate $\Sigma_d$ for $\alpha_{\text{turb}} = [1 \times 10^{-3}, 2 \times 10^{-4}, 4 \times 10^{-5}]$ and compare them with the observed dust surface density profile in Figure 7. Note that the peak dust surface densities are always normalized to 0.35 g cm<sup>-2</sup> as in the observations.

When $\alpha_{\text{turb}}$ is set to $2 \times 10^{-4}$, the observed width of the dust ring is well reproduced. Such a weak level of turbulence near the midplane is also inferred by multiple observations (e.g., K. M. Flaherty et al. 2018; K. Flaherty et al. 2020; M. Villenave et al. 2022; T. C. Yoshida et al. 2022) and might be common among disks. Note that the assumption of the midplane temperature does not affect the widths significantly, although the predicted dust peak slightly shifts toward the disk center for the higher-temperature case.

With the turbulence strength and Stokes parameter, we can estimate the dust fragmentation threshold velocity $v_{\text{frag}}$. Assuming that the maximum grain size is determined by the turbulence-induced collisional fragmentation, $v_{\text{frag}}$ is expressed as

$$ v_{\text{frag}} = \sqrt{3St_{\text{max}}\alpha_{\text{turb}}} c_s, \quad (19) $$

where $St_{\text{max}}$ is the Stokes parameter for the maximum dust size $a_{\text{max}}$, and $c_s$ is the sound speed (e.g., T. Birnstiel et al. 2009; T. Ueda et al. 2019, 2020; F. Zagaria et al. 2023). Using the constrained values, $v_{\text{frag}}$ is estimated to be ~1 m s<sup>-1</sup> at the dust ring. This low fragmentation velocity is roughly consistent with recent observational results (H. Jiang et al. 2024; T. Ueda et al. 2024; T. C. Yoshida et al. 2025) and may be explained by grains covered by less sticky material (G. Musiolik et al. 2016; S. Okuzumi & R. Tazaki 2019; S. Arakawa & S. Krijt 2021). Another possibility is that the maximum dust size can be limited by the bouncing, not the fragmentation, which also may explain the small dust size (e.g., C. Dominik & C. P. Dullemond 2024). We note that this estimate also depends on the accuracy of $a_{\text{max}}$, which could be improved by high-resolution observations in longer wavelengths.

## 4.3. Comparison with Other Targets

As part of the exoALMA large program, several independent attempts were made to measure disk masses: C. Longarini et al. (2025) used rotation curves to determine dynamical masses; L. Trapman et al. (2025) forward-modeled the $N_2H^+$ and C<sup>18</sup>O; and, finally, G. Rosotti et al. (2025) leveraged the CO emission height as a proxy for disk mass. Unfortunately, the J1604 disk is excluded from their samples because of its low inclination, which makes kinematical measurements difficult. However, it may be interesting to compare the J1604 disk with other exoALMA targets (R. Teague et al. 2025).

The estimated mass within $50 < r < 110$ au of the J1604 disk is a lower limit of the total mass since the CO emission is detected up to $r \sim 265$ au (J. Stadler et al. 2023) while the contribution from the inner cavity may be negligible as the <sup>12</sup>CO emission is not strong (R. Teague et al. 2025). Assuming that the gas surface density is ~10 g cm<sup>-2</sup> and zero at $r \sim 110$ au and 265 au, respectively, and the profile linearly decreases as a function of radius, we can estimate the gas mass beyond $r \sim 110$ au to be ~0.09 M<sub>☉</sub> and the total gas mass to be ~0.1–0.2 M<sub>☉</sub>. We note that this value may be an upper limit as the gas surface density might decrease more steeply. Therefore, we can roughly estimate the total gas mass to be 0.05–0.2 M<sub>☉</sub> with the lower and upper limits.

This overlaps the mass estimates of other samples from the kinematical analysis (C. Longarini et al. 2025) and is toward the higher end of values from the molecular line analysis (L. Trapman et al. 2025). If the actual mass is closer to the upper limit, the J1604 disk is one of the most massive disks in the sample.
In this study, the CO/H₂ ratio in the ring of the J1604 disk is found to be similar to or higher than the ISM value. Among the samples, T Tauri disks are found to have typically depleted in CO 3–30 times compared to the ISM (G. Rosotti et al. 2025; L. Trapman et al. 2025), which is different from the J1604 disk. This might imply that there is radial variation in the CO/H₂ ratio; in the inner region, CO/H₂ is high due to thermal desorption and additional supply from radially drifted ice, while it may be depleted in the outer region (e.g., K. Zhang et al. 2019, 2021).

### 4.4. Caveats

In our analysis, we assumed that CO molecules are present in the gas phase even in the disk midplane and that the CO/H₂ ratio is vertically constant. If the CO molecules were frozen out near the midplane, the derived gas surface density would be underestimated.

Regarding the gas-to-dust ratio, it should be noted that there might be a large uncertainty in the dust opacity model. For instance, the porous dust model with the maximum dust size of 1 mm in T. Birnstiel et al. (2018) has a ~6 times smaller absorption opacity at a 3 mm wavelength. Therefore, if we adopted the porous dust model, the gas-to-dust surface density ratio at the ring peak would become 10–70. This value satisfies the favorable condition for streaming instability (A. N. Youdin & J. Goodman 2005). On the other hand, it has been recently proposed that higher absorption opacity models such as L. Ricci et al. (2010) better match observations (L. Delussu et al. 2024). In this case, the gas-to-dust ratio would become even larger than the results based on the DSHARP model.

Furthermore, we note that the spatial resolution of our observations is limited. Therefore, the width of the gas surface density bump could be narrower in reality, and caution needs to be exercised when comparing with the dust ring. We note, however, that the gas ring is spatially resolved at our resolution, which means that this effect is not very significant.

In the CO spectra fitting, we assumed that the stacked spectrum is identical to the intrinsic one without stacking. However, these two could differ in reality due to the inclination effect, although we approximated this effect by adopting the effective velocity resolution. This might result in a potential error in the atmospheric temperature (Appendix G), which is effectively constrained from the narrow line core. However, we emphasize that the atmospheric temperature is almost independent of the gas surface density (Appendix D).

## 5. Summary

Thanks to the high-sensitivity observations of the exoALMA large program, we detected pressure-broadened ¹²CO line wings in the J1604 disk. This is the second detection of pressure broadening in protoplanetary disks, opening up further application of the method demonstrated by T. C. Yoshida et al. (2022) and enabling us to robustly derive the gas surface density profile. Notably, this method can achieve much higher accuracy by at least a factor of ~30 compared to a conventional method that utilizes a CO isotopologue line intensity. We modeled the ¹²CO emission line by solving the radiative transfer equation and fitted it to the azimuthally averaged spectra. Based on the analysis, we conclude the following:

1. We derived the gas surface density and midplane pressure profiles by using the pressure-broadening effect. The pressure profile has a clear peak at the dust ring radius (~82 au). This is direct evidence of radial dust trapping at the gas pressure maximum.

2. Assuming that the disk is gravitationally stable and the CO/H₂ ratio is within 10⁻³–10⁻⁶, the peak gas surface density is constrained to be 18–44 g cm⁻². The disk gas mass within 50 < r < 110 au and within the ring are estimated to be 0.02–0.06 M☉ and ~0.05–0.1 M☉, respectively, suggesting that there is a sufficient mass for planet formation to occur. The estimated mass range is comparable with the other disks in the exoALMA sample.

3. The gas surface density profile exhibits a drop at r ~ 60 au. If we attribute this to a planetary-induced gap, we can estimate the planetary mass to be 0.1–2 M_jup.

4. We analyzed the multiband continuum observations with a simple parametric model and derived the dust surface density profile. Comparing the gas and dust surface densities, we find that the ring is not very enriched in dust. This suggests that there is significant hidden dust mass and/or less efficient dust trapping.

5. From the dust ring width compared to the gas ring width, we constrained the turbulence strength in the ring to be α_turb ~ 2 × 10⁻⁴, which is consistent with previous studies in other disks.

6. Assuming that the maximum dust size is regulated by the turbulence-induced collisional fragmentation, the fragmentation threshold velocity $v_{frag}$ is estimated to be ~1 m s⁻¹, which suggests that the dust grains are fragile.

We note that the CO/H₂ ratio remains uncertain in this study. However, we can eliminate this uncertainty and even determine the CO/H₂ ratio if the CO optical depth is known. Observations of optically thin CO isotopologue lines are required. The methodology used in this Letter should be applicable to other disks with similar characteristics (such as inclination, disk mass, and temperature), providing a new window for high-precision studies of the physical properties of planet-forming environments.

## Acknowledgments

This Letter makes use of the following ALMA data: ADS/JAO. ALMA #2013.1.01020.S, #2015.1.00819.S, #2015.1.00964.S, #2017.A.01255.S, and #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. T.C.Y. is supported by a Grant-in-Aid for JSPS Fellows, JP23KJ1008. 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. 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 Côte d'Azur. M.M. is supported by a Grant-in-Aid from the Japan Society for the Promotion of Science (KAKENHI: No. JP18H05441). J.B. acknowledges support from NASA XRP grant No. 80NSSC23K1312. 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). P.C. acknowledges support by the Italian Ministero dell'Istruzione, Università e Ricerca through the grant Progetti Premiali 2012—iALMA (CUP C52I13000140001) and by the ANID BASAL project FB210003. 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). M. Fukagawa is supported by a Grant-in-Aid from the Japan Society for the Promotion of Science (KAKENHI: No. JP22H01274). 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. Support for A.F.I. was provided by NASA through the NASA Hubble Fellowship grant No. HST-HF2-51532.001-A awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS5-26555. C.L. has received funding from the European Union's Horizon 2020 research and innovation program under the Marie Skłodowska-Curie grant agreement No. 823823 (DUSTBUSTERS) and by the UK Science and Technology research Council (STFC) via the consolidated grant ST/W000997/1. C.P. acknowledges Australian Research Council funding via FT170100040, DP18010423, DP220103767, and DP240103290. D.J.P. acknowledges Australian Research Council funding via DP18010423, DP220103767, and DP240103290. G.R. acknowledges funding from the Fondazione Cariplo, grant No. 2022-1217, and the European Research Council (ERC) under the European Union's Horizon Europe Research & Innovation Programme under grant agreement No. 101039651 (DiscEvol). 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. H.-W.Y. acknowledges support from National Science and Technology Council (NSTC) in Taiwan through grant NSTC 113-2112-M-001-035 and from the Academia Sinica

Career Development Award (AS-CDA-111-M03). G.W.-F. acknowledges support from the European Research Council (ERC) under the European Union Horizon 2020 research and innovation program (grant agreement No. 815559 (MHDiscs)). G. W.-F. was granted access to the HPC resources of IDRIS under the allocation A0120402231 made by GENCI. 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. Support for B.Z. was provided by The Brinson Foundation.

Facility: ALMA.

Software: astropy (Astropy Collaboration et al. 2013, 2018), emcee (D. Foreman-Mackey et al. 2013), eddy (R. Teague 2019).

## Appendix A
### Data Reduction of the Band 3 Image

We obtained the visibility data from the ALMA Science Archive (#2015.1.00819.S; PI L. Ricci). The J1604 disk was observed on 2016 January 20 in Band 3 with 41 antennas. The baseline lengths range from 15 to 331 m, and the on-source integration time was 2.8 minutes. All the spectral windows are used in the Time Division Mode.

We first ran the calibration script provided by the observatory using the Common Astronomical Software Application (CASA; J. P. McMullin et al. 2007) pipeline v4.5.2. CASA modular v6.5.5 was used throughout the remaining analysis. We split out only the data observing the J1604 disk from the calibrated data set. We then CLEANed the data and performed iterative self-calibration. Since the relative bandwidth with respect to the observing frequency is large, we adopted the MTMFS mode for CLEAN. The Briggs weighting with robust = 0.5 was applied. Three rounds of phase self-calibration with solution intervals of the durations of the observation, 90 s, and 30 s were applied. The final rms noise level is measured to be 0.05 mJy, and the resultant beam size is $2.7' \times 1.7'$.

## Appendix B
### Multiband Continuum Modeling Results

We compare the observed and modeled continuum images in Figure 8. Since the absolute flux uncertainty dominates the spatial noise, we also plot the relative error between the observed and modeled images.
<figure><img src="image_11.png" /><figcaption>**Figure 8.** Observed continuum images, best-fit model images, and relative error between the observations and model images. The white ellipses on the left column indicate the observational beam size. The error maps are masked by the $3\sigma$ noise level measured in the observation images.</figcaption></figure>

## Appendix C
CO Line Stacking

In Figure 9, we present shifted spectra for each annulus as a function of velocity from the line center and azimuthal angle.

Figure 10 shows radially averaged azimuthally stacked spectra at $r = 0.4 - 0.8$ AU. Here, the azimuthal stacking is performed only over one quarter of an annulus. In the right panel, we show the residual of observed spectra after subtracting the fiducial model. There is no significant signal in the wing region (outside of the blue-shaded range) of the residual. In the line core, some residuals remain. This is due to the azimuthal asymmetry of

the peak intensity, which might be caused by the shadow cast by the center region (J. Stadler et al. 2023).

When we stacked the spectra, we fixed the disk inclination angle to 6.0° according to J. Stadler et al. (2023) and fitted other geometrical parameters (Section 3.1). To test if the wing feature persists for different possible inclination angles, we changed the inclination to 1.0°, 4.0°, 8.0°, and 11.0° and generated the teardrop plots again (Figure 11). Here, we fixed other parameters such as the position angle to the fiducial values used in the main text. As a result, no significant difference is observed, which means that our analysis is robust for the uncertainty in the geometrical parameters.
<figure><img src="image_12.png" /><figcaption>**Figure 9.** Shifted spectra for each annulus. In each panel, the spectra are plotted as a function of velocity from the line center and azimuthal angle in the disk coordinate. The numbers in the upper-left corner indicate the radius of the annulus in arcseconds.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 10.** (Left) Radially averaged azimuthally stacked spectra at $r = 0''4 - 0''8$. The stacked azimuthal range is limited to one-quarter of an annulus, and the center azimuthal angle is shown in the legend. The gray solid lines indicate the fiducial model spectrum. (Right) Residual of observed spectra after subtracting the fiducial model. Note that the spectra are vertically shifted for visibility purposes. The blue-shaded area indicates the optically thick line cores.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>**Figure 11.** Teardrop plots for different inclinations. The wing feature persists despite 5° differences in the input inclination. Unlike Figure 2, the velocity convolution is not applied for completeness.</figcaption></figure>
## Appendix D
### Corner Plots

A representative marginal posterior distribution of the CO line fitting at *r* = 0.6 is presented in Figure 12.

<figure><img src="image_15.png" /><figcaption>Figure 12. Corner plot of the CO line fitting at *r* = 0.6.</figcaption></figure>
## Appendix E
### The Case of the Gaussian Profile

Figure 13 indicates the results when we adopt the Gaussian function as an opacity profile instead of the Voigt function. There is a significant residual in the line wing region in contrast with the case of the Voigt function.

<figure><img src="image_16.png" /><figcaption>Figure 13. The same as Figure 2 but the Gaussian function is adopted for the opacity profile.</figcaption></figure>

## Appendix F
### Emitting Layer

Figure 14 shows the local specific intensity field as a function of (z, v) at r = 80 au in the best-fit model. While the

emission arises from the atmosphere at the line core ($v < 0.5$ km s$^{-1}$), the emission in the wing ($v > 0.5$ km s$^{-1}$) originates from the region near the midplane where most of the gas inhabits. This is because the wings are optically thin.

<figure><img src="image_17.png" /><figcaption>Figure 14. Local specific intensity field along the disk vertical axis in a unit of brightness temperature.</figcaption></figure>
Appendix G
Temperature Structure

By fitting the parametric temperature model, we recon-
structed the 2D temperature structure, as shown in Figure 15.
We plot the kinetic temperature traced by the optically thick
CO line center in the black solid line in Figure 15. It is shown
that the emitting layer of the CO line center is z/r ~ 0.4-0.6,
which is on the higher side among measured values in other

disks by channel map analysis and is roughly consistent (e.g.,
C. Pinte et al. 2018; C. J. Law et al. 2021, 2022). However,
since we rely on the parametric model, there may be some
uncertainty when comparing with more model-independent
measurements, even though the gas surface density measure-
ments do not depend on the vertical temperature structure as the
line wing emission is essentially optically thin, which is clear in
the corner plot (the upper-left panel in Figure 12).

<figure><img src="image_18.png" /><figcaption>**Figure 15.** Two-dimensional temperature structure constructed by using the best-fit parameters of the spectral model fitting.</figcaption></figure>
## Appendix H

### The C$^{18}$O J = 2–1 Line

The detection of pressure-broadened line wings in the $^{12}$CO line motivates us to check less optically thick isotopologue lines. As discussed in Section 4.1, an optically thin isotopologue line such as $^{13}$C$^{17}$O and $^{13}$C$^{18}$O can solve the degeneracy between the CO column density and the H$_2$ density. However, there is no available observation of them. Nonetheless, it would be worth confirming that more abundant isotopologues than $^{13}$C$^{17}$O and $^{13}$C$^{18}$O are detected and likely optically thick. Indeed, M. Galloway-Sprietsma et al. (2025) confirm that the $^{13}$CO J = 3–2 line is highly optically thick. The next less-abundant isotopologue is C$^{18}$O. R. Dong et al. (2017) and J. B. Bergner et al. (2019) observed the C$^{18}$O J = 2–1 line.

We downloaded the data sets published by R. Dong et al. (2017) and J. B. Bergner et al. (2019) from the ALMA Science Archive and ran pipeline calibration (IDs: #2013.1.01020.S, #2015.1.00964.S). We refer the readers to the original literature for details of each observation. We combined the measurement sets to achieve higher sensitivity. After subtracting the continuum emission, we created an image cube of the C$^{18}$O J = 2–1 line using the CASA command tclean with the robust parameter of 0.5 and a channel width of 0.2 km s$^{-1}$. We then smoothed the image cube with imsmooth and achieved a beam size of 0.26 × 0.26. The resultant rms noise level is 8.8 mJy beam$^{-1}$ in one channel. The integrated intensity and peak intensity maps are generated using the Python script better-moment (R. Teague & D. Foreman-Mackey 2018) by specifying the velocity channel of 4.0 < v < 5.2 km s$^{-1}$.

The integrated intensity map in Figure 16 shows ring-like emission, which is consistent with R. Dong et al. (2017) and J. B. Bergner et al. (2019). The integrated intensity near the ring is ~40 mJy beam$^{-1}$ km s$^{-1}$. In the right panel of Figure 16, we present the radial profile of the peak brightness

temperature using the full Planck function,

$$T_B' = \frac{h\nu}{k_B} \left[ \ln \left( \frac{2h\nu^3}{c^2 I_{C^{18}O}} + 1 \right) \right]^{-1}, \quad (H1)$$

where v is the frequency of the transition, and $I_{C^{18}O}$ is the peak intensity. The peak brightness temperature exceeds 30 K around the ring at r ~ 80 au and reaches ~35 K at the peak.

We found that the integrated intensity near the ring is exceptionally high as a C$^{18}$O J = 2–1 line among T Tauri disks when comparing similar radial locations (r ~ 80 au). Indeed, this has already been seen in J. B. Bergner et al. (2019). With their spatial resolution (~0.8), the peak integrated intensity around r ~ 80 au is ~180 mJy beam$^{-1}$ km s$^{-1}$. This is much higher than that in the other eight T Tauri disks in their sample but comparable with that in the Herbig disks around HD 163296 and MWC 480. In addition, the MAPS project observed the same line in five disks with a better spatial resolution of 0.15 (K. I. Öberg et al. 2021; K. Zhang et al. 2021). If we scale the integrated intensity in our map to the 0.15 resolution, we obtain ~13 mJy beam$^{-1}$ km s$^{-1}$. This value is again much higher than similar regions in T Tauri disks, IM Lup, GM Aur, and AS 209 by a factor of 2–6, but comparable with the Herbig disks (HD 163296 and MWC 480).

The peak brightness temperature of ~35 K suggests that the C$^{18}$O J = 2–1 line is highly optically thick assuming that the midplane temperature is ≤35 K. We note that the limited spatial resolution could underestimate the brightness temperature, implying that the actual brightness temperature is even higher. We can also compare it with the brightness temperature profiles of the same line in the MAPS sample (C. J. Law et al. 2021). Among these sources, only the very inner region (r < 50 au) of the Herbig disks (HD 163296 and MWC 480) show brightness temperatures higher than ~30 K.

<figure><img src="image_19.png" /><figcaption>Figure 16. (Left) Integrated intensity map of the C$^{18}$O J = 2–1 line. The white dotted circle indicates the dust ring location. The white solid circle shows the beam size. (Right) Radial profile of the peak brightness temperature. The black bar indicates the beam size.</figcaption></figure>
ORCID iDs

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

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

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

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

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

Munetake Momose https://orcid.org/0000-0002-3001-0897

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

References

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

Andrews, S. M., Wilner, D. J., Hughes, A. M., et al. 2012, ApJ, 744, 162

Arakawa, S., & Krijt, S. 2021, ApJ, 910, 130

Armitage, P. J., & Kley, W. 2019, From Protoplanetary Disks to Planet
Formation (Berlin: Springer)

Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., Günther, H. M.,
et al. 2018, AJ, 156, 123

Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A,
558, A33

Barenfeld, S. A., Carpenter, J. M., Ricci, L., & Isella, A. 2016, ApJ, 827, 142

Bergin, E. A., Cleeves, L. I., Gorti, U., et al. 2013, Natur, 493, 644

Bergner, J. B., Öberg, K. I., Bergin, E. A., et al. 2019, ApJ, 876, 25

Birnstiel, T., Dullemond, C. P., & Brauer, F. 2009, A&A, 503, L5

Birnstiel, T., Dullemond, C. P., & Brauer, F. 2010, A&A, 513, A79

Birnstiel, T., Dullemond, C. P., Zhu, Z., et al. 2018, ApJL, 869, L45

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

Carpenter, J. M., Ricci, L., & Isella, A. 2014, ApJ, 787, 42

Cleeves, L. I., Bergin, E. A., Bethell, T. J., et al. 2011, ApJL, 743, L2

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

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

Delussu, L., Birnstiel, T., Miotello, A., et al. 2024, A&A, 688, A81

De Zeeuw, P., Hoogerwerf, R. v., De Bruijne, J. H., Brown, A., & Blaauw, A.
1999, AJ, 117, 354

Doi, K., & Kataoka, A. 2023, ApJ, 957, 11

Dominik, C., & Dullemond, C. P. 2024, A&A, 682, A144

Dong, R., van der Marel, N., Hashimoto, J., et al. 2017, ApJ, 836, 201

Dullemond, C. P., Birnstiel, T., Huang, J., et al. 2018, ApJL, 869, L46

Dullemond, C. P., van Zadelhoff, G. J., & Natta, A. 2002, A&A, 389, 464

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

Flaherty, K. M., Hughes, A. M., Teague, R., et al. 2018, ApJ, 856, 117

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

Francis, L., Johnstone, D., Herczeg, G., Hunter, T. R., & Harsono, D. 2020, AJ,
160, 270

Gaia Collaboration, Vallenari, A., Brown, A. G. A., et al. 2023, A&A, 674, A1

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

Gordon, I. E., Rothman, L. S., Hargreaves, R. J., et al. 2022, JQSRT, 277,
107949

Guidi, G., Isella, A., Testi, L., et al. 2022, A&A, 664, A137

Hensley, B. S., & Draine, B. T. 2021, ApJ, 906, 73

Izquierdo, A. F., Stadler, J., Bae, J., et al. 2025, ApJL, 984, L8

Jiang, H., Macías, E., Guerra-Alvarado, O. M., & Carrasco-González, C. 2024,
A&A, 682, A32

Kanagawa, K. D., Muto, T., Tanaka, H., et al. 2015, ApJL, 806, L15

Keppler, M., Teague, R., Bae, J., et al. 2019, A&A, 625, A118

Krijt, S., Bosman, A. D., Zhang, K., et al. 2020, ApJ, 899, 134

Krijt, S., Schwarz, K. R., Bergin, E. A., & Ciesla, F. J. 2018, ApJ, 864, 78

Lankhaar, B., & Teague, R. 2023, A&A, 678, A17

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

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

Liu, H. B. 2019, ApJL, 877, L22

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

Loomis, R., Teague, R., Benisty, M., et al. 2025, ApJL, 984, L7

Luhman, K. L., & Mamajek, E. E. 2012, ApJ, 758, 31

Macías, E., Guerra-Alvarado, O., Carrasco-González, C., et al. 2021, A&A,
648, A33

Manara, C. F., Natta, A., Rosotti, G. P., et al. 2020, A&A, 639, A58

Mathews, G. S., Pinte, C., Duchêne, G., Williams, J. P., & Ménard, F. 2013,
A&A, 558, A66

Mayama, S., Akiyama, E., Panić, O., et al. 2018, ApJL, 868, L3

McMullin, J. P., Waters, B., Schiebel, D., Young, W., & Golap, K. 2007, in
ASP Conf. Ser. 376, Astronomical Data Analysis Software and Systems
XVI, ed. R. A. Shaw, F. Hill, & D. J. Bell (San Francisco, CA: ASP), 127

Minissale, M., Aikawa, Y., Bergin, E., et al. 2022, ESC, 6, 597

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

Miotello, A., van Dishoeck, E. F., Kama, M., & Bruderer, S. 2016, A&A,
594, A85

Miyake, K., & Nakagawa, Y. 1993, Icar, 106, 20

Musiolik, G., Teiser, J., Jankowski, T., & Wurm, G. 2016, ApJ, 827, 63

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

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

Okuzumi, S., & Tazaki, R. 2019, ApJ, 878, 132

Paneque-Carreño, T., Miotello, A., van Dishoeck, E. F., et al. 2023, A&A,
669, A126

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

Pinilla, P., Birnstiel, T., Ricci, L., et al. 2012, A&A, 538, A114

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

Preibisch, T., & Zinnecker, H. 1999, AJ, 117, 2381

Qi, C., Öberg, K. I., Espaillat, C. C., et al. 2019, ApJ, 882, 160

Qi, C., Öberg, K. I., Wilner, D. J., et al. 2013, Sci, 341, 630

Ricci, L., Testi, L., Natta, A., et al. 2010, A&A, 512, A15

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

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

Rybicki, G. B., & Lightman, A. P. 1991, Radiative Processes in Astrophysics
(New York: John Wiley & Sons)

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

Schwarz, K. R., Bergin, E. A., Cleeves, L. I., et al. 2016, ApJ, 823, 91

Sierra, A., Lizano, S., Macías, E., et al. 2019, ApJ, 876, 7

Sierra, A., Pérez, L. M., Zhang, K., et al. 2021, ApJS, 257, 14

Stadler, J., Benisty, M., Izquierdo, A., et al. 2023, A&A, 670, L1

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

Stammler, S. M., Dąrkowska, T., Birnstiel, T., et al. 2019, ApJL, 884, L5

Tan, Y., Skinner, F. M., Samuels, S., et al. 2022, ApJS, 262, 40

Teague, R. 2019, JOSS, 4, 1220

Teague, R., Bae, J., Bergin, E. A., Birnstiel, T., & Foreman-Mackey, D. 2018a,
ApJL, 860, L12
Teague, R., Bae, J., Birnstiel, T., & Bergin, E. A. 2018b, ApJ, 868, 113

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

Teague, R., & Foreman-Mackey, D. 2018, RNAAS, 2, 173

Toomre, A. 1964, ApJ, 139, 1217

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

Ueda, T., Flock, M., & Okuzumi, S. 2019, ApJ, 871, 10

Ueda, T., Kataoka, A., & Tsukagoshi, T. 2020, ApJ, 893, 125

Ueda, T., Tazaki, R., Okuzumi, S., Flock, M., & Sudarshan, P. 2024, NatAs, 8, 1148

van der Marel, N. 2023, EPJP, 138, 225

van der Marel, N., Birnstiel, T., Garufi, A., et al. 2021, AJ, 161, 33

Villenave, M., Stapelfeldt, K. R., Duchêne, G., et al. 2022, ApJ, 930, 11

Yang, B., Stancil, P. C., Balakrishnan, N., & Forrey, R. C. 2010, ApJ, 718, 1062

Yoshida, T. C., Nomura, H., Tsukagoshi, T., et al. 2025, ApJ, 980, 50

Yoshida, T. C., Nomura, H., Tsukagoshi, T., Furuya, K., & Ueda, T. 2022, ApJL, 937, L14

Youdin, A. N., & Goodman, J. 2005, ApJ, 620, 459

Yu, H., Teague, R., Bae, J., & Öberg, K. 2021, ApJL, 920, L33

Zagaria, F., Clarke, C. J., Booth, R. A., Facchini, S., & Rosotti, G. P. 2023, ApJL, 959, L15

Zhang, K., Bergin, E. A., Blake, G. A., Cleeves, L. I., & Schwarz, K. R. 2017, NatAs, 1, 0130

Zhang, K., Bergin, E. A., Schwarz, K., Krijt, S., & Ciesla, F. 2019, ApJ, 883, 98

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

Zhang, K., Bosman, A. D., & Bergin, E. A. 2020a, ApJL, 891, L16

Zhang, K., Isella, A., Carpenter, J. M., & Blake, G. A. 2014, ApJ, 791, 42

Zhang, K., Schwarz, K. R., & Bergin, E. A. 2020b, ApJL, 891, L17

Zhu, Z., Zhang, S., Jiang, Y.-F., et al. 2019, ApJL, 877, L18