exoALMA. XVI. Predicting Signatures of Large-scale Turbulence in Protoplanetary Disks

Marcelo Barraza-Alfaro, Mario Flock, William Béthune, Richard Teague, Jaehan Bae, Myriam Benisty, Gianni Cataldi, Pietro Curone, Ian Czekala, Stefano Facchini, et al.

► To cite this version:

Marcelo Barraza-Alfaro, Mario Flock, William Béthune, Richard Teague, Jaehan Bae, et al. exoALMA. XVI. Predicting Signatures of Large-scale Turbulence in Protoplanetary Disks. The Astrophysical Journal Letters, 2025, 984 (1), pp.L21. 10.3847/2041-8213/adc42d . hal-05064351

HAL Id: hal-05064351

https://hal.science/hal-05064351v1

Submitted on 12 May 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. XVI. Predicting Signatures of Large-scale Turbulence in Protoplanetary Disks

Marcelo Barraza-Alfaro¹, Mario Flock², William Béthune³, Richard Teague¹, Jaehan Bae⁴, Myriam Benisty⁵,², Gianni Cataldi⁶, Pietro Curone⁷,⁸, Ian Czekala⁹, Stefano Facchini⁷, Daniele Fasano⁵, Misato Fukagawa⁶, Maria Galloway-Sprietsma⁴, Himanshi Garg¹⁰, Cassandra Hall¹¹,¹²,¹³, Jane Huang¹⁴, John D. Ilee¹⁵, Andrés F. Izquierdo⁴,¹⁶,¹⁷,<fnref n="27" />, Kazuhiro Kanagawa¹⁸, Eric W. Koch¹⁹, Geoffroy Lesur²⁰, Cristiano Longarini²¹,⁷, Ryan A. Loomis²², Ryuta Orihara¹⁸, Christophe Pinte²⁰,¹⁰, Daniel J. Price¹⁰, Giovanni Rosotti⁷, Jochen Stadler⁵, Gaylor Wafflard-Fernandez²⁰, Andrew J. Winter⁵,², Lisa Wölfer¹, Hsi-Wei Yen²³, Tomohiro C. Yoshida⁶,²⁴, and Brianna Zawadzki²⁵,²⁶

¹Department of Earth, Atmospheric, and Planetary Sciences, Massachusetts Institute of Technology, Cambridge, MA 02139, USA; mbarraza@mit.edu

²Max-Planck-Institute for Astronomy (MPIA), Königstuhl 17, 69117 Heidelberg, Germany

³DPHY, ONERA, Université Paris-Saclay, 91120, Palaiseau, France

⁴Department of Astronomy, University of Florida, Gainesville, FL 32611, USA

⁵Université Côte d'Azur, Observatoire de la Côte d'Azur, CNRS, Laboratoire Lagrange, France

⁶National Astronomical Observatory of Japan, 2-21-1 Osawa, Mitaka, Tokyo 181-8588, Japan

⁷Dipartimento di Fisica, Università degli Studi di Milano, Via Celoria 16, 20133 Milano, Italy

⁸Departamento de Astronomía, Universidad de Chile, Camino El Observatorio 1515, Las Condes, Santiago, Chile

⁹School of Physics & Astronomy, University of St. Andrews, North Haugh, St. Andrews KY16 9SS, UK

¹⁰School of Physics and Astronomy, Monash University, Clayton VIC 3800, Australia

¹¹Department of Physics and Astronomy, The University of Georgia, Athens, GA 30602, USA

¹²Center for Simulational Physics, The University of Georgia, Athens, GA 30602, USA

¹³Institute for Artificial Intelligence, The University of Georgia, Athens, GA 30602, USA

¹⁴Department of Astronomy, Columbia University, 538 W. 120th Street, Pupin Hall, New York, NY 10027, USA

¹⁵School of Physics and Astronomy, University of Leeds, Leeds, LS2 9JT, UK

¹⁶Leiden Observatory, Leiden University, P.O. Box 9513, NL-2300 RA Leiden, The Netherlands

¹⁷European Southern Observatory, Karl-Schwarzschild-Str. 2, D-85748 Garching bei München, Germany

¹⁸College of Science, Ibaraki University, 2-1-1 Bunkyo, Mito, Ibaraki 310-8512, Japan

¹⁹Center for Astrophysics | Harvard & Smithsonian, 60 Garden St., Cambridge, MA 02138, USA

²⁰Univ. Grenoble Alpes, CNRS, IPAG, 38000 Grenoble, France

²¹Institute of Astronomy, University of Cambridge, Madingley Road, CB3 0HA, Cambridge, UK

²²National Radio Astronomy Observatory, 520 Edgemont Rd., Charlottesville, VA 22903, USA

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

²⁴Department of Astronomical Science, The Graduate University for Advanced Studies, SOKENDAI, 2-21-1 Osawa, Mitaka, Tokyo 181-8588, Japan

²⁵Department of Astronomy, Van Vleck Observatory, Wesleyan University, 96 Foss Hill Drive, Middletown, CT 06459, USA

²⁶Department of Astronomy & Astrophysics, 525 Davey Laboratory, The Pennsylvania State University, University Park, PA 16802, USA
Received 2024 November 30; revised 2025 March 12; accepted 2025 March 14; published 2025 April 28

Abstract

Turbulent gas motions drive planet formation and protoplanetary disk evolution. However, empirical constraints on turbulence are scarce, halting our understanding of its nature. Resolving signatures of the large-scale perturbations driven by disk instabilities may reveal clues on the origin of turbulence in the outer regions of planet-forming disks. We aim to predict the observational signatures of such large-scale flows, as they would appear in high-resolution Atacama Large Millimeter/submillimeter Array observations of CO rotational lines, such as those conducted by the exoALMA Large Program. Post-processing 3D numerical simulations, we explored the observational signatures produced by three candidate (magneto)hydrodynamical instabilities to operate in the outer regions of protoplanetary disks: the vertical shear instability (VSI), the magnetorotational instability (MRI), and the gravitational instability (GI). We found that exoALMA-quality observations should capture signatures of the large-scale motions induced by these instabilities. Mainly, flows with ring, arc, and spiral morphologies are apparent in the residuals of synthetic velocity centroid maps. A qualitative comparison between our predictions and the perturbations recovered from exoALMA data suggests the presence of two laminar disks and a scarcity of ring- and arc-like VSI signatures within the sample. Spiral features produced by the MRI or the GI are still plausible in explaining observed disk perturbations. Supporting these scenarios requires further methodically comparing the predicted perturbations and the observed disks' complex dynamic structure.

Unified Astronomy Thesaurus concepts: Protoplanetary disks (1300); Planet formation (1241); Hydrodynamical simulations (767); Radiative transfer simulations (1967)

1. Introduction

Turbulence is a ubiquitous phenomenon in astrophysical fluids (A. Brandenburg & Å. Nordlund 2011). In protoplanetary disks, gas turbulence can transport angular momentum radially, which, together with angular momentum extraction by

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

<footnote n="27">²⁷ NASA Hubble Fellowship Program Sagan fellow.</footnote>
disk winds, is needed to explain the accretion of matter onto
young stars (N. J. Turner et al. 2014). In the planet formation
process, turbulence can play a major role by inducing mixing/
diffusion and regulating the development of substructures
(G. Lesur et al. 2023), which directly affects the growth and
evolution of dust particles (T. Birnstiel 2024) and the chemical
inventory of the disk within planet-forming regions. Therefore,
turbulence is critical in the initial stages of planet formation
(J. Drążkowska et al. 2023; A. Miotello et al. 2023) and will
potentially impact the composition of the forming planets.
However, the dominant mechanism driving turbulence in
protoplanetary disks and its role in disk evolution and planet
formation are still uncertain and under active debate.

Turbulent motions span a wide range of scales. In a
simplified picture, the large scales inject energy into the
system, which cascades down to smaller turbulent eddies (see,
e.g., W. Lyra & O. M. Umurhan 2019). In protoplanetary disks,
turbulent length scales larger than or comparable to the disk's
local pressure scale height (roughly 10 au at 100 au from the
star) could be resolvable with Atacama Large Millimeter/
submillimeter Array (ALMA) observations of molecular lines.
Conversely, smaller eddies within the microscales are unresol-
vable with ALMA and can take the form of line broadening
(J. B. Simon et al. 2015). Contemporary methods to directly
constrain turbulence levels in the disk's outer region have focused
on detecting turbulent line broadening in ALMA observations
of molecular lines. These measurements of "nonthermal" line
broadening have reported evidence of nonzero subsonic turbu-
lence in two disks (S. Guilloteau et al. 2012; K. Flaherty et al.
2020; K. Flaherty et al. 2024; T. Paneque-Carreño et al. 2024) and
have set upper limits indicating relatively weak turbulence in four
other disks (K. M. Flaherty et al. 2015, 2017, 2018; R. Teague
et al. 2016; K. Flaherty et al. 2020). However, these methods are
model dependent, subject to observational limitations and
degeneracies, requiring stringent constraints on the disk structure
(e.g., A. M. Hughes et al. 2011; K. Flaherty et al. 2020;
G. P. Rosotti 2023; C. Hardiman et al. 2025, in preparation).
Further, these constraints provide limited information about the
nature of turbulence in the probed outer disk regions.

High spatial resolution ALMA observations of molecular
lines have revealed velocity perturbations in protoplanetary
disks (Disk Dynamics Collaboration et al. 2020; C. Pinte et al.
2023). Such perturbations, observed as deviations from the
background disk rotation, have opened the possibility to hunt
for signatures of large-scale motions driven by turbulence
mechanisms (C. Hall et al. 2020; M. Barraza-Alfaro et al.
2021). Identifying signatures from disk instabilities has the
potential to unravel the origin of turbulence in the disk's outer
regions. Moreover, since disk instabilities require specific
conditions to develop and be sustained, their detection can
constrain the disk's physical properties. Nevertheless, finding
evidence of large-scale turbulent motions in disks has been
challenging, as it requires high spectral and spatial resolution.
Such quality has been delivered by the exoALMA survey
(R. Teague et al. 2025), providing a detailed view of the disk
kinematic structure of 15 protoplanetary disks. These new
observations allow us to hunt for signposts of turbulence in the
large-scale flow predicted by 3D numerical simulations of disk
instabilities.

Various instabilities are expected to be active and to drive
turbulence in the disk outer regions (G. Lesur et al. 2023).
Gravitational instabilities (GIs) are predicted to develop in

massive protoplanetary disks (A. Toomre 1964, K. Kratter &
G. Lodato 2016). The magnetorotational instability (MRI) can
operate in ionized regions of magnetized protoplanetary disks
(S. A. Balbus & J. F. Hawley 1991; J. F. Hawley & S. A. Balbus
1991), while purely hydrodynamic (HD) instabilities can
dominate in the weakly ionized regions (H. Klahr et al. 2018;
W. Lyra & O. M. Umurhan 2019). Among the hydro instabilities,
the vertical shear instability (VSI) can operate in fast-cooling
regions of the outer disk (V. Urpin & A. Brandenburg 1998;
R. P. Nelson et al. 2013; T. Pfeil & H. Klahr 2019; C. Cui &
X.-N. Bai 2020, 2022), while the convective overstability
(H. Klahr & A. Hubbard 2014; W. Lyra 2014) and zombie
vortex instability (P. S. Marcus et al. 2015, 2016) are unlikely to
be dominant in the bulk of the disk outward tens of au from the
star (M. G. Malygin et al. 2017; T. Pfeil & H. Klahr 2019). These
instabilities give rise to disk substructures that are distinct from
each other and from other substructure-forming processes, such
as planet-disk interactions (J. Bae et al. 2023). Therefore, large-
scale turbulent motions detected by molecular line emission of
protoplanetary disks may unravel the instability at play from its
characteristic observational features.

This work explores the kinematic signatures of large-scale
gas motions produced by instability candidates to drive
turbulence in the outer regions of planet-forming disks:
the VSI, MRI, and GI. We present our results in the context
of the exoALMA Large Program, addressing its potential
to unravel large-scale turbulence signatures. The paper is
organized as follows: We describe the methods for the
numerical simulations in Section 2, and we detail the methods
used for the radiative transfer predictions, mock observations,
and calculation of observables in Section 3. We show the
simulations' results in Section 4 and their post-processing
results in Section 5. We discuss the findings and caveats of our
work in Section 6, and we draw our conclusions in Section 7.

## 2. Methods for Numerical Simulations

The numerical simulations presented in this work were
conducted with upgraded versions of the grid-based code
PLUTO<fnref n="28" /> (A. Mignone et al. 2007; see details of upgrades in
M. Flock et al. 2015, 2020; W. Béthune et al. 2021). We post-
processed outputs from three 3D global HD (magnetohydro-
dynamic (MHD)) simulations, one representative of each of the
disk instabilities explored. Particularly, we analyzed an
individual simulation snapshot for each of the studied
instabilities taken after the disk had reached a quasi-steady
state, which we considered representative of their gas
dynamics.

### 2.1. VSI-unstable Disk Simulation

The disk gas dynamics induced by the VSI is explored with
outputs of the 3D global HD simulation with radiative energy
transport presented in M. Flock et al. (2020). The simulation
was conducted with PLUTO version 4.2, using the setup of the
numerical simulations performed in M. Flock et al. (2017). The
radiation HD (henceforth RHD) equations are solved using a
hybrid stellar irradiation and flux-limited diffusion method
(M. Flock et al. 2013; see further details in M. Flock et al.
2020, Section 2). The parameters were chosen to simulate a
disk around a typical T Tauri star with a stellar mass of

<footnote n="28">²⁸ https://plutocode.ph.unito.it/</footnote>
$M_\star = 0.5 M_\odot$. The simulation domain covers 20–100 au in the radial domain, 0.7 rad in colatitude ($\pi/2 \pm 0.35$ rad, with the midplane at $\theta = \pi/2$), and $2\pi$ rad in azimuth. The simulation of the nonlinear evolution of the VSI was run with a grid with 1024, 512, and 2044 cells in $r, \theta$, and $\phi$, respectively, with a mesh logarithmically spaced in the radial direction and uniformly spaced in colatitude and azimuth. However, in order to reduce the cost of the radiative transfer calculations, we post-processed outputs interpolated to a coarser grid, half the grid resolution of the original in each direction, that is, 512, 256, and 1022 cells in $r, \theta$, and $\phi$, respectively. The domain in colatitude covers approximately 7.3 pressure scale heights, with approximately 70 cells per pressure scale height. Buffer zones were applied in the innermost 2 au and outermost 3 au of the radial domain. For the boundary conditions, a modified outflow condition was used in the $r$ and $\theta$ directions. The condition prevents gas inflow at the domain boundaries. Additionally, at the boundaries in colatitude the gas density was extrapolated logarithmically into the ghost zones (see Section 2 in M. Flock et al. 2020). For a summary of the simulation parameters, see Table 1 (see also Table 1 in M. Flock et al. 2020).

## 2.2. MRI-unstable Disk Simulation

The post-processing of the gas dynamics of an MRI-turbulent protoplanetary disk is performed with outputs from the MHD simulation presented in M. Flock et al. (2015). In particular, the 3D global nonideal MHD stratified simulation labeled D2G_e-2 ran for a disk with a dust-to-gas mass ratio of $10^{-2}$. Only the 2D ($R, \theta$) disk temperature structure was obtained using radiative transfer calculations, and it remained fixed via a prescribed locally isothermal equation of state (details in Section 2 of M. Flock et al. 2015). The nonideal MHD equations were solved using the MHD module of PLUTO (A. Mignone et al. 2007) including ohmic diffusion. Our simulation neglects the Hall effect and ambipolar diffusion. Similarly to the ohmic diffusion, ambipolar diffusion is a dissipative effect and can further damp the MRI, especially in the disk's outer regions (X.-N. Bai 2015; J. B. Simon et al. 2018), and it can drive additional annular substructures (M. W. Kunz & G. Lesur 2013; J. B. Simon & P. J. Armitage 2014; C. Cui & X.-N. Bai 2021). Distinctively, the Hall effect plays a role in transporting magnetic flux in the disk (X.-N. Bai & J. M. Stone 2017) but is not dissipative and can lead to higher magnetic stress, driving accretion (G. Lesur et al. 2014; W. Béthune et al. 2017; X.-N. Bai & J. M. Stone 2017). Overall, the Hall effect typically stabilizes small scales against the MRI (e.g., S. A. Balbus & C. Terquem 2001) and also facilitates large-scale structure formation (M. W. Kunz & G. Lesur 2013; W. Béthune et al. 2016; L. Krapp et al. 2018). Therefore, the inclusion of ambipolar diffusion and the Hall effect will modify the morphology of disk flow structures, increasing the complexity of the disk kinematic structure, which can make it difficult to identify signatures of MRI turbulence via gas kinematics. Moreover, ambipolar diffusion can also weaken the MRI-induced velocity perturbations in the outer disk regions, which is not optimal for exploring the observability of MRI turbulence via ALMA CO kinematics. Nevertheless, the Hall effect is predicted to be relevant in the midplane region of the outer disk (X.-N. Bai 2015), which is challenging to probe with ALMA molecular line observations.

Similar to the RHD simulation, the MHD simulation parameters are set to simulate a typical T Tauri system, with a central star of $M_\star = 0.5 M_\odot$ and a disk mass of $0.085 M_\star$. The radial and azimuthal domains of the MHD simulation are equal to the set for the RHD simulation, with $0.72$ rad in colatitude. In this case, the domain in $\theta$ covers approximately 6.4 disk pressure scale heights, $3.2H$ at each disk hemisphere. The MHD simulation grid was set to $256 \times 128 \times 512$ in $r, \theta$, and $\phi$, respectively ($\approx 20$ cells per $H$). The boundary condition in radius and colatitude continues the gradient of the gas density, pressure, azimuthal velocity, and resistivity. For the velocities in the $r$ and $\theta$ directions, a zero gradient condition was used, with a linear damping of the normal velocity component in the case of gas moving toward the simulation domain. In addition, the gas density is prevented from increasing in the $\theta$ ghost cells (see Section 2 in M. Flock et al. 2015). Further details on the MHD simulation's parameters are presented in Table 1 in Appendix B (see also Table 1 in M. Flock et al. 2015).

## 2.3. GI-unstable Disk Simulation

To explore the case of a disk with active GI, we post-processed an output from a self-gravitating HD (SGHD) simulation presented in W. Béthune et al. (2021). In particular, we used the outputs of the model labeled M5B10 in W. Béthune et al. (2021), which corresponds to a disk with an initial disk-to-star mass ratio of 1/5 and a dimensionless constant cooling parameter $\beta$ of 10, setting the timescale for radiative cooling ($t_{\text{cool}} = \beta \Omega_*^{-1}$, with $\Omega_*$ the local Keplerian frequency; see Equation (5) in W. Béthune et al. 2021). The radial domain of the SGHD simulation extends from 5 to 160 au (1–32 code units of length, set to 5 au), covered with 518 cells following a logarithmic spacing. The domain in $\phi$ covers the full $2\pi$ rad in the azimuthal direction with 512 uniformly sampled cells. In colatitude, the domain is covered by 96 cells in total, extending from $\theta = \pi/2 - 0.35$ rad to $\theta = \pi/2 + 0.35$ rad, where the grid has a finer sampling in the 0.1 rad covering the midplane region ($\pi/2 \pm 0.05$ rad), meshed with 32 uniformly spaced cells. The regions at the disk's upper layers are covered by a coarser mesh, with grid spacings stretched in $\theta$ (see Section 2.5.1 in W. Béthune et al. 2021). Near the midplane the simulation has approximately nine cells per $H$. For the simulation boundaries in colatitude ($\theta$), symmetric and antisymmetric boundary conditions are set for the gas density and velocities, respectively. In the radial direction, a symmetric boundary condition is set for $\rho$ and $v_\phi/r$, while an antisymmetric condition is used for $v_\theta$ and $r^2v_r$. In addition, a small pressure is imposed in the ghost cells for all boundaries (see Section 2.5.4. in W. Béthune et al. 2021). The parameters used for the SGHD simulation are presented in Table 1 (see also Table 1 in W. Béthune et al. 2021).

## 2.4. Considerations for Post-processing

Turbulence levels are expected to be strongest at the uppermost layers in disks with active VSI or MRI (see, e.g., M. Flock et al. 2015; J. B. Simon et al. 2015; T. Pfeil & H. Klahr 2021). Therefore, we focus our analysis and discussion on the surface disk layers, which are the regions probed by $^{12}$CO (3–2) line emission observed by exoALMA (within ALMA Band 7), with typical $Z/R$ values of ~0.27 (M. Galloway-Sprietsma et al. 2025). Deeper layers of the disk (closer to the disk midplane), as traced by $^{13}$CO (3–2) and CS (7–6) exoALMA observations (with typical $Z/R$ values of
~0.17 and ~0.18, respectively; M. Galloway-Sprietsma et al. 2025), will trace regions of weaker turbulence, in particular for the VSI and the MRI cases (M. Flock et al. 2011; J. B. Simon et al. 2018; T. Pfeil & H. Klahr 2021). In addition, the signal-to-noise ratios from the CS (7–6) line observations are not sufficient to fully exploit the spectral resolution achieved by exoALMA (R. Teague et al. 2025), crucial to capture the subtle gas motions driven by disk instabilities.

Due to the inclusion of radiation HD in the RHD simulation (M. Flock et al. 2020) and the use of a disk temperature obtained from radiative transfer calculations in the MHD simulation (M. Flock et al. 2020), it is not straightforward to adapt the code unit of length to rescale the simulation domain to model a larger disk, more suited to compare with the exoALMA sources. In particular, the RHD and MHD simulations’ disks temperature structures are set to model the radial range between 20 and 100 au, while the majority of exoALMA disks are significantly larger than 100 au in size. Since the size of disk substructures in the gas increases with disk radius (see, e.g., J. Bae et al. 2023), proportional to the disk pressure scale height, structures at larger radius will be easier to resolve spatially. However, if our RHD and MHD simulations are rescaled to larger disk sizes, the size of substructures in the synthetic predictions would be underestimated. Therefore, to account for the difference in disk sizes between our models and the exoALMA disks when exploring the observability of kinematic structures (Section 5.2), we convolved the synthetic images by a smaller synthesized beam than that of the fiducial exoALMA images (see Section 3.2). Our convolved images have a spatial resolution such that the substructures we resolve in our predictions can be compared to the ones we observe in the outer regions of the exoALMA sources, allowing a fair comparison between synthetic observations and exoALMA data (see Section 5.2). Finally, to compare the same radial region between the three simulations, we only explore the region between 20 and 100 au of the GI simulation, truncating the innermost 15 au and outermost 60 au in the post-processing.

### 3. Methods for Radiative Transfer Models

We post-processed the numerical simulations’ outputs with the Monte Carlo radiative transfer code RADMC-3D<fnref n="29" /> (C. P. Dullemond et al. 2012) version 2.0. We build our RADMC-3D input files using scripts based in FARGO2R-ADMC3D<fnref n="30" /> (C. Baruteau et al. 2019, 2021) but built for outputs from PLUTO (see further details in M. Barraza-Alfaro et al. 2021, 2024). The adopted framework is analogous to the exoALMA benchmark, presented and validated in J. Bae et al. (2025).

#### 3.1. Radiative Transfer Setup

Our disk model considered only one hemisphere of the disk (the front half of the disk as oriented in the sky, also referred to as the upper part of the disk; for a schematic view, see Figure 2 in C. Pinte et al. 2019). We set the number density of ¹²CO to zero for the lower half of the disk, having only the contribution from the upper disk ¹²CO layer to the line emission in our predictions. This choice allows a more direct comparison with

model residuals obtained from observations, which, with a complex line fitting, separate the line emission contribution from the upper and lower sides (see, e.g., S. Casassus et al. 2021; A. F. Izquierdo et al. 2022; A. F. Izquierdo et al. 2025). Since we are only interested in studying the perturbations obtained from the upper layer of the disk, not including the lower side in our models avoids complications in the analysis that would come from the mixing of line emission from both upper and lower ¹²CO layers. We adopted a constant value throughout the upper half of the disk of 10⁻⁴ for the fraction of ¹²CO relative to H₂ (n(¹²CO)/n(H₂)), including a highly simplified CO photodissociation and freeze-out prescriptions. Relevant toward the midplane region of the disk, CO freeze-out was included by reducing the CO number density by a factor 10⁻⁵ in regions where the gas temperature falls below 21 K, assuming that CO is then frozen onto dust grains (e.g., K. R. Schwarz et al. 2016; J. Bae et al. 2021). Important at the disk surface, photodissociation of ¹²CO was implemented by reducing its number density n(¹²CO) by a factor of 10⁻⁵ above the disk heights, where the vertically integrated gas column density, computed from the surface toward the midplane, has a value lower than 10²¹ cm⁻² (e.g., R. Visser et al. 2009; J. B. Simon et al. 2015; K. M. Flaherty et al. 2015; J. Bae et al. 2021; C. Baruteau et al. 2021). The freeze-out prescription allows us to recover a geometrically thin layer of ¹²CO (3–2) emission, while photodissociation helps us predict the emission layer at a height similar to empirical constraints, as obtained from resolved ¹²CO (3–2) observations of protoplanetary disks (typical Z/R ~ 0.27 in exoALMA disks; M. Galloway-Sprietsma et al. 2025; see also, e.g., C. J. Law et al. 2023; T. Paneque-Carreño et al. 2023). In addition, including photodissociation avoids possible artifacts from the edge of the simulation domain in colatitude in the line ray-tracing.

For simplicity and to have the same temperature structure for all disk models, we assumed a parametric temperature profile that takes into account radial and vertical gradients (K. Flaherty et al. 2020; see also, e.g., E. Dartois et al. 2003; K. A. Rosenfeld et al. 2013; K. M. Flaherty et al. 2015; M. Galloway-Sprietsma et al. 2025):

$$T_{\text{mid}} = T_{\text{mid0}} \left( \frac{R}{150 \text{ au}} \right)^q, \qquad (1)$$

$$T_{\text{atm}} = T_{\text{atm0}} \left( \frac{R}{150 \text{ au}} \right)^q, \qquad (2)$$

$$T_{\text{gas}} = \begin{cases} T_{\text{atm}} + (T_{\text{mid}} - T_{\text{atm}}) \left[ \cos \frac{\pi Z}{2Z_q} \right]^2, & \text{if } Z < Z_q \\ T_{\text{atm}}, & \text{otherwise} \end{cases} \qquad (3)$$

$$Z_q = Z_{q0} \left( \frac{R}{150 \text{ au}} \right)^{1.3}, \qquad (4)$$

where R and Z are the disk cylindrical disk coordinates, corresponding to the disk radius and vertical height, respectively. Z_q denotes the height from the disk midplane where the gas temperature reaches its maximum, set by T_atm. We assumed the parameters of the fiducial model for DM Tau presented in K. Flaherty et al. (2020), that is, T_atm0 = 24.68 K, T_mido = 14.3 K, q = -0.371, and Z_q0 = 70 au (see also M. Galloway-Sprietsma et al. 2025). However, the selection of our model temperature structure and parameters is not to reproduce a specific target but to consider a disk structure that

<footnote n="29">²⁹ https://www.itau.uni-heidelberg.de/dullemond/software/radmc-3d/</footnote>

<footnote n="30">³⁰ https://github.com/charango/fargo2radmc3d</footnote>
follows empirical constraints from ALMA observations of molecular lines. Two-dimensional views of our disk models' gas temperature and ¹²CO number densities are shown in Figure C4 in Appendix C.

For the line ray-tracing of the synthetic images, we assumed a distance to the source of 140 pc (approximately the average distance of the nearest 13 exoALMA sources; R. Teague et al. 2025). We assumed a disk position angle (PA) of 90° east of north, where the south part of the disk is nearest to the observer and the disk rotation is counterclockwise. We computed images for five disk inclinations ranging from $i = 5°$ to $i = 60°$, in order to cover the range of inclinations of the exoALMA disks (R. Teague et al. 2025). Finally, to predict the raw (unconvolved) kinematic signatures from large-scale turbulence (presented in Section 5.1), we considered the case of achieving ALMA's highest spectral resolution available in Band 7 at the frequency of the ¹²CO $J = 3-2$ rotational transition (≈345.796 GHz). Therefore, we produced data cubes with a 26 m s⁻¹ resolution (shown in the left 2.5 panels of Figure 2 and in Figure 3).

## 3.2. Synthetic ExoALMA Observations

For the case of predictions at exoALMA resolution (presented in Section 5.2), we binned the raw data cubes in frequency down to a velocity resolution of 104 m s⁻¹ using the CASA³¹imrebin task (J. P. McMullin et al. 2007), matching approximately the velocity resolution of the fiducial exoALMA data cubes for ¹²CO (3–2) (R. A. Loomis et al. 2025; R. Teague et al. 2025).

We explored the effect of the spatial resolution of ALMA observations by convolving our raw radiative transfer synthetic images with a circular Gaussian beam. Due to the difference in the disk sizes between our models and the observed disk sizes, we consider a Gaussian beam with an FWHM of 35 mas, while the exoALMA fiducial images have a ~150 mas spatial resolution (R. Teague et al. 2025). Our disk models have a radial extent of 100 au, while the average value of the flux-based disk radius enclosing 95% of the ¹²CO (3–2) emission, $R_{12CO}^{95}$, is 411 au for the exoALMA sample (M. Galloway-Sprietsma et al. 2025). Therefore, with our predictions' spatial resolution, the synthetic observations maintain a similar number of independent beams across the disk as the fiducial exoALMA images for the average exoALMA disk size (roughly ~40 beams across the disk diameter). This results in a comparable amount of resolved disk substructures to those present in the exoALMA observations. Before convolving the images with the Gaussian beam, we include noise in our models as white noise, resulting in an rms of ~0.16 mJy beam⁻¹ after convolution. The target rms value is set to be equivalent to the achieved by the fiducial images of the exoALMA survey ¹²CO (3–2) observations (~3 mJy beam⁻¹ for a larger beam of 150 mas and 100 m s⁻¹ velocity resolution; R. Teague et al. 2025).

## 3.3. Line Analysis Tools

To compute the observables from the raw and convolved synthetic images, we used the Python package BETTERMOMENTS<fnref n="32" /> (R. Teague & D. Foreman-Mackey 2018). We applied the percentiles method from BETTERMOMENTS to collapse each data cube and obtain the line centroid and line

width at every image pixel. The percentile method collapses the cube by taking the percentiles of the intensity-weighted velocity, for which the median returns the line centroid and the half of the range between the 16th and 84th percentiles returns an accurate line width measurement. To obtain the ¹²CO temperature, we computed the line peak maximum at each image pixel.

We performed the analysis of the observables obtained from the synthetic convolved images using the Python packages EDDY<fnref n="33" /> (Extracting Disk Dynamics; R. Teague 2019a) and GOFISH<fnref n="34" /> (R. Teague 2019b). We used EDDY to obtain the best-fit modified Keplerian disk model for an elevated surface to our synthetic ¹²CO (3–2) line centroid maps. The modified Keplerian disk rotation profile includes an additional radial-dependent component that can account for super-Keplerian or sub-Keplerian deviations from a purely Keplerian disk (see details in R. Teague et al. 2022; also included in Appendix A). The flexible rotation profile helps to capture the contributions from the gravitational potential of the disk and the pressure support in the disk's outer edge (see Figure 14 in R. Teague et al. 2022), giving an improved description of the background bulk rotation of the disk. The modified Keplerian model is particularly important in properly describing the super-Keplerian background rotation of the self-gravitating disk simulation.

To obtain our best-fit disk models, we computed the posterior distribution of our free parameters with eddy (see details in Appendix A). By subtracting the best-fit model of the projected disk bulk rotation from the perturbed disk line centroid map, the observable kinematic substructures driven by the disk instabilities can be explored in the residuals (presented in Section 5).

To obtain the line width and temperature maps residuals, we used GOFISH to subtract the disk azimuthal average in the ¹²CO (3–2) line peak (temperature) and line width maps, which takes into account the elevated emission surface used as input into the EDDY fitting. These residual maps retrieve the nonaxisymmetric substructures present in the disk, induced by deviations in the traced temperature and line broadening relative to the disk background.

# 4. Results from Numerical Simulations

In the following, we present the outputs of the set of numerical simulations of turbulent protoplanetary disks post-processed with radiative transfer calculations (see Section 3). We highlight the distinct morphology of the perturbations driven by the VSI, MRI, and GI in the gas density and velocities.

Figure 1 shows the perturbations driven by the VSI, MRI, and GI at 2.5 pressure scale heights ($H$) from the disk midplane. The layer shown is chosen as a representative snapshot of the gas dynamics at the upper layers, as traced by ¹²CO (3–2). We present the normalized gas density relative to a radial power-law fit to its azimuthal average in logarithmic scale ($\log_{10}(\rho/(\rho)^{\text{fit}}_{\phi})$), gas radial velocity ($v_r$), meridional (or vertical) velocity ($v_{\theta}$), and azimuthal velocity deviations relative to a radial power-law fit to its azimuthal average ($v_{\phi} - (v_{\phi})^{\text{fit}}_{\phi}$). The fitting of a radial power law to the azimuthal average of the azimuthal velocity, for which we use scipy.

<footnote n="31">³¹ https://casa.nrao.edu/</footnote>

<footnote n="32">³² https://github.com/richteague/bettermoments</footnote>

<footnote n="33">³³ https://github.com/richteague/eddy</footnote>

<footnote n="34">³⁴ https://github.com/richteague/gofish</footnote>
`optimize.curve_fit`, allows a superb subtraction of the disk bulk rotation, extracting the perturbations in $v_φ$. We also show the perturbations at the disk midplane in Figure C1 in Appendix C for completeness.

In the VSI snapshots (first row of Figure 1), we display the characteristic dynamics from the nonlinear state of a VSI-unstable disk, in which the instability has reached a saturated state. In the gas density, VSI drives subtle perturbations, where the maximum density contrast within annuli at $Z = 2.5H$ has a median value of ≈15% (computed within 35 and 85 au). Such weak density perturbations are likely not relevant in observations of line emission of protoplanetary disks. However, in the gas velocity the VSI drives substantial perturbations in the three velocity components. In particular, the VSI induces a characteristic corrugated flow pattern in the meridional direction, which displays a high degree of axisymmetry (see, e.g., L. Flores-Rivera et al. 2020; M. Barraza-Alfaro et al. 2021). These velocity perturbations have an average rms velocity at $Z \sim 2.5H$ of ~1% the local Keplerian velocity. In CO kinematic observations, the apparent difference in axisymmetry of the velocity perturbations among velocity components becomes relevant when exploring their projection into the line of sight, in which inclined disks will have a larger contribution from the radial and azimuthal velocities (see Section 5).

For the MRI simulation (middle row of Figure 1), the snapshots also reveal that the MRI-driven gas density structure is dominated by a radial gap and overdensity induced by the presence of an outer MRI dead-zone edge, a radial transition in the disk turbulence levels (M. Flock et al. 2017). A steep transition in turbulence levels develops in the simulation owing to the inclusion of ohmic diffusion. However, the location of the MRI dead-zone outer edge and its dynamics can be substantially affected if the Hall effect and ambipolar diffusion are also taken into account in the nonideal MHD simulation (e.g., G. Lesur et al. 2014; see also Sections 2 and 6.3). In particular, the Hall effect could “revive” regions of the dead zone (M. W. Kunz & G. Lesur 2013; G. Lesur et al. 2014).

Azimuthally asymmetric perturbations in the gas density at $Z \sim 2.5H$ remain below 65% relative to the background. In the velocities, the perturbations driven by the MRI are larger than VSI-driven ones, with an average rms velocity at $Z \sim 2.5H$ of ~3.6% the local Keplerian velocity, roughly 3.6 times larger than the velocity perturbations driven by the VSI. Distinctively, MRI-induced perturbations show radial flows that are substantially stronger than those in the meridional and azimuthal directions, where the latter two have similar magnitudes. At $Z \sim 2.5H$, the mean rms velocity of the radial perturbations of the MRI is approximately a factor 1.46 larger than the mean rms velocity of the ones in the meridional and azimuthal directions. In addition, its morphology shows a spiral structure, composed of multiple spiral arms (high-order mode spirals are usually dominant for the MRI; see Section 5.1.2 in J. Bae et al. 2023 and references therein). The MRI-driven spiral structure is particularly defined in the radial and meridional velocity components (see also Figure C2 for a polar deprojection of the simulation fields). In the azimuthal velocity, rings of super- and sub-Keplerian velocities are present, induced at the MRI dead-zone outer edge. Overall, the velocity structure driven by MRI shows a turbulent-like morphology, with localized flows and “mesoscale” eddies (with scales smaller than the gas pressure scale height) distorting the large-scale spiral structures and non-

Keplerian rings. These turbulent features might complicate the identification of MRI-driven large-scale coherent structures (e.g., with filament detection methods; E. W. Koch & E. W. Rosolowsky 2016; A. F. Izquierdo et al. 2022, 2023, 2025) and could lead to misidentifying localized deviations when searching for embedded planets' signposts, further discussed in Section 5.

GI-driven perturbations (bottom row of Figure 1) show a characteristic large-scale spiral structure, exhaustively studied in previous works (e.g., G. Lodato & W. K. M. Rice 2004; P. Cossins et al. 2009; D. Forgan et al. 2011; R. Dong et al. 2015; C. Hall et al. 2018; W. Béthune et al. 2021). These spiral structures are present in the three velocity components, with common pitch angles (see Figure C2). From Figure 1, we observe that the density perturbations induced by GI are substantially larger than those driven by the VSI or the MRI. Therefore, the GI-driven density perturbations have a higher chance to create substructures in the $^{12}$CO (3–2) line peak and line width, which we address in further detail in Section 5. Similarly, the velocity perturbations driven by the MRI and the GI are substantially stronger than those driven by the VSI. In our sample of simulations, the average rms velocities between 25 and 95 au at $Z \sim 2.5H$ are 29, 111, and $142 \text{ m s}^{-1}$ for the VSI, MRI, and GI, respectively. Therefore, the GI is more likely to create observable features in $^{12}$CO (3–2) kinematic observations. For completeness, we show the entire domain of the SGHD simulation in Figure C3 in Appendix C.

Differences in magnitude between the velocity perturbations at the disk midplane and upper layers for the MRI and the VSI cases are evident when comparing the Cartesian view of the perturbations at $Z \sim 2.5H$, shown in Figure 1, and the perturbations at $Z \sim 0$, provided in Figure C1 in Appendix C. However, in the case of perturbations driven by GI, the radial and azimuthal motions are still vigorous near the midplane, showing a decrease in magnitude toward the midplane for the meridional velocity component only, consistent with previous findings (J.-M. Shi & E. Chiang 2014). In our sample of simulations, the average rms velocity moving down from $Z \sim 2.5H$ to $Z \sim 1.5H$ shows a decrease of ~47% for the VSI, ~62% for the MRI, but only ~18% for the GI. The weak dependency of the GI-induced rms turbulent velocities with disk height aligns with previous studies (e.g., D. Forgan et al. 2012; J.-M. Shi & E. Chiang 2014). Therefore, mapping deeper layers with molecular line observations should still reveal GI perturbations (see also C. Longarini et al. 2021) but is less likely to reveal VSI- or MRI-driven large-scale motions.

Finally, we stress that our simulations of unstable disks are for particular sets of parameters, which control the resulting strength of the instabilities. In conditions closer to the ideal MHD limit, the MRI can be stronger (e.g., M. Flock et al. 2015). Conversely, ambipolar diffusion can further decrease the strength of the MRI (e.g., J. B. Simon et al. 2015; see also Section 2). In the case of GI, the strength of the perturbations tightly depends on the disk mass and its cooling properties (W. Béthune et al. 2021). For the VSI, a steep disk thermal stratification enhances the vertical shear, driving stronger VSI turbulence (H.-G. Yun et al. 2025a, 2025b). Therefore, our observability study focuses on the differences in morphology of the kinematic signatures rather than their strength. Nevertheless, we comment on the variations in magnitude of the velocity perturbations and their impact on the observable signatures. Still, the explored observable kinematic features
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Cartesian view of the disk perturbations at $Z \sim 2.5H$ for the HD (top), MHD (middle), and SGHD (bottom) simulations, with active VSI, MRI, and GI, respectively. Shown are, from left to right, perturbations in the gas density, radial velocity, meridional velocity, and azimuthal velocity. The color bar limits are adjusted for each instability.</figcaption></figure>

must be taken as a reference from a particular set of simulations, for which the post-processed snapshots were considered to be representative of the overall dynamics of the studied disk instabilities.

# 5. Results from Radiative Transfer Models

## 5.1. Kinematic Signatures in the Raw Images

In the left 2.5 panels of Figure 2, we show a set of raw channel maps computed from our radiative transfer models for a disk inclination of 30°. From top to bottom, we show the raw channels of disks unstable to the VSI, MRI, and GI. In all channels, perturbations are visible as spurs or wiggles at the edges of the emission areas, driven by the large-scale velocity perturbations. While in our models there are apparent differences between the signatures of the MRI and the GI in the channel maps, these signatures alone may not be sufficient to robustly distinguish their origin. In particular, the differences in the channel maps may become less distinguishable in the

case of disks with stronger MRI or weaker GI than in our simulations. For mid-inclination disks, although both instabilities, the MRI and the GI, show clear perturbations in all channel maps, the vigorous radial gas motions of the GI create substantial perturbations in the channels mapping near the systemic velocity, for which the GI produces its most studied signature, the “GI wiggle” along the disk semiminor axis (C. Hall et al. 2020; C. Longarini et al. 2021; J. P. Terry et al. 2022; J. Speedie et al. 2024), whose quantitative analysis may hint at the GI origin of the signatures. However, in disks with low inclination, for which the radial velocity does not significantly contribute to the line-of-sight velocity, both the MRI and GI may produce similar features in the velocity channel maps. In any case, robust identification of the origin instability from the channel maps would require a framework to quantitatively study their properties across the full $^{12}$CO line. Currently, the identification of the origin of the observed perturbations in the channel maps requires dedicated numerical simulations and radiative transfer predictions, as in the case of
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Selected channel maps from ¹²CO (3–2) predictions of disks unstable to the VSI, MRI, and GI, for a disk inclination of 30°. The left 2.5 panels show the raw product, while the right 2.5 panels include the effect of the spatial resolution (35 mas FWHM Gaussian beam), spectral resolution (104 m s⁻¹ channel spacing), and sensitivity of exoALMA. The models only consider the upper (or front) CO layer.</figcaption></figure>

massive planets driving a localized signature (see, e.g., C. Pinte et al. 2023) or GI signatures (T. Paneque-Carreño et al. 2021; J. Speedie et al. 2024).

In Figure 3, we show the perturbations in the line centroid maps, predicted for our raw synthetic images assuming a perfect background rotation subtraction. To obtain these ideal residuals, we generated a radiative transfer model in which we subtract the disk rotation in the input azimuthal velocity. In order to achieve a superb background rotation subtraction, we fit a power-law profile to the azimuthal average of the azimuthal velocity at each radius and height (r and θ; see also Section 2). The fitted 2D background rotation of the disk is then subtracted from the original azimuthal velocity field. The line centroid velocity map of such a model directly gives the disk velocity residuals relative to its bulk rotation, showing the expected signatures from the perturbations driven by the disk instabilities. From top to bottom in Figure 3, we show the residuals from our VSI-unstable, MRI-unstable, and GI-unstable disk models. From left to right, we show the effect of increasing disk inclination from 5° to 60°, approximately covering the range of inclinations of the exoALMA sample of protoplanetary disks (see Table 1 in R. Teague et al. 2025).

For low to moderate disk inclinations ($5° \le i \le 30°$), we recover the quasi-axisymmetric rings and arcs in the residual maps from the VSI, predicted from the presence of strong meridional motions (M. Barraza-Alfaro et al. 2021), while for the MRI and the GI spiral-like structures dominate the line centroid residuals (see also C. Hall et al. 2020). For disk inclinations larger than 30°, the characteristic morphology of the residuals from the VSI is lost, and spiral-like arc structures are seen. The lower degree of azimuthal symmetry of VSI-induced non-Keplerian residuals for higher disk inclinations is

a result of the larger contributions of the radial and azimuthal velocity components to the projected line-of-sight velocity. In particular, coherent radial velocity flows will have their strongest contribution to the line of sight along the semiminor axis. In contrast, the contribution will be most substantial along the semimajor axis for the azimuthal velocity component (see R. Teague et al. 2019). The distinct radial and azimuthal velocity structure induced by the VSI (see Figure 1), combined with the modulation of their contribution to the line-of-sight velocity as function of the disk PA, and partially mixed with the azimuthally extended meridional flows, gives rise to the observed pattern in the non-Keplerian residuals of the mid-inclination VSI-unstable disk. Therefore, the observed pattern in the residuals of the inclined VSI-unstable disk is mostly produced by the varying modulations of the radial and velocity components into the line of sight for different disk PAs.

In the case of the MRI, higher disk inclinations strengthen the contribution from the sub- and super-Keplerian rings induced at the MRI dead-zone outer edge (see Section 4). In addition, nonaxisymmetric structures seen as localized perturbations are apparent in the residuals of the mid-inclination disks ($30° \le i \le 60°$). Due to the projection effects, the azimuthal velocity deviations induced by non-Keplerian rings are more prominent along the disk’s major axis, with opposite signs at the east and west sides of the disk. Similarly, the spiral-like radial velocity flows produced by the MRI are better traced for mid-inclination disks, inducing filamentary substructures. Still, the radial and azimuthal velocity perturbations are mixed when projected into the line of sight, in which their contributions vary depending on the disk PA, resulting in intricate non-Keplerian velocity residuals.

For the GI, the spiral velocity perturbations have similar morphologies for all disk inclinations, due to the predominance
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Line centroid map residuals of the disk velocity relative to its rotation from raw ¹²CO (3–2) predictions of disks unstable to the VSI, MRI, and GI, assuming a perfect background subtraction of the disk rotation. The raw images have a 26 m s⁻¹ velocity resolution. From left to right, different disk inclinations are shown. The models consider only the upper (or front) disk hemisphere, and the disk rotation is counterclockwise in the sky. The color bar limits have been adjusted for each instability.</figcaption></figure>

of spiral structures in all three of the gas velocity components. Interestingly, fewer spiral arms are apparent for larger disk inclinations, resulting from the minor contribution to the line-of-sight velocity of the spatially narrower meridional spiral flows. Additionally, the residuals' strength is fairly consistent for all disk inclinations along the northern part of the disk minor axis, which can be attributed to the strong radial velocity of the spirals driven by the GI (C. Hall et al. 2020; C. Longarini et al. 2021; J. Speedie et al. 2024). For the disk with an inclination of 60°, we attribute to projection effects along the line of sight the observed stronger residuals along the northern disk minor axis relative to the southern side of the disk, mainly due to the inclined viewing angle combined with the elevation of the emission surface being traced.

## 5.2. Signatures in the Synthetic Observations

In the right 2.5 panels of Figure 2, we show the effect of exoALMA observation resolution in the ¹²CO (3–2) channel maps (see details in Section 3.2), computed from our radiative transfer models for a disk inclination of 30°. Due to the effect of the synthetic observation resolution, the spurs and wiggles induced by the instabilities are blurred. In particular, the perturbations driven by the VSI and the MRI are challenging to identify, which may limit its observability in a direct inspection of the channel maps (e.g., C. Pinte et al. 2025).

In Figure 4, we show the kinematic signatures retrieved from synthetic observations with a spatial and spectral resolution comparable to that of the exoALMA survey. Following an observer's approach, the gas velocity perturbations are extracted by subtracting the line-of-sight velocity map of a smooth

modified Keplerian disk model from the line-of-sight velocity map obtained from the unstable disks' synthetic observations (see details in Section 3.3). The resulting residuals unravel the kinematic signatures observable in exoALMA quality data. Relative to the unconvolved ideal case, presented in Figure 3, the convolution by the Gaussian beam blurs the finer substructures driven by the instabilities, due to their spatial scales being smaller than the angular resolution (35 mas). Still, large-scale ring-like, arc-like, and spiral-like substructures are detectable and well resolved with the exoALMA spatial resolution. Thus, the superb spectral resolution achieved by exoALMA, equal to a velocity resolution of ~100 m s⁻¹ in the ¹²CO (3–2) fiducial images, is sufficient to clearly resolve in velocity the perturbations induced by the instabilities, capturing even subtle gas motions of a few tens of meters per second.

In the predictions of the VSI, subtle perturbations with a quasi-axisymmetric ring or an azimuthally extended arc-like morphology are recovered for low to intermediate inclinations. For higher disk inclinations ($i \ge 45°$), the residuals have a more significant degree of asymmetry, resulting in arc-like features with smaller azimuthal extent. Additionally, the near side of the disk with an inclination of $i = 60°$ shows substantial systematic errors arising from the EDDY model not accounting for radiative transfer effects in the line emission projected into the line of sight (also apparent in the MRI disk model with $i = 60°$). In particular, differences in the height of the probed ¹²CO layer between the north and south regions of the disk are not accounted for in the EDDY model.

For the MRI case, spiral-like structures are apparent for low disk inclinations ($i \le 15°$), while only filaments tracing spiral
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Line centroid map residuals of the disk velocity relative to its rotation from ¹²CO (3–2) predictions of turbulent protoplanetary disks, including the effects of the spatial resolution (35 mas FWHM Gaussian beam), the spectral resolution (104 m s⁻¹ channel spacing), and the sensitivity of exoALMA. We show the resulting velocity deviations from a modified Keplerian model for a geometrically thick disk, fitted to the line centroid map using EDDY. The disk models consider only the upper (or front) ¹²CO layer, and the disk rotation is counterclockwise in the sky. From top to bottom, we show unstable disks with active VSI, MRI, and GI. From left to right, different disk inclinations are shown. The beam size is displayed in the lower left corner of each panel. A 5σ clipping was applied to all moment maps. The color bar limits have been adjusted for each instability.</figcaption></figure>

segments are recovered in our mid-inclination models (30° ≤ i ≤ 60°). We observe that, in some cases, spiral-like substructures that are separated features in the raw predictions (Figure 3) are merged by the limited spatial resolution, appearing as one coherent large-scale structure, which can allow its observability within the capabilities of ALMA. For the disks with inclinations of i = 45° and i = 60°, the dominant residuals are the azimuthal velocity deviations induced at the edge of the dead zone, prominent along the semimajor axis of the disk. However, we stress that spatially resolved substructures such as the ones in our predictions are only resolvable in the exoALMA disks' outermost regions. Therefore, the MRI dead-zone outer edge should be located farther from the star in the exoALMA disks than in our model to produce resolvable features. An MRI dead-zone outer edge at a larger radius is a possibility, as its location can substantially vary depending on the systems' properties (T. N. Delage et al. 2022). Last, localized features become more apparent with increasing disk inclination, with a prominent “Doppler flip” resolved for the disk with i = 60°, which we discuss in Section 6.2.

The velocity perturbations recovered from our GI predictions are more prominent than those of the VSI or the MRI and, therefore, have a higher chance of detection within the exoALMA sample. The difference comes from the GI-driven spiral perturbations having a larger spatial scale and substantially higher velocities in our models. Our predictions indicate that large-scale spiral structures triggered by GI should be apparent if present in the sample, reaching magnitudes of ~100 m s⁻¹. These results align with estimations by analytical

models of GI-induced spirals (C. Longarini et al. 2021) and are consistent with findings of GI signatures in AB Aur (J. Speedie et al. 2024). For high disk inclinations (i ≥ 45°), our fitted disk rotation model to the GI-unstable disk requires the inclusion of a component that accounts for the super-Keplerian nature of massive self-gravitating disks (see Section 3.3). If the modification to the purely Keplerian rotation model is not included, a global dipole pattern dominates the line centroid map residuals. The need for this correction is consistent with other studies showing that, for fairly massive disks (with masses above 5% of the stellar mass), the effect of disk self-gravity in the disk rotation is detectable with high-resolution ALMA observations of molecular lines, allowing its kinematic mass measurement (B. Veronesi et al. 2021, 2024; S. M. Andrews et al. 2024; C. Longarini et al. 2025).

We also explored the nonaxisymmetric signatures in the line peak intensity (or brightness temperature) and width (see method details in Section 3.3). The residuals in temperature and line width are shown in Figures C5 and C6 in Appendix C, respectively. The GI is the only instability that drives significant nonaxisymmetric substructures in the temperature residuals. These deviations are likely due to changes in the elevation of the ¹²CO-emitting layer resulting from our simplified photodissociation prescription. In the interarm regions, photodissociation removes CO from regions deeper within the disk. Because of the assumption of a smooth axisymmetric temperature structure in our model, the variations in CO emission height are reflected in brightness temperature deviations, traced by the peak intensity of the optically thick
$^{12}$CO line. Still, similar features may be observable in GI-unstable protoplanetary disks and possibly amplified owing to GI spiral arm heating from shocks and compression (e.g., P. Cossins et al. 2009), ignored in our radiative transfer model. In the case of the MRI, previous numerical simulations have predicted that it can drive nonaxisymmetric temperature deviations (e.g., J. Ross & H. N. Latter 2018), ignored in our predictions owing to the locally isothermal equation of state used in our MHD simulation and the assumption of an axisymmetric disk temperature in the post-processing. However, MRI-driven temperature deviations found in previous works usually remain small (see, e.g., M. Flock et al. 2013) and are mainly due to local current sheets (see, e.g., Figure 9 in M. Flock et al. 2013). Moreover, temperature deviations driven by the MRI are subject to uncertainties related to radiative cooling times and efficiency. Models considering a temperature structure from a radiative MHD simulation of an MRI-active disk will be able to properly assess the detectability of MRI-induced temperature perturbations in ALMA CO rotational line emission observations.

Coherent nonaxisymmetric line width residuals with spiral-
like morphology are also prominent in the GI case, possibly
originating from a combined effect of temperature broadening,
turbulent broadening, and opacity broadening at overdensities
(A. Hacar et al. 2016). However, for the line width map residuals,
only the disk with an inclination of 5° is not substantially affected
by artifacts resulting from our procedure. The subtraction of an
azimuthal average does not consider projection effects in the line
width map, which requires a comprehensive disk model (e.g., as
provided by the Discminer modeling framework; A. F. Izquierdo
et al. 2021; A. F. Izquierdo et al. 2025). Our predictions of
observable GI-induced spiral signatures in the line width map are
in agreement with ALMA observations and supporting models
exhibiting evidence of GI acting in the AB Aurigae planet-
forming disk (J. Speedie et al. 2024).

While the observable nonaxisymmetric residuals induced by
the MRI or the large-scale spirals driven by the GI could
be confused by planetary signatures, the different scenarios can
be disentangled by using information of the peak and width of
the emission line, in combination with the use of multiple
molecular tracers, which we expand in Section 6.2. Overall, our
results for the VSI and the GI cases are consistent with previous
findings presented in M. Barraza-Alfaro et al. (2021) and
C. Hall et al. (2020), while we show for the first time that MRI-
driven large-scale velocity perturbations can be detectable with
ALMA CO kinematics.

6. Discussion

6.1. *Turbulence Signatures in ExoALMA Disks*

Observations of the kinematic structure of protoplanetary
disks from the exoALMA program have revealed that
signatures of disk perturbations are common (R. Teague
et al. 2025). These perturbations are present as “kinks” in
channel maps and coherent substructures in the line centroid,
line peak, and line width maps (e.g., C. Pinte et al. 2018;
R. Teague et al. 2022). In particular, signatures in channel
maps and residuals in the line centroid maps of $^{12}$CO (3–2) line
observations unravel the disk gas velocity perturbations in
detail, hinting at the physical processes occurring in the disk.

The nearest protoplanetary disk in the sample, V4046 Sgr,
shows no evidence of large-scale signatures in its outer regions

(R. Teague et al. 2025), displaying smooth channel maps
(C. Pinte et al. 2025). Similarly, smooth outer regions are
present in the PDS 66 disks (C. Pinte et al. 2025; R. Teague
et al. 2025). These results indicate that V4046 Sgr and PDS 66
have nearly laminar outer regions. Therefore, these regions
have the specific conditions not to develop or sustain any disk
instability: they have low surface densities such that they are
stable to the GI; they have a low ionization fraction such that
the gas is not susceptible to magnetic effects, that is, stable to
the MRI; and they have cooling timescales sufficiently large to
damp VSI motions due to the disk vertical buoyancy
(timescales comparable to or larger than the local dynamical
timescale $\Omega_{\text{Kep}}^{-1}$; R. P. Nelson et al. 2013; M.-K. Lin &
A. N. Youdin 2015). These results are consistent with previous
turbulence level upper limits, constraining weak turbulence in
the outer regions of V4046 Sgr (K. Flaherty et al. 2020).

The quasi-axisymmetric ring- and arc-like non-Keplerian
motions predicted from large-scale VSI-driven turbulence are
present in two disks from the exoALMA sample (see, e.g.,
A. F. Izquierdo et al. 2025). In particular, LkCa 15 shows arcs
of non-Keplerian motions in the line centroid map residuals
(A. F. Izquierdo et al. 2025) and azimuthal velocity deviations
in the rotation curves (A. F. Izquierdo et al. 2025; J. Stadler
et al. 2025). However, its inclination of ~50° makes a
challenging case to robustly identify VSI motions owing to the
minor contribution of the meridional velocity to the line-of-
sight velocity (see Section 5). Similarly, the outermost regions
of RXJ 1604 show one ring-like structure in its projected
velocity residuals (J. Stadler et al. 2023). However, its apparent
connection to a tightly wound spiral might suggest a planetary
origin (J. Stadler et al. 2023). In addition, predictions of VSI
motions result in multiple concentric rings or arcs of Doppler-
shifted velocities, not a single ring-like substructure. The
overall lack of VSI motions in the exoALMA disks could
indicate that most disks have relatively long cooling timescales
in their outer regions, which halts the development of the VSI.
Alternatively, it could suggest that the gas dynamics of the disk
in the sample is dominated by other disk instabilities (e.g., the
MRI or the GI) or planet-disk interactions.

Spiral structures are present in several exoALMA disks,
including CQ Tau, HD 135344B, and MWC 758 (L. Wölfer
et al. 2021; A. F. Izquierdo et al. 2025; T. E. Hilder et al. 2025;
L. Wölfer et al. 2025). GI- or MRI-driven motions could be
responsible for some of the observed spiral-like velocity
residuals. In particular, GI-induced substructures could also
explain some of the observed spiral features in temperature and
line width (A. F. Izquierdo et al. 2025; L. Wölfer et al. 2025).
However, the strong nonaxisymmetric features or low disk
inclinations impede the accurate measurement of the kinematic
disk mass of the targets with the most prominent large-scale
spirals (including MWC 758, CQ Tau, HD 135344B, HD
143006, and J1604; C. Longarini et al. 2025). Therefore, their
stability against GI is still an open question. The subset of disks
with kinematic mass measurements is GI stable (C. Longarini
et al. 2025; also validated by an alternative method in L. Trapman
et al. 2025). Therefore, only the MRI could still be active within
that subset of disks, driving some of the filamentary spiral-like
structures observed.

The only disk in the exoALMA sample with previous
measurements of nonzero turbulence, DM Tau, shows evidence
of large-scale flows in the channel maps (e.g., C. Pinte et al.
2025; C. Hardiman et al. 2025, in preparation). Its disk mass
<figure><img src="image_8.png" /><figcaption>**Figure 5.** $^{12}$CO (3–2) line centroid map residuals and selected channel maps from synthetic predictions of an MRI-unstable disk for a disk with an inclination of 60°, highlighting a localized feature. The images include the effect of the spatial resolution (35 mas FWHM Gaussian beam, shown in the lower left corner), spectral resolution (104 m s$^{-1}$ channel spacing), and sensitivity equivalent to that of exoALMA.</figcaption></figure>

measured via kinematics indicates that the DM Tau is GI stable (C. Longarini et al. 2025), indicating the possibility of active MRI, previously proposed as the source of the observed nonthermal broadening (K. Flaherty et al. 2020). Nevertheless, efforts in constraining the level of turbulence in DM Tau via line broadening from the exoALMA data have shown that these measurements are still degenerated even for high-quality data (C. Hardiman et al. 2025, in preparation). Particularly, the level of microscale turbulence needed to fit the data is highly sensitive to the model temperature structure and the prescription used for the parameter controlling the turbulent broadening (C. Hardiman et al. 2025, in preparation; see also K. M. Flaherty et al. 2015). Therefore, new methods to measure turbulence that lift these degeneracies and can account for resolved substructures are needed to robustly measure the turbulent broadening in DM Tau (C. Hardiman et al. 2025, in preparation). Further, analysis of the residual maps and models of higher complexity are required to obtain conclusions regarding the nature of turbulence in DM Tau, possibly confirming the MRI origin of its large-scale perturbations.

The disks around HD 34282 (L. Wölfer et al. 2025) and J1615 (M. Galloway-Sprietsma et al. 2025) show an intricate velocity structure. Similar to the case of LkCa 15, their inclinations of ~58° and ~47°, respectively, make it difficult to interpret their line centroid residuals, due to the mixed contributions of the different disk velocity components to the line-of-sight velocity. A forthcoming paper will present a quantitative analysis of the 2D residual maps and their comparison to predictions from large-scale turbulence for the complete exoALMA sample.

Finally, several exoALMA targets show evidence of dynamical interactions with massive planets or companions (C. Pinte et al. 2025). The perturbations driven by massive planets can drive strong signatures (e.g., L. Wölfer et al. 2021; A. F. Izquierdo et al. 2025), complicating the identification of substructures from other physical mechanisms. Therefore, planet–disk interactions can make it more difficult to characterize signatures of turbulent motions (e.g., M. Barraza-Alfaro et al. 2024). Follow-up higher spatial resolution observations exploiting the spectral resolution available with ALMA, combined with simultaneous predictions of velocity, temperature, and line width of disk perturbations, will lead to the robust characterization of large-scale turbulence signatures in exoALMA disks.

## 6.2. Localized Features from MRI-driven Perturbations

Localized features (in radius and azimuth) are often interpreted as the presence of embedded massive planets, due to their strong influence around their location. An embedded massive planet can manifest in the observables as a major “kink” spanning a set of channel maps (e.g., C. Pinte et al. 2019; J. Bae et al. 2025; C. Pinte et al. 2025), a “Doppler flip” (dipole pattern of blue- and redshifted velocities) in the line centroid residuals from Keplerian rotation (e.g., S. Pérez et al. 2018; S. Casassus et al. 2022; R. Teague et al. 2022; J. Bae et al. 2025), and a local excess in the line width around its location (A. F. Izquierdo et al. 2025). Interestingly, our $^{12}$CO (3–2) predictions of kinematic signatures from MRI-driven motions show a localized feature (highlighted in Figure 5). In particular, the line centroid map residuals for the model with an inclination of $i = 60^{\circ}$ show a localized signature at roughly $\varphi = 45^{\circ}$ west of north, seen in the rightmost panel of the middle row of Figures 3 and 4. Moreover, the localized velocity perturbation creates apparent wiggles in a limited range of channel maps. While we do not thoroughly compare such features with planet-driven velocity deviations, our predictions indicate that MRI-driven turbulence can create localized features that could be misinterpreted as embedded planets. However, we observe that MRI-driven perturbations substantially change in location, morphology, and strength with disk height. Therefore, the MRI origin of localized features can be ruled out if the detected signatures in channel maps and line centroid residuals are consistent among different tracers, as they probe distinct disk layers. In addition, our synthetic predictions do not show an excess of the $^{12}$CO (3–2) line width colocated with the MRI-induced localized feature (see Figure C6 in Appendix C), different from the predictions of signatures induced by embedded massive planets. Yet the line width map residuals of the mid-inclination disks are affected by artifacts, limiting the observability of the velocity-induced broadening from the localized feature. Our results highlight the importance of observing and analyzing localized perturbations in multiple molecular tracers, as performed by the exoALMA large program. Furthermore, our predictions stress the crucial role of alternative observables to the line centroid to robustly identify planet-driven signatures, such as the line peak and line width extracted with line profile tomography (A. F. Izquierdo et al. 2025). Still, a systematic study directly comparing the
observability of localized kinematic signatures from MRI and embedded planets in molecular line observations is needed to confirm the hinted features that could be used to distinguish between scenarios.

## 6.3. Limitations of Our Numerical Methods

Global 3D HD (MHD) simulations that include additional physics could provide new insights into the expected kinematic structure of magnetized protoplanetary disks. In particular, our post-processed nonideal MHD simulations ignore the effect of the Hall effect and ambipolar diffusion. Ambipolar diffusion can suppress the MRI activity in the outer disk (see G. Lesur et al. 2023, and references therein), while both ambipolar diffusion and the Hall effect can induce additional large-scale substructures (e.g., W. Béthune et al. 2016; C. Cui & X.-N. Bai 2021). These nonideal MHD effects can also modify the extent of the MRI dead zone and influence the gas dynamics interior to the MRI dead zone and near its edges (M. W. Kunz & G. Lesur 2013; G. Lesur et al. 2014; T. N. Delage et al. 2022; K. Iwasaki et al. 2024). Moreover, the strength of the ambipolar diffusion can define whether the MRI or the VSI dominates in outer regions of fast-cooling GI-stable disks (C. Cui & X.-N. Bai 2022). In addition, we neglect the magnetic effects on the GI gas dynamics in our SGHD simulation. In magnetized disks, the GI turbulent motions can drive a dynamo effect amplifying the magnetic field, resulting in a saturated state with less defined spiral arms (H. Deng et al. 2020; W. Béthune & H. Latter 2022), possibly more difficult to identify and characterize in CO kinematic observations.

The simulations presented in this paper are fine-tuned for the detailed study of disk instabilities; therefore, they are not tailored toward specific protoplanetary disks. Additionally, our models are set to explore a limited radial region of the disk (from 20 to 100 au), while the observed disks can extend up to several hundreds of au in size (e.g., M. Galloway-Sprietsma et al. 2025; C. Longarini et al. 2025). Further studies on identifying kinematic signatures from turbulence in the exoALMA sample will require source-specific numerical simulations and modeling, including information on each disk's structure. For example, we could use the disk temperature structure, stellar dynamical mass, disk size, disk mass, and other disk properties constrained from the exoALMA data (see, e.g., M. Galloway-Sprietsma et al. 2025; C. Longarini et al. 2025). In particular, the local pressure scale height, set by the disk temperature and increasing with disk radius, will define the radial size of substructures (J. Bae et al. 2023); therefore, it will set their observability within the spatial resolution of exoALMA. These new constraints on disk properties from the exoALMA observations add valuable information to guide a parameter space exploration in 3D HD simulations, which will substantially help us conclude whether the disk has the conditions to develop and sustain the steady state of the disk instabilities. By confirming that disk instabilities can drive the observed substructures unraveled via line profile tomography, we will ultimately have a robust detection of a turbulence mechanism active in planet-forming disks (e.g., J. Speedie et al. 2024).

## 6.4. Effect of Systematic Errors in the Characterization of Turbulence-driven Substructures

We have presented predictions of coherent kinematic substructures driven by large-scale motions from disk instabilities.

The observational signatures are extracted by subtracting a smooth background model for the disk rotation in the observed line centroid maps. However, such an approach is sensitive to various systematic errors, and careful extraction of the residuals is a must. First, small errors in the fitting parameters of a (sub- or super-) Keplerian model to the disk line centroid map obtained from rotational line observations can lead to global modulations (H.-W. Yen & P.-G. Gu 2020; B. J. Norfolk et al. 2022). In particular, patterns with characteristic dipole or quadrupole patterns will appear in the line centroid residual maps, prominent at the innermost and outermost regions of the disk. Such an artifact in the residual maps can introduce limitations on extracting coherent substructures, or errors in any analysis of the residual map properties. Ideally, methods to determine disk regions affected by systematic errors have to be applied during analysis to avoid considering them when interpreting the observations. In this direction, the data analysis tools developed within the exoALMA collaboration provide a significant step in advancing toward robust frameworks to extract observables and their residuals accurately (see, e.g., A. F. Izquierdo et al. 2025; T. E. Hilder et al. 2025). In addition, a comprehensive exploration of the kinematic analysis tools has been performed as part of exoALMA for post-processed 3D numerical simulations of planet-disk interactions, whose benchmarking will be fundamental for future studies of kinematic signatures (J. Bae et al. 2025).

## 6.5. Influence of Embedded Planets

The simulations presented in this work do not include the presence of embedded planets or companions in the disk. The presence of massive planets strongly modifies the disk's dynamical structure, inducing large-scale spiral arms, rings of sub- and super-Keplerian velocities, meridional flows, and vortices (for a review see S. Paardekooper et al. 2023, and references therein). Additional kinematic substructures driven by planets can alter the moment map residuals (A. F. Izquierdo et al. 2025), affecting the extraction of information from active turbulence. In particular, it has been demonstrated that massive planets can suppress observational signatures of the GI (S. Rowther et al. 2020) and the VSI (M. H. R. Stoll et al. 2017; M. Hammer & M.-K. Lin 2023; A. Ziampras et al. 2023; M. Barraza-Alfaro et al. 2024). Similarly, ideal MHD simulations have shown that a massive planet embedded in the disk can reduce the magnetic stresses, possibly weakening the MRI-driven turbulence along the gap region (R. P. Nelson & J. C. B. Papaloizou 2003; A. L. Uribe et al. 2011; M. Cilibrasi et al. 2023). Moreover, the inclusion of nonideal MHD effects and magnetic winds in magnetized planet-disk interaction simulations results in a complex disk dynamical structure (Y. Aoyama & X.-N. Bai 2023; G. Wafflard-Fernandez & G. Lesur 2023; X. Hu et al 2025), potentially hiding kinematic signatures of the MRI. Still, to properly address the influence of planets in the disk structure, a precise knowledge of the disk's physical properties and the location and masses of embedded planets is needed (see, e.g., C. Pinte et al. 2025). To this end, the constraints obtained from exoALMA will help to conduct dedicated simulations, exploring the kinematic signatures resulting from the interplay of turbulence-driving mechanisms and embedded planets.

## 7. Summary and Conclusions

We conducted a study of the large-scale turbulence signatures from numerical simulations of unstable protoplanetary disks. In
particular, we explored how coherent substructures are mani-
fested in the residuals from a background disk model in
exoALMA quality data of ¹²CO (3–2) line emission. Post-
processing 3D global HD (MHD) simulations with radiative
transfer calculations, we examined three disk instability
candidates to operate in the outer regions of planet-forming
disks: the VSI, MRI, and GI. From building a gallery of
expected outcomes, we put our results in the context of the
exoALMA survey, discussing the qualitative resemblance of the
perturbations extracted in exoALMA disks and our predictions
of substructures driven by disk instabilities. We summarize our
results in the following:

1. Our results show that the large-scale perturbations driven by the VSI, MRI, and GI are detectable by exoALMA quality ¹²CO (3–2) observations. Among the instabilities, the VSI shows the weakest perturbations with dominant ring- and arc-like substructures in the velocity centroid map residuals from a (sub-)Keplerian disk, decreasing in their degree of axisymmetry with increasing disk inclination.

2. Velocity perturbations driven by the MRI manifest as spiral-like substructures. Additionally, for mid-inclination disks, non-Keplerian rings at the MRI dead-zone outer edge are apparent. Similarly, the GI drives large-scale coherent substructures with spiral morphology, present in the residuals of line centroid, line peak, and line width.

3. We predict that the MRI can induce localized kinematic signatures in mid-inclination disks. These signatures may resemble the signposts induced by embedded massive planets. However, the MRI localized signature is observed only in a layer of the disk. Therefore, observations of multiple molecular tracers can help distinguish between scenarios.

4. By qualitatively comparing the overall morphology of perturbations from our predictions and the results with that from exoALMA observations, a robust case of VSI acting in the outermost regions of the disk is not yet observed.

5. The nondetection of spurs in the channel maps of PDS 66 and V4046 Sgr implies that their outer regions are nearly laminar, therefore suggesting that they are weakly ionized (stable to the MRI), have low gas surface densities (stable to the GI), and have relatively long cooling timescales (stable to the VSI). At the other end, dynamic exoALMA disks, such as MWC 758, HD 135344B, and CQ Tau, show apparent substructures with spiral morphologies, indicating that spiral-inducing mechanisms dominate (e.g., the MRI, the GI, or massive planets).

6. Spirals driven by the MRI or the GI could explain some of the substructures observed in exoALMA targets. However, these disks show evidence of dynamical interactions with embedded massive planets, resulting in a complex kinematic structure. Further work addressing the presence of disk instabilities from a methodical study of the statistical properties of residual maps is required to disentangle between scenarios.

Acknowledgments

This paper makes use of the following ALMA data: ADS/
JAO.ALMA#2021.1.01123.L. ALMA is a partnership of ESO
(representing its member states), NSF (USA) and NINS

(Japan), together with NRC (Canada), MOST and ASIAA
(Taiwan), and KASI (Republic of Korea), in cooperation with
the Republic of Chile. The Joint ALMA Observatory is
operated by ESO, AUI/NRAO and NAOJ. The National Radio
Astronomy Observatory is a facility of the National Science
Foundation operated under cooperative agreement by Asso-
ciated 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.

The authors thank the anonymous referee, whose feedback
significantly improved the manuscript. The authors acknowledge
the MIT Office of Research Computing and Data, the MIT
SuperCloud, and Lincoln Laboratory Supercomputing Center for
providing HPC resources that have contributed to the research
results reported within this paper. M.F. has received funding from
the European Research Council (ERC) under the European
Union's Horizon 2020 research and innovation program (grant
agreement No. 757957). J.B. acknowledges support from NASA
XRP grant No. 80NSSC23K1312. M.B., D.F., J.S., and A.J.W.
have received funding from the European Research Council
(ERC) under the European Union's Horizon 2020 research and
innovation program (PROTOPLANETS, grant agreement No.
101002188). Computations by D.F. and J.S. have been performed
on the "Mesocentre SIGAMM" machine, hosted by Observatoire
de la Cote d'Azur. P.C. acknowledges support by the Italian
Ministero dell'Istruzione, Università e Ricerca through the grant
Progetti Premiali 2012 - iALMA (CUP C52I13000140001) and
by the ANID BASAL project FB210003. S.F. is funded by the
European Union (ERC, UNVEIL, 101076613) and acknowledges
the financial contribution from PRIN-MUR 2022YP5ACE.
M.F. is supported by a grant-in-aid from the Japan Society for
the Promotion of Science (KAKENHI: No. JP22H01274).
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. G.L. and G.W.-F. acknowl-
edge support from the European Research Council (ERC) under
the European Union's Horizon 2020 research and innovation
program (grant agreement No. 815559 (MHDiscs)). 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 Austra-
lian Research Council funding via FT170100040, DP18010423,
DP220103767, and DP240103290. D.P. acknowledges Australian
Research Council funding via DP18010423, DP220103767, and
DP240103290. G.R. acknowledges funding from the Fondazione
Cariplo, grant No. 2022-1217, and the European Research
Council (ERC) under the European Union's Horizon Europe
Research & Innovation Programme under grant agreement No.
101039651 (DiscEvol). 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. H.-W.Y.
acknowledges support from the 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. Support for B.Z. was provided by the Brinson Foundation. Views and opinions expressed by ERC-funded scientists are, however, those of the author(s) only and do not necessarily reflect those of the European Union or the European Research Council. Neither the European Union nor the granting authority can be held responsible for them.

**Software:** PLUTO (A. Mignone et al. 2007), RADMC-3D (C. P. Dullemond et al. 2012), fargo2radmc3d (C. Baruteau et al. 2019), CASA (J. P. McMullin et al. 2007), betterments (R. Teague & D. Foreman-Mackey 2018), emcee (D. Foreman-Mackey et al. 2013), eddy (R. Teague 2019a), GoFish (R. Teague 2019b), Numpy (C. R. Harris et al. 2020), Scipy (P. Virtanen et al. 2020), VTK (W. Schroeder et al. 2006), Matplotlib (J. D. Hunter 2007), CMasher (E. van der Velden 2020).

# Appendix A
## Eddy Disk Model and Fitting Procedure

In our EDDY disk models described in Section 3.3, we include a modified Keplerian rotation profile that includes a radially varying deviation from purely Keplerian rotation, parameterized by a mass component (R. Teague et al. 2022):

$$v_{\phi} = \sqrt{\frac{G(M_{*} + M_{d}(R))R^{2}}{(R^{2} + z^{2})^{3/2}}}, \quad (A1)$$

$$M_{d}(R) = M_{\text{disk}} \times \frac{R^{2-\gamma} - R_{\text{in}}^{2-\gamma}}{R_{\text{out}}^{2-\gamma} - R_{\text{in}}^{2-\gamma}}, \quad (A2)$$

where R is the cylindrical radius, z is the disk surface height, and $M_*$ is the dynamical mass of the central star. The parameters $R_{\text{in}}$ and $R_{\text{out}}$ describe the disk inner and outer edge, respectively. $M_{\text{disk}}$ represents the total disk mass, and $\gamma$ is the exponent of the underlying power-law surface density profile, assumed to be of the form $\Sigma(r) = \Sigma_0 \times R^{-\gamma}$.

The parameterization of the emission surface in our models is given by

$$z(R^*) = z_0 \times \left(\frac{R^*}{1''}\right)^{\psi} \times \exp\left(-\left(\frac{R^*}{R_{\text{taper}}}\right)^{q_{\text{taper}}}\right), \quad (A3)$$

where the cylindrical radius $R^*$ and scale radius of the exponential taper $R_{\text{taper}}$ are implemented as angular distances defined in arcseconds. The disk aspect ratio of the emission surface is set by $z_0$, $\psi$ describes the surface flaring, and $R_{\text{taper}}$ and $q_{\text{taper}}$ describe an exponential taper to the emission height at the disk outermost region. Finally, the fitting is done in the projected velocity into the line of sight $v_0$, described as

$$v_0 = v_{\phi} \cos \phi \sin i + v_{\text{LSR}}, \quad (A4)$$

where $\phi$ is the polar angle, $i$ is the disk inclination, and $v_{\text{LSR}}$ is the systemic velocity. The free parameters in our calculations are $M_*$, the emission surface parameters [$z_0$; $\psi$; $R_{\text{taper}}$; $q_{\text{taper}}$], the disk inclination $i$, the disk mass parameters [$M_{\text{disk}}$; $\gamma$], and the center of the image [$x_0$, $y_0$]. We fixed the disk PA to the input model value, set the systemic velocity $v_{\text{LSR}}$ to zero in our models, fixed the distance to the source to 140 pc, and set the

inner and outer disk edges considered in the disk mass component $r_{\text{in}}$ and $r_{\text{out}}$ to 20 and 100 au, respectively.

To obtain our best-fit disk models, we computed the posterior distribution of our free parameters with Markov Chain Monte Carlo sampling implemented with the `emcee` package (D. Foreman-Mackey et al. 2013) in `eddy`. We set 80 walkers initialized with the disk radiative transfer input model parameters for $M_*$ and $M_{\text{disk}}$, while [$\gamma = 1$; $x_0 = 0$, $y_0 = 0$] were used as initial values. For the emission surface parameters, an approximated location was obtained by roughly reproducing the layer where the optical depth reaches unity in our model (see Figure C4 in Appendix C). We set 2000 steps for the walkers to burn in and 4000 additional steps to sample the posterior distribution.

# Appendix B
## Parameters of the Simulations

In Table 1 of this Appendix we provide a summary of the main parameters of the numerical simulations used in our study and described in Section 2.

<table><caption>**Table 1**
Summary of the Simulations' Parameters</caption><thead><tr><th>Label</th><th>RHD [VSI]</th><th>MHD [MRI]</th><th>SGHD [GI]</th></tr></thead><tbody><tr><td>L<sub>*</sub></td><td>0.95 L<sub>☉</sub></td><td>0.95 L<sub>☉</sub></td><td>...</td></tr><tr><td>R<sub>*</sub></td><td>2.0 R<sub>☉</sub></td><td>2.0 R<sub>☉</sub></td><td>...</td></tr><tr><td>M<sub>*</sub></td><td>0.5 M<sub>☉</sub></td><td>0.5 M<sub>☉</sub></td><td>0.5 M<sub>☉</sub></td></tr><tr><td>T<sub>*</sub></td><td>4000 K</td><td>4000 K</td><td>...</td></tr><tr><td>Disk total gas mass</td><td>0.05 M<sub>*</sub></td><td>0.085 M<sub>*</sub></td><td>0.2 M<sub>*</sub></td></tr><tr><td>Σ<sub>0</sub> at R = 100 au</td><td>6.0 g cm<sup>-2</sup></td><td>5.94 g cm<sup>-2</sup></td><td>...</td></tr><tr><td>Spectral type</td><td>K</td><td>K</td><td>...</td></tr><tr><td>Dust-to-gas mass ratio</td><td>10<sup>-3</sup></td><td>10<sup>-2</sup></td><td>...</td></tr><tr><td>Maximum dust size</td><td>10 µm</td><td>10 µm</td><td>...</td></tr><tr><td>Minimum dust size</td><td>0.01 µm</td><td>0.01 µm</td><td>...</td></tr><tr><td>Dust size slope</td><td>-3.5</td><td>-3.5</td><td>...</td></tr><tr><td>Δr</td><td>20–100 au</td><td>20–100 au</td><td>5–160 au</td></tr><tr><td>Δθ</td><td>0.7 rad</td><td>0.72 rad</td><td>0.7 rad</td></tr><tr><td>Δφ</td><td>2π rad</td><td>2π rad</td><td>2π rad</td></tr><tr><td>N<sub>r</sub></td><td>1024</td><td>256</td><td>518</td></tr><tr><td>N<sub>θ</sub></td><td>512</td><td>128</td><td>96</td></tr><tr><td>N<sub>φ</sub></td><td>2044</td><td>512</td><td>512</td></tr><tr><td>Reference H/R</td><td>0.1</td><td>0.1</td><td>0.028</td></tr><tr><td>Cells per H</td><td>70</td><td>20</td><td>9*</td></tr><tr><td>⟨α⟩</td><td>1.5 × 10<sup>-4</sup></td><td>~ 3 × 10<sup>-3</sup></td><td>9.7 × 10<sup>-2</sup></td></tr></tbody></table>

**Note.** See further details in M. Flock et al. (2015, 2020) and W. Béthune et al. (2021) for the MHD, RHD, and SGHD simulations, respectively. Disk- and time-averaged values of the effective $\alpha$ viscosity (N. I. Shakura & R. A. Sunyev 1973) are reported. An asterisk indicates cells per pressure scale height near the disk midplane.

# Appendix C
## Additional Figures

In this Appendix we present additional views of our simulations' snapshots, the radiative transfer model structure, and alternative observables to the line centroid. In Figures C1 and C2, we provide a Cartesian view of the simulations at the disk midplane, and a Polar view of the simulations at $Z \sim 2.5H$, respectively. In Figure C3, we show the full domain of the SGHD simulation, post-processed to study the signatures driven by GI. In Figure C4, we show the thermal structure and
<figure><img src="image_9.png" /><figcaption>**Figure C1.** Cartesian view of perturbations at the disk midplane for the HD (top), MHD (middle), and SGHD (bottom) simulations, with active VSI, MRI, and GI, respectively. From left to right, we show gas density perturbations, radial velocity, meridional velocity, and azimuthal velocity perturbations. The color bar limits are adjusted for each instability.</figcaption></figure>

number density structure of the radiative transfer models,
indicating the height at which the optical depth reaches unity.
Finally, in Figures C5 and C6 we show the observable

residual maps of the line peak intensity and line width,
relative to an azimuthal disk average (described in
Section 5.2).
<figure><img src="image_10.png" /><figcaption>**Figure C2.** Polar view of the disk perturbations at $Z \sim 2.5H$ for the HD (top), MHD (middle), and SGHD (bottom) simulations, with active VSI, MRI, and GI, respectively. From left to right, we show perturbations in density, radial velocity, meridional velocity, and azimuthal velocity. The color bar limits are adjusted for each instability.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure C3.** Cartesian view of perturbations at $Z \sim 2.5H$ for the SGHD simulation with active GI, showing the simulation full domain. From left to right, we show gas density perturbations, radial velocity, meridional velocity, and azimuthal velocity perturbations. The dotted lines enclose the radial region of the simulation considered in our study.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>**Figure C4.** Gas temperature and azimuthally averaged number density of $^{12}$CO of our disk radiative transfer models. Overlaid in blue, a solid line traces the height at which the optical depth reaches unity assuming a face-on disk. The red dashed-dotted line shows the initial guess for the $^{12}$CO (3–2) emission surface in our EDDY models. The solid black line indicates the simulations’ domain in colatitude, while the gray dashed lines show spacings along Z equal to the disks’ pressure scale height.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure C5.** Temperature residuals from $^{12}$CO (3–2) predictions of inclined disks unstable to VSI, MRI, and GI. The residuals are extracted by computing the disks’ eighth moment maps and then subtracting their azimuthal averages using GOFISH. The models only consider the upper (or front) CO layer. The images have a 104 mas$^{-1}$ velocity resolution and a 35 mas angular resolution and include the noise level achieved by exoALMA. From left to right, different disk inclinations are shown. The beam size is shown in the lower left corner of each panel.</figcaption></figure>
<figure><img src="image_14.png" /></figure>

**Figure C6.** Line width residuals from ¹²CO (3–2) predictions of inclined disks unstable to VSI, MRI, and GI. The residuals are extracted by computing the disks’ line width maps using the percentiles method in EDDY and then subtracting their azimuthal averages using GOFISH. The models only consider the upper (or front) CO layer. The images have a 104 m s⁻¹ velocity resolution and a 35 mas angular resolution and include the noise level achieved by exoALMA. From left to right, different disk inclinations are shown. The beam size is shown in the lower left corner of each panel. The color bar limits have been adjusted for each instability.

## ORCID iDs

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

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

William Béthune https://orcid.org/0009-0001-2471-0336

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

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

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

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

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

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

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

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

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

Eric W. Koch https://orcid.org/0000-0001-9605-780X

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

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

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

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

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

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

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

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

## References

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

Aoyama, Y., & Bai, X.-N. 2023, ApJ, 946, 5

Bae, J., Flock, M., Izquierdo, A., et al. 2025, ApJL, 984, L12

Bae, J., Isella, A., Zhu, Z., et al. 2023, in ASP Conf. Ser. 534, Protostars and Planets VII, ed. S. Inutsuka et al. (San Francisco, CA: ASP), 423

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

Bai, X.-N. 2015, ApJ, 798, 84

Bai, X.-N., & Stone, J. M. 2017, ApJ, 836, 46

Balbus, S. A., & Hawley, J. F. 1991, ApJ, 376, 214

Balbus, S. A., & Terquem, C. 2001, ApJ, 552, 235

Barraza-Alfaro, M., Flock, M., & Henning, T. 2024, A&A, 683, A16

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

Baruteau, C., Barraza, M., Pérez, S., et al. 2019, MNRAS, 486, 304

Baruteau, C., Wafflard-Fernandez, G., Le Gal, R., et al. 2021, MNRAS, 505, 359

Béthune, W., & Latter, H. 2022, A&A, 663, A138

Béthune, W., Latter, H., & Kley, W. 2021, A&A, 650, A49

Béthune, W., Lesur, G., & Ferreira, J. 2016, A&A, 589, A87

Béthune, W., Lesur, G., & Ferreira, J. 2017, A&A, 600, A75

Birnstiel, T. 2024, ARA&A, 62, 157

Brandenburg, A., & Nordlund, Å. 2011, RPPh, 74, 046901
Casassus, S., Cárcamo, M., Hales, A., Weber, P., & Dent, B. 2022, ApJ, 933, L4

Casassus, S., Christiaens, V., Cárcamo, M., et al. 2021, MNRAS, 507, 3789

Cilibrasi, M., Flock, M., & Szulágyi, J. 2023, MNRAS, 523, 2039

Cossins, P., Lodato, G., & Clarke, C. J. 2009, MNRAS, 393, 1157

Cui, C., & Bai, X.-N. 2020, ApJ, 891, 30

Cui, C., & Bai, X.-N. 2021, MNRAS, 507, 1106

Cui, C., & Bai, X.-N. 2022, MNRAS, 516, 4660

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

Delage, T. N., Okuzumi, S., Flock, M., Pinilla, P., & Dzyurkevich, N. 2022, A&A, 658, A97

Deng, H., Mayer, L., & Latter, H. 2020, ApJ, 891, 154

Disk Dynamics Collaboration, Armitage, P. J., & Bae, J. 2020, arXiv:2009.04345

Dong, R., Hall, C., Rice, K., & Chiang, E. 2015, ApJ, 812, L32

Dražkowska, J., Bitsch, B., Lambrechts, M., et al. 2023, in ASP Conf. Ser. 534, Protostars and Planets VII, ed. S. Inutsuka et al. (San Francisco, CA: ASP), 717

Dullemond, C. P., Juhasz, A., Pohl, A., et al., 2012 RADMC-3D: A multi-purpose radiative transfer tool, Astrophysics Source Code Library, ascl:1202.015

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

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

Flaherty, K. M., Hughes, A. M., Rose, S. C., et al. 2017, ApJ, 843, 150

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

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

Flock, M., Dzyurkevich, N., Klahr, H., Turner, N. J., & Henning, T. 2011, ApJ, 735, 122

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

Flock, M., Nelson, R. P., Turner, N. J., et al. 2017, ApJ, 850, 131

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

Flock, M., Turner, N. J., Nelson, R. P., et al. 2020, ApJ, 897, 155

Flores-Rivera, L., Flock, M., & Nakatani, R. 2020, A&A, 644, A50

Foreman-Mackey, D., Conley, A., Meierjurgen Farr, W., et al., 2013 emcee: The MCM. Hammer, Astrophysics Source Code Library, ascl:1303.002

Forgan, D., Armitage, P. J., & Simon, J. B. 2012, MNRAS, 426, 2419

Forgan, D., Rice, K., Cossins, P., & Lodato, G. 2011, MNRAS, 410, 994

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

Guilloteau, S., Dutrey, A., Wakelam, V., et al. 2012, A&A, 548, A70

Hacar, A., Alves, J., Burkert, A., & Goldsmith, P. 2016, A&A, 591, A104

Hall, C., Dong, R., Teague, R., et al. 2020, ApJ, 904, 148

Hall, C., Rice, K., Dipierro, G., et al. 2018, MNRAS, 477, 1004

Hammer, M., & Lin, M.-K. 2023, MNRAS, 525, 123

Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Natur, 585, 357

Hawley, J. F., & Balbus, S. A. 1991, ApJ, 376, 223

Hilder, T. E., Casey, A. R., Price, D. J., et al. 2025, ApJL, 984, L13

Hu, X., Li, Z.-Y., Bae, J., & Zhu, Z. 2025, MNRAS, 536, 1374

Hughes, A. M., Wilner, D. J., Andrews, S. M., Qi, C., & Hogerheijde, M. R. 2011, ApJ, 727, 85

Hunter, J. D. 2007, CSE, 9, 90

Iwasaki, K., Tomida, K., Takasao, S., Okuzumi, S., & Suzuki, T. K. 2024, PASJ, 76, 616

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

Izquierdo, A. F., Facchini, S., Rosotti, G. P., van Dishoeck, E. F., & Testi, L. 2022, ApJ, 928, 2

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

Klahr, H., & Hubbard, A. 2014, ApJ, 788, 21

Klahr, H., Pfeil, T., & Schreiber, A. 2018, in Handbook of Exoplanets, ed. H. J. Deeg & J. A. Belmonte, 138 (Cham: Springer International Publishing), 2251

Koch, E. W., & Rosolowsky, E. W., 2016 FilFinder: Filamentary structure in molecular clouds, Astrophysics Source Code Library, ascl:1608.009

Krapp, L., Gressel, O., Benítez-Llambay, P., et al. 2018, ApJ, 865, 105

Kratter, K., & Lodato, G. 2016, ARA&A, 54, 271

Kunz, M. W., & Lesur, G. 2013, MNRAS, 434, 2295

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

Lesur, G., Flock, M., Ecolano, B., et al. 2023, in ASP Conf. Ser. 534, Astronomical Society of the Pacific Conference Series, ed. S. Inutsuka et al. (San Francisco, CA: ASP), 465

Lesur, G., Kunz, M. W., & Fromang, S. 2014, A&A, 566, A56

Lin, M.-K., & Youdin, A. N. 2015, ApJ, 811, 17

Lodato, G., & Rice, W. K. M. 2004, MNRAS, 351, 630

Longarini, C., Lodato, G., Toci, C., et al. 2021, ApJL, 920, L41

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

Lyra, W. 2014, ApJ, 789, 77

Lyra, W., & Umurhan, O. M. 2019, PASP, 131, 072001

Malygin, M. G., Klahr, H., Semenov, D., Henning, T., & Dullemond, C. P. 2017, A&A, 605, A30

Marcus, P. S., Pei, S., Jiang, C.-H., & Barranco, J. A. 2016, ApJ, 833, 148

Marcus, P. S., Pei, S., Jiang, C.-H., et al. 2015, ApJ, 808, 87

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

Mignone, A., Bodo, G., Massaglia, S., et al. 2007, ApJS, 170, 228

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

Nelson, R. P., & Papaloizou, J. C. B. 2003, MNRAS, 339, 993

Norfolk, B. J., Pinte, C., Calcino, J., et al. 2022, ApJL, 936, L4

Paardekooper, S., Dong, R., Duffell, P., et al. 2023, in ASP Conf. Ser. 534, Astronomical Society of the Pacific Conference Series, ed. S. Inutsuka et al. (San Francisco, CA: ASP), 685

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

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

Paneco-Carreño, T., Pérez, L. M., Benisty, M., et al. 2021, ApJ, 914, 88

Pérez, S., Casassus, S., & Benitez-Llambay, P. 2018, MNRAS, 480, L12

Pfeil, T., & Klahr, H. 2019, ApJ, 871, 150

Pfeil, T., & Klahr, H. 2021, ApJ, 915, 130

Pinte, C., Price, D. J., Ménard, F., et al. 2018, ApJL, 860, L13

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

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

Pinte, C., van der Plas, G., Ménard, F., et al. 2019, NatAs, 3, 1109

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

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

Ross, J., & Latter, H. N. 2018, MNRAS, 477, 3329

Rotherw, S., Meru, F., Kennedy, G. M., Nealon, R., & Pinte, C. 2020, ApJL, 904, L18

Schroeder, W., Martin, K., & Lorensen, B. 2006, The Visualization Toolkit (4th ed.; New York: Kitware)

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

Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337

Shi, J.-M., & Chiang, E. 2014, ApJ, 789, 34

Simon, J. B., & Armitage, P. J. 2014, ApJ, 784, 15

Simon, J. B., Bai, X.-N., Flaherty, K. M., & Hughes, A. M. 2018, ApJ, 865, 10

Simon, J. B., Hughes, A. M., Flaherty, K. M., Bai, X.-N., & Armitage, P. J. 2015, ApJ, 808, 180

Speedie, J., Dong, R., Hall, C., et al. 2024, Natur, 633, 58

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

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

Stoll, M. H. R., Picogna, G., & Kley, W. 2017, A&A, 604, A28

Teague, R. 2019a, JOSS, 4, 1220

Teague, R. 2019b, JOSS, 4, 1632

Teague, R., Bae, J., Andrews, S. M., et al. 2022, ApJ, 936, 163

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

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

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

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

Terry, J. P., Hall, C., Longarini, C., et al. 2022, MNRAS, 510, 1671

Toomre, A. 1964, ApJ, 139, 1217

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

Turner, N. J., Fromang, S., Gammie, C., et al. 2014, in Protostars and Planets VI, ed. H. Beuther et al. (Tucson, AZ: Univ. of Arizona Press), 411

Uribe, A. L., Klahr, H., Flock, M., & Henning, T. 2011, ApJ, 736, 85

Urpin, V., & Brandenburg, A. 1998, MNRAS, 294, 399

van der Velden, E. 2020, JOSS, 5, 2004

Veronesi, B., Longarini, C., Lodato, G., et al. 2024, A&A, 688, A136

Veronesi, B., Paneque-Carreño, T., Lodato, G., et al. 2021, ApJL, 914, L27

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

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

Wafflard-Fernandez, G., & Lesur, G. 2023, A&A, 677, A70

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

Wölfer, L., Barraza-Alfaro, M., Teague, R., et al. 2025, ApJL, 984, L22

Yen, H.-W., & Gu, P.-G. 2020, ApJ, 905, 89

Yun, H.-G., Kim, W.-T., Bae, J., & Han, C. 2025a, ApJ, 980, 14

Yun, H. G., Kim, W.-T., Bae, J., & Han, C. 2025b, ApJ, 980, 15

Ziampras, A., Kley, W., & Nelson, R. P. 2023, A&A, 670, A135