# Gravitational instability in a planet-forming disk

Jessica Speedie, Ruobing Dong, Cassandra Hall, Cristiano Longarini, Benedetta Veronesi, Teresa Paneque-Carreño, Giuseppe Lodato, Ya-Wen Tang, Richard Teague, Jun Hashimoto

► To cite this version:

Jessica Speedie, Ruobing Dong, Cassandra Hall, Cristiano Longarini, Benedetta Veronesi, et al. Gravitational instability in a planet-forming disk. Nature, 2024, 633, pp.58-62. 10.1038/s41586-024-07877-0. insu-04838849

HAL Id: insu-04838849

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

Submitted on 8 Oct 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
# Gravitational instability in a planet-forming disk

https://doi.org/10.1038/s41586-024-07877-0

Received: 11 December 2023

Accepted: 25 July 2024

Published online: 4 September 2024

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

Check for updates

Jessica Speedie<fnref n="1" />,², Ruobing Dong<fnref n="1" />,²,³, Cassandra Hall³,⁴, Cristiano Longarini⁵,⁶,
Benedetta Veronesi⁷, Teresa Paneque-Carreño⁸,⁹, Giuseppe Lodato⁵, Ya-Wen Tang¹⁰,
Richard Teague¹¹ & Jun Hashimoto¹²,¹³,¹⁴

The canonical theory for planet formation in circumstellar disks proposes that planets are grown from initially much smaller seeds$^{1-5}$. The long-considered alternative theory proposes that giant protoplanets can be formed directly from collapsing fragments of vast spiral arms$^{6-11}$ induced by gravitational instability$^{12-14}$—if the disk is gravitationally unstable. For this to be possible, the disk must be massive compared with the central star: a disk-to-star mass ratio of 1:10 is widely held as the rough threshold for triggering gravitational instability, inciting substantial non-Keplerian dynamics and generating prominent spiral arms$^{15-18}$. Although estimating disk masses has historically been challenging$^{19-21}$, the motion of the gas can reveal the presence of gravitational instability through its effect on the disk-velocity structure$^{22-24}$. Here we present kinematic evidence of gravitational instability in the disk around AB Aurigae, using deep observations of $^{13}$CO and $^{18}$O line emission with the Atacama Large Millimeter/submillimeter Array (ALMA). The observed kinematic signals strongly resemble predictions from simulations and analytic modelling. From quantitative comparisons, we infer a disk mass of up to a third of the stellar mass enclosed within 1" to 5" on the sky.

We targeted the disk around AB Aurigae (hereafter AB Aur), a 2.5–4.4-Myr-old$^{25-28}$ Herbig Ae (ref. 29) star of intermediate mass ($M_* = 2.4 M_\odot$)$^{26,27,30}$ at a distance of 155.9 ± 0.9 pc (ref. 31). AB Aur is at a relatively late stage of protostellar evolution, classified as a Class II young stellar object$^{32,33}$. To investigate the velocity structure of the disk, we obtained deep ALMA Band 6 observations of molecular emission lines $^{13}$CO (J = 2–1) and $^{18}$O (J = 2–1) with high velocity resolution (channel widths of $v_{chan} = 42$ m s$^{-1}$ and 84 m s$^{-1}$, respectively). The observations were taken in two array configurations with baselines ranging from 14 to 2,216 m, reaching a total on-source integration time of 5.75 h. Imaging with a Briggs robust value of 0.5 provided image cubes with a spatial resolution or beam size of 0.237" × 0.175" (beam position angle (PA) = 1.2°), equivalent to 37 × 27 AU. We collapse the 3D image cubes into 2D moment maps to expose the velocity-integrated intensity (moment 0), intensity-weighted line-of-sight velocity ($v_{los}$, moment 1) and emission linewidth (moment 2). This collection is shown in Extended Data Fig. 1.

To reveal the spiral arms in the disk, we apply a high-pass filter$^{34}$ (see Methods) to the ALMA $^{13}$CO moment maps (Fig. 1b–d). In the filtered line-of-sight velocity (moment 1) map, we observe spiral-shaped disturbances in the gas velocity field throughout the disk (Fig. 1b). With the filtered velocity-integrated intensity (moment 0) and linewidth

(moment 2) maps, we visually highlight regions of peak density and temperature (Fig. 1c,d). Compression and shock heating are expected to lead to temperature enhancements (and thus localized line broadening) within gravitational instability (GI)-induced density spirals in self-regulating disks$^{13,23}$. The VLT/SPHERE H-band scattered-light image of AB Aur originally presented in ref. 35 is shown for comparison (Fig. 1a). Scattered light comes from the disk surface, tracing the distribution of (sub-)micrometre-sized dust usually well coupled with the gas. Previous simulations have shown that GI-induced density spirals are prominent in scattered light$^{16,36}$. At least seven spiral structures (S1–S7) have been previously identified in the H-band image$^{35,37}$, although not all occupy the same radial region and some may be branches of adjacent arms$^{38}$. The disk rotates anticlockwise (the spiral arms are trailing) and the south side is the near side, tilted towards us$^{38-40}$.

To provide a qualitative comparison with the ALMA observations, we run 3D smoothed-particle hydrodynamic (SPH) simulations of a gravitationally unstable disk (see Methods). The simulations were post-processed with radiative transfer and then further processed to have the same viewing angle, sensitivity, spectral and angular resolution as the AB Aur data. To place the disk comfortably within the gravitationally unstable regime ($M_{disk}/M_* \ge 0.1$), we set the total gas mass to

<footnote n="1">¹Department of Physics and Astronomy, University of Victoria, Victoria, British Columbia, Canada. ²Kavli Institute for Astronomy and Astrophysics, Peking University, Beijing, People's Republic of China. ³Department of Physics and Astronomy, The University of Georgia, Athens, GA, USA. ⁴Center for Simulational Physics, The University of Georgia, Athens, GA, USA. ⁵Università degli Studi di Milano, Milan, Italy. ⁶Institute of Astronomy, University of Cambridge, Cambridge, UK. ⁷Univ Lyon, Univ Lyon1, Ens de Lyon, CNRS, Centre de Recherche Astrophysique de Lyon UMR5574, Saint-Genis-Laval, France. ⁸Leiden Observatory, Leiden University, Leiden, The Netherlands. ⁹European Southern Observatory, Garching, Germany. ¹⁰Institute of Astronomy and Astrophysics, Academia Sinica, Taipei, Taiwan. ¹¹Department of Earth, Atmospheric, and Planetary Sciences, Massachusetts Institute of Technology, Cambridge, MA, USA. ¹²Astrobiology Center, National Institutes of Natural Sciences, Mitaka, Japan. ¹³Subaru Telescope, National Astronomical Observatory of Japan, Mitaka, Japan. ¹⁴Department of Astronomy, School of Science, Graduate University for Advanced Studies (SOKENDAI), Mitaka, Japan. ¹⁵e-mail: jspeedie@uvic.ca; rbdong@uvic.ca</footnote>
<figure><img src="image_3.png" /></figure>

**Fig. 1** | Global spirals in the AB Aur disk. **a,** VLT/SPHERE H-band scattered-light image of the AB Aur disk<sup>35</sup> tracing a spiral structure in (sub-)micrometre-sized dust grains. The labelled spirals S1–S7 are taken from previous works<sup>37,38</sup>. **b,** Filtered ALMA <sup>13</sup>CO intensity-weighted mean velocity (moment 1) map, revealing residual gas motion within the bulk flow. The synthesized beam is shown in the bottom-left corner as an ellipse. The inset zooms into the region around where S5 crosses the minor axis, highlighting converging flows on the

two sides of S5 indicated by arrows. **c,** Filtered ALMA <sup>13</sup>CO integrated intensity (moment 0) map, highlighting peaks in the gas density and/or temperature. **d,** Filtered ALMA <sup>13</sup>CO emission linewidth (moment 2) map, showing localized line broadening within the spiral arms. Insets in **c** and **d** zoom into the same region as the inset in **b**, showing enhanced gas density and/or temperature caused by the radially converging flows around S5.

0.3 times the mass of the star. For sustained spiral arms, we set the cooling timescale to ten times the local dynamical timescale (β = 10). The simulated GI disk shows spiral structures in all three moment maps, resembling those in the AB Aur disk (Extended Data Figs. 1 and 2). Overall, the AB Aur disk hosts a global architecture of spiral arms at 100-AU to 1,000-AU scales across all azimuths in multiwavelength observations tracing different disk components and quantities, strongly indicating ongoing GI.

One characteristic kinematic feature in the AB Aur disk can be found in the isovelocity curve at the systemic velocity $v_{\text{sys}}$ in the moment 1 map. Figure 2a shows a sinusoidal pattern at $v_{\text{los}} = v_{\text{sys}}$ (along the minor axis; white colour), more prominent towards the south. This signature, known as a 'minor-axis GI wiggle'<sup>22</sup>, has been predicted in hydrodynamic simulations<sup>22,24</sup> and analytic theory<sup>23</sup> as a clear kinematic signature of GI (Fig. 2b,c). It is one of a global set of GI wiggles in isovelocity curves that we observe throughout the AB Aur disk (Extended Data Fig. 3). These wiggles are generated by self-gravitating spiral arms, which constitute local minima in the gravitational potential field and induce corresponding oscillations in the gas velocity field. The synthetic moment 1 map of the SPH GI disk simulation shows a minor-axis GI wiggle with similar

morphology as the observed one (Fig. 2c), completely distinct from the linear pattern found in a disk undergoing Keplerian rotation with no radial motions (Fig. 2b,c insets).

Among all GI wiggles, the minor-axis GI wiggle has been known and targeted in past studies for its convenience in quantitative analysis<sup>23,24</sup>. Owing to projection effects, only the radial and vertical components of the disk-velocity field ($v_r$ or $v_z$) contribute to $v_{\text{los}}$ at the systemic velocity traced by this wiggle. In the case of GI-induced velocity perturbations, the $v_r$ contribution is expected to dominate<sup>22</sup>. As we show with 2D analytic calculations of gravitationally unstable disks (see Methods), a self-gravitating spiral arm induces radial motion convergent on itself, appearing as a wiggle in the moment 1 map at $v_{\text{sys}}$ where the spiral crosses the minor axis (see Extended Data Fig. 5). The filtered moment 1 map in Fig. 1b shows redshift and blueshift patterns corresponding to convergent flows towards spiral S5 (visible in both scattered-light and <sup>13</sup>CO moments 0 and 2; Fig. 1a,c,d), supporting the interpretation that the GI wiggle along the southern minor axis in Fig. 2a is generated by a self-gravitating spiral arm.

Having identified evidence of GI in disk kinematics and in the detections of spirals across several tracers and moment maps, we now
<figure><img src="image_4.png" /><figcaption>**Fig. 2** | Detection of the Gl wiggle in the AB Aur disk. **a**, ALMA $^{13}$CO intensity-weighted mean velocity (moment 1) map showing $v_{\text{los}}$ of gas in the AB Aur disk. The observations show the 'Gl wiggle' along the minor axis (arrow) predicted in ref. 22 as a clear kinematic signature of Gl. **b**, $v_{\text{los}}$ map of a gravitationally unstable disk at the inclination and PA of the AB Aur disk, computed with 2D analytic modelling$^{23}$. Self-gravitating spiral arms crossing the minor axis</figcaption></figure>

induce radial motion that appears as a wiggle (arrow). **c**, Synthetic ALMA $^{13}$CO moment 1 map of the 3D SPH Gl disk simulation, revealing the same Gl wiggle signature (arrow). The insets in **b** and **c** show corresponding images for Keplerian disks with no radial gas motion, in which the isovelocity curve at the systematic velocity appears as a straight line along the minor axis.

quantitatively analyse the Gl wiggle along the southern minor axis to constrain the disk mass. We extract the $^{13}$CO and $^{18}$O emission spectra along the southern disk minor axis (Fig. 3a,b) and detect the wiggle in position-velocity space (hereafter referred to as the ‘PV wiggle’), which is a different view of the position-position wiggle in Fig. 2a. Slicing the 3D image cubes this way more comprehensively exposes the gas velocity structure and enables us to quantify the perturbation in units of velocity. We measure the emission line centres by performing a quadratic fit to the spectrum in each spatial pixel of the image cube$^{41}$. This method achieves sub-spectral-resolution precision on the line centre and yields statistically meaningful and robust uncertainties$^{42}$. We find remarkably similar sinusoidal morphology between the PV wiggles in $^{13}$CO and $^{18}$O emission (Fig. 4a).

Theoretical studies have shown that the dynamical response of a
disk to its own self-gravity is sensitive to the disk-to-star mass ratio

and the cooling rate$^{23,24}$. Specifically, the amplitude of the induced
radial-velocity perturbations is proportional to ($M_{\text{disk}}/M_{\text{⋆}}$)$^2$ and $\beta^{-1/2}$
(equations (11) and (18) in the Methods). This allows us to use the
observed minor-axis PV wiggle to infer the disk mass once we make
assumptions on the disk cooling rates. Following ref. 23, we use a sta-
tistical metric to quantify the ‘magnitude’ of the minor-axis PV wiggle,
defined as the standard deviation of the line-centre velocities over a
radial range. Bounded by the inner central cavity and outer edge of
recovered $C^{18}$O emission, our radial range spans 1" to 5" (155 to 780 AU).
We find a magnitude of $37.4 \pm 2.9 \ \text{m s}^{-1}$ for the southern minor-axis PV
wiggle in $^{13}$CO and $44.2 \pm 1.3 \ \text{m s}^{-1}$ in $C^{18}$O (Fig. 4b). For comparison,
the gravitationally unstable disk in the SPH simulation has a south-
ern minor-axis PV wiggle in $^{13}$CO emission with quantitatively similar
altitude and sinusoidal morphology (Fig. 3c) and a magnitude of
$39.1 \pm 1.8 \ \text{m s}^{-1}$ (Extended Data Fig. 7a).

<figure><img src="image_5.png" /><figcaption>**Fig. 3** | The PV wiggle. Emission spectra (intensity as a function of velocity) extracted along the southern minor axis of the disk, plotted with distance from the star. The line centres are shown as yellow points. The insets show the corresponding line-centre map, with the circles delineating 1" radial increments. The yellow line along the southern minor axis is the narrow (0.5"-wide) wedge-shaped mask within which the spectra and line centres are extracted. In all PV diagram panels, the grey box in the bottom-left corner has horizontal width</figcaption></figure>

equal to the beam major axis and vertical height equal to the channel width.
**a,b**, ALMA observations of the AB Aur disk in $^{13}$CO (a) and $^{18}$O (b). **c**, Synthetic
ALMA $^{13}$CO observations generated from 3D SPH simulations of a gravitationally
unstable disk with a disk-to-star mass ratio of 0.3 and a cooling rate described
by $\beta = 10$. **d**, The simulated disk has its velocity structure artificially post-
processed to be Keplerian.
<figure><img src="image_6.png" /></figure>

**Fig. 4 |** *PV wiggle morphology, magnitude and constraints on the AB Aur disk mass. a, The ALMA $^{13}$CO and C$^{18}$O line centres along the southern minor axis from Fig. 3a,b, after quadratic detrending (see Methods). Uncertainties on the line centres are shown by yellow-shaded regions. For qualitative comparison, the minor-axis PV wiggle of the analytic GI model disk from Fig. 2b is shown in the background in light grey. b, The magnitude of the southern minor-axis PV wiggle in AB Aur is measured to be 37.4 ± 2.9 m s$^{-1}$ in $^{13}$CO and 44.2 ± 1.3 m s$^{-1}$ in C$^{18}$O. c, Map of the minor-axis PV wiggle magnitude of 3,600 analytic GI model*

disks, calculated for a 60 × 60 grid of disk-to-star mass ratios and cooling timescales. Each cell in the map represents the minor-axis PV wiggle magnitude from a different model. A yellow contour is drawn at each of the AB Aur $^{13}$CO and C$^{18}$O measured magnitude values and dashed lines represent the quoted uncertainties. White-shaded regions denote $1\sigma$, $2\sigma$ and $3\sigma$ departures from a Keplerian signal in $^{13}$CO (see Fig. 3d). Horizontal bars indicate independently derived $\beta$ ranges at a selection of $M_{\text{disk}}/M_*$ values (see Methods).

Quantifying the minor-axis PV wiggle magnitude as above, we perform comparisons against analytic models to identify the combinations of disk mass ($M_{\text{disk}}/M_*$) and cooling timescale ($\beta$) that satisfy the AB Aur observations. A proof of concept of this technique with the SPH simulation is shown in Extended Data Fig. 7b. Using the analytic modelling code giggle (https://doi.org/10.5281/zenodo.10205110) from ref. 23 (Methods), we calculate the minor-axis PV wiggle magnitude in gravitationally unstable disk models for 60 × 60 combinations of $M_{\text{disk}}/M_*$ and $\beta$, letting each vary within the ranges $0 \le M_{\text{disk}}/M_* \le 0.4$ and $10^{-2} \le \beta \le 10^2$. A demonstrative analytic curve for the minor-axis PV wiggle from the same model shown in Fig. 2b is underlaid in Fig. 4a for qualitative comparison. Figure 4c shows the resulting map of 60 × 60 analytic minor-axis PV wiggle magnitudes. Overlaying contours in this map at the magnitude values measured for the AB Aur $^{13}$CO and C$^{18}$O southern minor-axis PV wiggles, we find a disk mass in the gravitationally unstable regime with $0.1 \le M_{\text{disk}}/M_* \le 0.3$ for a cooling timescale of $0.1 < \beta < 10$. This result is robust to plausible variations in the analytic model parameter choices (Extended Data Fig. 8). This disk mass range is broadly consistent with the observed spiral morphology—a lower disk mass may result in a large number of more tightly wound spirals than we observe, and vice versa$^{12,43}$. To demonstrate that the implied cooling timescales are compatible with the constrained disk-mass values, Fig. 4c also shows the ranges of $\beta$ derived from independent radiative cooling prescriptions (see Methods).

The detection of GI in the disk around AB Aur, a Class II young stellar object$^{32,33}$, demonstrates that GI can take place during later evolutionary stages. This result, together with previous reports of several protoplanet candidates in and among spiral arms in the system$^{35,44–46}$ (Extended Data Fig. 9), provides a direct observational connection between GI and planet formation. Looking forward, the AB Aur system can be an ideal test bed for understanding how planet formation is facilitated by GI-induced spiral arms—whether by fragmentation into gas clumps enabled by rapid cooling$^{7–10}$ ($\beta \le 3$) or by dust collapse of solids concentrated within spiral arms sustained by slow cooling$^{47–50}$ ($\beta \ge 5$).

## Online content

Any methods, additional references, Nature Portfolio reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/s41586-024-07877-0.

1. Chiang, E. & Youdin, A. N. Forming planetesimals in solar and extrasolar nebulae. Annu. Rev. Earth Planet. Sci. **38**, 493–522 (2010).

2. Johansen, A. & Lambrechts, M. Forming planets via pebble accretion. Annu. Rev. Earth Planet. Sci. **45**, 359–387 (2017).

3. Ormel, C. W. Formation, Evolution, and Dynamics of Young Solar Systems. Astrophysics and Space Science Library Vol. 445 (eds Pessah, M. & Gressel, O.) 197–228 (Springer, 2017).

4. Liu, B. & Ji, J. A tale of planet formation: from dust to planets. Res. Astron. Astrophys. **20**, 164 (2020).

5. Drążkowska, J. et al. Planet formation theory in the era of ALMA and Kepler: from pebbles to exoplanets. In Protostars and Planets VII Vol. 534 of the Astronomical Society of the Pacific Conference Series (eds Inutsuka, S. et al.) 717 (ASP, 2023).

6. Boss, A. P. Giant planet formation by gravitational instability. Science **276**, 1836–1839 (1997).

7. Gammie, C. F. Nonlinear outcome of gravitational instability in cooling, gaseous disks. Astrophys. J. **553**, 174–183 (2001).

8. Rice, W. K. M. et al. Substellar companions and isolated planetary-mass objects from protostellar disc fragmentation. Mon. Not. R. Astron. Soc. **346**, L36–L40 (2003).

9. Zhu, Z., Hartmann, L., Nelson, R. P. & Gammie, C. F. Challenges in forming planets by gravitational instability: disk irradiation and clump migration, accretion, and tidal destruction. Astrophys. J. **746**, 110 (2012).

10. Deng, H., Mayer, L. & Helled, R. Formation of intermediate-mass planets via magnetically controlled disk fragmentation. Nat. Astron. **5**, 440–444 (2021).

11. Cadman, J., Rice, K. & Hall, C. AB Aurigae: possible evidence of planet formation through the gravitational instability. Mon. Not. R. Astron. Soc. **504**, 2877–2888 (2021).

12. Lodato, G. & Rice, W. K. M. Testing the locality of transport in self-gravitating accretion discs. Mon. Not. R. Astron. Soc. **351**, 630–642 (2004).

13. Cossins, P., Lodato, G. & Clarke, C. J. Characterizing the gravitational instability in cooling accretion discs. Mon. Not. R. Astron. Soc. **393**, 1157–1173 (2009).

14. Dipierro, G., Lodato, G., Testi, L. & de Gregorio Monsalvo, I. How to detect the signatures of self-gravitating circumstellar discs with the Atacama Large Millimeter/sub-millimeter Array. Mon. Not. R. Astron. Soc. **444**, 1919–1929 (2014).

15. Kratter, K. & Lodato, G. Gravitational instabilities in circumstellar disks. Annu. Rev. Astron. Astrophys. **54**, 271–311 (2016).

16. Dong, R., Hall, C., Rice, K. & Chiang, E. Spiral arms in gravitationally unstable protoplanetary disks as imaged in scattered light. Astrophys. J. Lett. **812**, L32 (2015).

17. Hall, C. et al. Directly observing continuum emission from self-gravitating spiral waves. Mon. Not. R. Astron. Soc. **458**, 306–318 (2016).
18. Hall, C. et al. The temporal requirements of directly observing self-gravitating spiral waves in protoplanetary disks with ALMA. Astrophys. J. 871, 228 (2019).

19. Paneque-Carreño, T. et al. Spiral arms and a massive dust disk with non-Keplerian kinematics: possible evidence for gravitational instability in the disk of Elias 2-27. Astrophys. J. 914, 88 (2021).

20. Veronesi, B. et al. A dynamical measurement of the disk mass in Elias 227. Astrophys. J. Lett. 914, L27 (2021).

21. Stapper, L. M. et al. Constraining the gas mass of Herbig disks using CO isotopologues. Astron. Astrophys. 682, A149 (2024).

22. Hall, C. et al. Predicting the kinematic evidence of gravitational instability. Astrophys. J. 904, 148 (2020).

23. Longarini, C. et al. Investigating protoplanetary disk cooling through kinematics: analytical GI wiggle. Astrophys. J. Lett. 920, L41 (2021).

24. Terry, J. P. et al. Constraining protoplanetary disc mass using the GI wiggle. Mon. Not. R. Astron. Soc. 510, 1671-1679 (2022).

25. van den Ancker, M. E. et al. HIPPARCOS data on Herbig Ae/Be stars: an evolutionary scenario. Astron. Astrophys. 324, L33-L36 (1997).

26. DeWarf, L. E., Sepinsky, J. F., Guinan, E. F., Ribas, I. & Nadalin, I. Intrinsic properties of the young stellar object SU Aurigae. Astrophys. J. 590, 357-367 (2003).

27. Beck, T. L. & Bary, J. S. A search for spatially resolved infrared rovibrational molecular hydrogen emission from the disks of young stars. Astrophys. J. 884, 159 (2019).

28. Garufi, A. et al. The SPHERE view of the Taurus star-forming region. Astron. Astrophys. 685, A53 (2024).

29. Rodríguez, L. F. et al. An ionized outflow from AB Aur, a Herbig Ae Star with a transitional disk. Astrophys. J. Lett. 793, L21 (2014).

30. Guzmán-Diaz, J. et al. Homogeneous study of Herbig Ae/Be stars from spectral energy distributions and Gaia EDR3. Astron. Astrophys. 650, A182 (2021).

31. Gaia Collaboration. Gaia Data Release 3. Summary of the content and survey properties. Astron. Astrophys. 674, A1 (2023).

32. Henning, T., Burkert, A., Launhardt, R., Leinert, C. & Stecklum, B. Infrared imaging and millimetre continuum mapping of Herbig Ae/Be and FU Orionis stars. Astron. Astrophys. 336, 565-586 (1998).

33. Bouwman, J., de Koter, A., van den Ancker, M. E. & Waters, L. B. F. M. The composition of the circumstellar dust around the Herbig Ae stars AB Aur and HD 163296. Astron. Astrophys. 360, 213-226 (2000).

34. Pérez, L. M. et al. Spiral density waves in a young protoplanetary disk. Science 353, 1519-1521 (2016).

35. Boccaletti, A. et al. Possible evidence of ongoing planet formation in AB Aurigae. A showcase of the SPHERE/ALMA synergy. Astron. Astrophys. 637, L5 (2020).

36. Dong, R., Vorobyov, E., Pavlyuchenkov, Y., Chiang, E. & Liu, H. B. Signatures of gravitational instability in resolved images of protostellar disks. Astrophys. J. 823, 141 (2016).

37. Hashimoto, J. et al. Direct imaging of fine structures in giant planet-forming regions of the protoplanetary disk around AB Aurigae. Astrophys. J. Lett. 729, L17 (2011).

38. Fukagawa, M. et al. Spiral structure in the circumstellar disk around AB Aurigae. Astrophys. J. Lett. 605, L53-L56 (2004).

39. Lin, S.-Y. et al. Possible molecular spiral arms in the protoplanetary disk of AB Aurigae. Astrophys. J. 645, 1297-1304 (2006).

40. Perrin, M. D. et al. The case of AB Aurigae's disk in polarized light: is there truly a gap? Astrophys. J. Lett. 707, L132-L136 (2009).

41. Teague, R. & Foreman-Mackey, D. A robust method to measure centroids of spectral lines. Res. Notes AAS 2, 173 (2018).

42. Teague, R. Statistical uncertainties in moment maps of line emission. Res. Notes AAS 3, 74 (2019).

43. Lodato, G. & Rice, W. K. M. Testing the locality of transport in self-gravitating accretion discs — II. The massive disc case. Mon. Not. R. Astron. Soc. 358, 1489-1500 (2005).

44. Oppenheimer, B. R. et al. The solar-system-scale disk around AB Aurigae. Astrophys. J. 679, 1574-1581 (2008).

45. Tang, Y.-W. et al. Planet formation in AB Aurigae: imaging of the inner gaseous spirals observed inside the dust cavity. Astrophys. J. 840, 32 (2017).

46. Currie, T. et al. Images of embedded Jovian planet formation at a wide separation around AB Aurigae. Nat. Astron. 6, 751-759 (2022).

47. Rice, W. K. M., Lodato, G., Pringle, J. E., Armitage, P. J. & Bonnell, I. A. Accelerated planetesimal growth in self-gravitating protoplanetary discs. Mon. Not. R. Astron. Soc. 355, 543-552 (2004).

48. Longarini, C., Armitage, P. J., Lodato, G., Price, D. J. & Ceppi, S. The role of the drag force in the gravitational stability of dusty planet-forming disc - II. Numerical simulations. Mon. Not. R. Astron. Soc. 522, 6217-6235 (2023).

49. Booth, R. A. & Clarke, C. J. Collision velocity of dust grains in self-gravitating protoplanetary discs. Mon. Not. R. Astron. Soc. 458, 2676-2693 (2016).

50. Rowther, S. et al. The role of drag and gravity on dust concentration in a gravitationally unstable disc. Mon. Not. R. Astron. Soc. 528, 2490-2500 (2024).

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Springer Nature or its licensor (e.g. a society or other partner) holds exclusive rights to this article under a publishing agreement with the author(s) or other rightsholder(s); author self-archiving of the accepted manuscript version of this article is solely governed by the terms of such publishing agreement and applicable law.

© The Author(s), under exclusive licence to Springer Nature Limited 2024
# Methods

## Further information on the source

AB Aur is accreting from the disk at a rate of $\dot{M} \approx 10^{-7} M_{\odot} \text{yr}^{-1}$ (ref. 51), within the range expected for modest GI-driven accretion ($10^{-7} - 10^{-6} M_{\odot} \text{yr}^{-1}$)$^{52}$. This accretion rate, taken together with the current age $t_0 = 2.5-4.4$ Myr (refs. 25-28), implies a high 'latent disk mass': $M_{\text{disk}}^{\text{latent}} = \dot{M}(t_0) \times t_0 = 0.25-0.44 M_{\odot}$, or $M_{\text{disk}}^{\text{latent}}/M_{\star} \approx 0.1-0.2$. $M_{\text{disk}}^{\text{latent}}$ provides an accretion-rate-based assessment of disk mass, assuming a constant stellar accretion rate $\dot{M}$ and we are observing the system midway through the lifetime of the disk$^{53,54}$. This is a conservative estimate as the accretion rate at earlier epochs is probably higher$^{55}$. In millimetre continuum observations, the disk shows a dust ring at about 1" and a cavity inside$^{56}$, probably caused by the trapping of millimetre-sized dust at a pressure bump. The dust ring is located inside the main spirals in both the scattered light and gas emission. Late infall from above or below the main disk plane$^{56-59}$ is probably encouraging GI by providing a source of mass to maintain a high $M_{\text{disk}}/M_{\star}$ value$^{18,38}$.

## ALMA observations

We observed AB Aur with ALMA in April, May and September 2022 under ALMA programme ID 2021.1.00690.S (PI: R. Dong). Measurements were taken with the Band 6 receivers$^{60}$ in array configurations C-3 (two execution blocks) and C-6 (six execution blocks). In total, the eight execution blocks reached an on-source integration time of 5.75 h, making this the longest fine-kinematics ($v_{\text{chan}} < 100 \text{ m s}^{-1}$) programme towards a single protoplanetary disk so far. Extended Data Table 1 provides details of the observations. We centred one spectral window (SPW) at the $^{13}\text{CO}/J=2-1$ molecular emission line transition rest frequency (220.3986 GHz), covering a bandwidth of 58.594 MHz with 1,920 channels, resulting in the highest achievable spectral resolution of 41.510 m s$^{-1}$ after default spectral averaging with $N=2$ by Hanning smoothing within the correlator data processor. A second SPW was centred at the $^{18}\text{O}/J=2-1$ rest frequency (219.5603 GHz), covering the same bandwidth with half as many channels (960 channels; owing to sharing a baseband with another SPW), achieving a 83.336 m s$^{-1}$ spectral resolution. To enable self-calibration, our correlator setup sampled the continuum in another SPW centred at 233.012 GHz with 128 channels each 15.625 MHz in width, obtaining the full available 2.0-GHz bandwidth. Using the continuum data, all execution blocks were aligned to a common phase centre in the $u-v$ plane. We performed a series of phase-only self-calibration iterations and avoided combining by SPW in the first two rounds to remove any potential per-SPW phase offsets. We also carried out one round of amplitude and phase self-calibration. Finally, we applied the phase-centre realignments and calibration-gain tables (that we generated with the continuum data) to the line data. We performed continuum subtraction in the $u-v$ plane using the uvcontsub task.

All imaging was performed with the CASA tclean task. We used the multiscale deconvolution algorithm$^{61}$ with (Gaussian) deconvolution scales [0.02", 0.1", 0.3", 0.6", 1.0"). We did not image with a Keplerian mask so as not to restrict our ability to observe non-Keplerian emission. After experimentation with CASA's auto-multithresh masking algorithm$^{62}$, we used an imaging strategy similar to PHANGS-ALMA$^{63}$, in which we clean conservatively, with a broad mask (usemask='pb' and pbmask=0.2), forcing frequent major cycles. (The $^{13}\text{CO}$ robust 0.5 cube underwent 198 major cycles and the $^{18}\text{O}$ cube underwent 76.) To achieve frequent major cycles, we set the maximum number of minor-cycle iterations per channel to cycleniter=80, the minor cycle threshold to max_psf_sidelobe_level=3.0 and minpsffraction=0.5, and the maximum assigned clean component to gain=0.2 times the peak residual. We used a Briggs robust weighting scheme and generated two sets of image cubes; one with a robust value of 0.5 and a second with robust 1.5. The corresponding beam sizes for $^{13}\text{CO}$ are 237 × 175 mas, 1.2° for robust 0.5 and 390 × 274 mas, -1.4° for robust 1.5. We imaged

with a field of view out to the primary beam full width at half maximum (FWHM; 38") with 0.02" pixels (9 or 12 pixels per synthesized beam minor or major axis, respectively). We imaged in LSRK velocity channels at 42 m s$^{-1}$ for $^{13}\text{CO}$ and 84 m s$^{-1}$ for $^{18}\text{O}$, respectively (nearly native channel spacing). The CLEAN threshold was set to five times the r.m.s. noise measured in 20 line-free channels of the dirty image cube. We applied JyM correction$^{64,65}$ and primary-beam correction. The r.m.s. noise in the resulting $^{13}\text{CO}$ cubes imaged with robust 0.5 and robust 1.5 is 2.0 mJy per beam and 1.2 mJy per beam, respectively, and 0.6 mJy per beam in the $^{18}\text{O}$ cube imaged with robust 1.5.

We used the robust 0.5 image cubes for our position-position analysis (moment maps; Figs. 1 and 2) and the robust 1.5 cubes for our PV analysis (PV diagrams and line centres; Figs. 3 and 4). We made the moment 0, 1 and 2 maps using the bettermoments$^{41,66}$ methods collapse zeroth, collapse first and collapse percentiles, respectively. We note that we calculate our 'moment 2' maps as the average of the redshifted and blueshifted linewidths about the intensity-weighted median line centre (that is, as the average of the wpdVr and wpdVb maps). Mathematically, this is a different approach to find the linewidth than the classic moment 2 approach, although in our case, we find that the two yield nearly identical outcomes. We applied sigma clipping at five times the r.m.s. noise and performed no spectral smoothing.

## Geometric properties

We used the Python package eddy$^{67}$ to infer geometric properties of the disk, namely to constrain the disk centre $x_0, y_0$, the disk inclination $i$, the PA, the $v_{\text{sys}}$ and the dynamical stellar mass $M_{\star}$. We performed Markov chain Monte Carlo to fit the $^{18}\text{O}$ moment 1 map (Extended Data Fig. 4) with a geometrically thin Keplerian disk rotation profile:

$$v_0 = \sqrt{\frac{GM_*}{r}} \sin i \cos\phi + v_{\text{sys}} \quad (1)$$

in which $r$ is the disk radius, $\phi$ is the azimuthal angle around the disk and $G$ is the gravitational constant. Following convention, we fix the inclination to the value found from fitting the continuum, $i=23.2^{\circ}$ (refs. 45,56), and the distance to 155.9 pc (refs. 31,68). We assumed flat priors for all values and spatially downsampled the rotation map to the beam FWHM before the likelihood calculation so that only spatially independent pixels were considered. The calculation of the posterior distributions was run with 128 walkers and an initial burn-in period of 10,000 steps before the posterior distributions were sampled for a further 10,000 steps. The resulting posterior distributions were $x_0 = -5 \pm 7$ mas, $y_0 = -17 \pm 7$ mas, PA = 236.7 $\pm$ 0.3$^{\circ}$, $M_* = 2.23 \pm 0.02 M_{\odot}$ and $v_{\text{sys}} = 5,858 \pm 5 \text{ m s}^{-1}$, for which we report the uncertainties represented by the 16th and 84th percentiles about the median value. The last three values are consistent with constraints from previous observations$^{26,27,45,56,69}$.

## Hydrodynamic simulations and synthetic ALMA observations

We performed 3D global SPH simulations with the PHANTOM code$^{70}$ using 1 million SPH particles. We assumed a central star mass of $2.4 M_{\odot}$ (refs. 26,27), represented by a sink particle$^{71}$ with accretion radius set to 60 AU. The initial inner and outer disk radii were set to $r_{\text{in,SPH}} = 80$ AU and $r_{\text{out,SPH}} = 500$ AU, respectively. We set the initial gas mass to $0.7 M_{\odot}$, corresponding to $M_{\text{disk}}/M_{\star} = 0.29$. The surface density profile follows $\Sigma \propto r^{-p}$ (in which the power-law index $p=1.0$) and the sound speed profile follows $c_s \propto r^{-q}$ (in which $q=0.25$). The initial disk aspect ratio was set to $H/r = 0.05$ at 80 AU. We set $\alpha_{\text{SPH}}$ such that $\alpha_{\text{min}} \le \alpha_{\text{SPH}} \le \alpha_{\text{max}}$, with $\alpha_{\text{min}} = 0.001$ and $\alpha_{\text{max}} = 1.0$, with the value of $\alpha_{\text{SPH}}$ set by the Cullen and Dehnen$^{72}$ switch that increases viscosity only in the case of converging flows. This results in a Shakura-Sunyaev viscosity of $\alpha_{\text{SS}} \approx 0.01$ throughout the disk.

We assumed an adiabatic equation of state, with heating from compressional P dV work and shock heating. The disk cools by Gammie
cooling⁷ (also known as β cooling), for which the cooling timescale is proportional to the local dynamical time by the factor $β$, such that $t_{\text{cool}}(r) = \beta\Omega^{-1}(r)$, in which $\Omega(r) = (GM_*/r^3)^{1/2}$ is the Keplerian frequency. We set $\beta = 10$, a typical value used or found in simulations¹⁹,²²⁻²⁴. We let the simulation evolve for five orbital periods of the outermost particle, at which point the disk settles into a state in which the Toomre $Q$ parameter is between 1 and 2 from $r_{\text{in,SPH}}$ to $1.1r_{\text{out,SPH}}$.

We computed the disk thermal structure and ¹³CO ($J=2-1$) model line cubes using the Monte Carlo radiative transfer code MCFOST⁷,³,⁴. We assumed that the ¹³CO molecule is in local thermodynamic equilibrium with its surroundings and that the dust is in thermal equilibrium with the gas ($T_{\text{gas}} = T_{\text{dust}}$). We set the ¹³CO/H₂ abundance to $7 \times 10^{-7}$ (refs. 22,75) and we used approximately $10^7$ photon packets to calculate $T_{\text{dust}}$. We performed Voronoi tessellation on 990,972 SPH particles, which corresponded to 99% of the mass in the simulation. We set the total dust mass to 1% of the total SPH gas mass and used a dust-grain population with 50 logarithmic bins ranging in size from 0.1 µm to 3.0 mm. The dust optical properties are computed using Mie theory. The central star was represented as a sphere of radius 2.5 $R_\odot$ radiating isotropically at an effective temperature $T_{\text{eff}} = 9,770$ K, set to match AB Aur⁴,⁶,⁷,⁸. The disk was given an inclination of 23.2°, a PA of 236.7° (for which PA is measured east of north to the redshifted major axis) and placed at a distance of 155.9 pc, all consistent with the AB Aur system.

We used the same PHANTOM simulation to create both the GI and Keplerian model line cubes shown in Figs. 2c and 3. We created the Keplerian counterpart with MCFOST, using the flags -no vr and -no vz to force the radial and vertical velocities to be zero, and -vphi Kep to force the azimuthal velocities to be Keplerian. Both ¹³CO model line cubes were generated with MCFOST, binned at the observed spectral resolution of 42 m s⁻¹ and gridded in the image plane to have 2,048 × 2,048 pixels of size 0.02". We assumed a turbulent velocity of 0.05 km s⁻¹.

We generated synthetic ALMA image cubes from the ¹³CO model line cubes using syndisk (https://github.com/richteague/syndisk) to match the properties of the observed AB Aur ¹³CO image cubes (robust 0.5 and 1.5). In the latter case, the model line cube was convolved with a beam of size 0.390" × 0.274" and PA -1.4°. Correlated noise was added with an r.m.s. of 1.2 mJy per beam. The model data were then smoothed with a Hanning spectral response function with a resolution of 42 m s⁻¹. Effects associated with interferometric or spatial filtering are not captured by this process and our synthetic ALMA image cubes are effectively fully sampled in the u-v plane. The synthetic cubes were collapsed into moment maps following the same procedure as the AB Aur data (Extended Data Fig. 1).

## Analytic modelling

We analytically compute the velocity fields of gravitationally unstable disks using the giggle (https://doi.org/10.5281/zenodo.10205110) package developed by Longarini et al.²³. Working in 2D polar coordinates $(r, \phi)$, giggle considers a geometrically thin disk with surface density profile $\Sigma_0 \propto r^{-p}$ and inclination $i$, centred on a star of mass $M_*$. It computes the projected line-of-sight velocity field as:

$$v_{\text{los}} = (v_r \sin\phi + v_\phi \cos\phi) \sin i + v_{\text{sys}}, \quad (2)$$

in which $v_r$ and $v_\phi$ are the radial and azimuthal components of the disk-velocity field, respectively. The basic state of the disk (that is, considering only the gravitational potential contribution from the central star) is assumed to be Keplerian: $v_r = 0$ and $v_\phi = v_{\text{kep}}$. The scheme of the model is to determine the perturbations in $v_r$ and $v_\phi$ generated by GI by taking into account the extra gravitational contribution from the disk, which is initialized as marginally unstable and imprinted with global spiral density perturbations. The model computes the velocity field under the assumption that the disk is self-regulated. This state is imposed by assuming a balance between heating (by compression and shocks within the spiral arms) and cooling (by radiative processes).

As such, the amplitude of the spiral density perturbations $A_{\Sigma_{\text{spir}}}/\Sigma_0$ is saturated to a finite value proportional to the cooling timescale $\beta$ (refs. 13,79) is:

$$\frac{A_{\Sigma_{\text{spir}}}}{\Sigma_0} = \chi \beta^{-1/2}, \quad (3)$$

in which the proportionality factor $\chi$ is of order unity¹³,²³. The imprinted spiral density perturbation is assumed to be small relative to the background surface density, so that all the relevant quantities (density $\Sigma$, gravitational potential $\Phi$, velocities $v_r$ and $v_\phi$ and enthalpy $h$) can be written as a linear sum of the basic state and the perturbation:

$$X(r, \phi) = X_0(r) + X_{\text{spir}}(r, \phi). \quad (4)$$

The spiral perturbation in density is given the form:

$$\Sigma_{\text{spir}}(r, \phi) = \Re[A_{\Sigma_{\text{spir}}} e^{j(m\phi + \psi(r))}], \quad (5)$$

in which $j = \sqrt{-1}$ (as we are using $i$ to represent the disk inclination) and $m$ is the azimuthal wavenumber. The 'shape function' $\psi(r)$ is described by $m$ and the spiral pitch angle $\alpha_{\text{pitch}}$ as:

$$\psi(r) = \frac{m}{\tan\alpha_{\text{pitch}}} \log r, \quad (6)$$

which is related to the radial wavenumber $k$ by $d\psi/dr = k$. The spiral density perturbation necessarily introduces a corresponding perturbation to the gravitational potential:

$$\Phi_{\text{spir}}(r, \phi) = -\frac{2\pi G}{|k|} \Sigma_{\text{spir}}(r, \phi). \quad (7)$$

The negative proportionality $\Phi_{\text{spir}} \propto -\Sigma_{\text{spir}}$ is the definition of self-gravitating spiral arms. As a result, corresponding perturbations in the azimuthal and radial velocities are driven:

$$v_r(r, \phi) = \Re[A_{v_r}(r) e^{j(m\phi + \psi(r))}], \quad (8)$$

$$v_\phi(r, \phi) = \Re[A_{v_\phi}(r) e^{j(m\phi + \psi(r))} + r\Omega], \quad (9)$$

in which we note $r\Omega \neq v_{\text{kep}}$ because the angular frequency $\Omega$ includes super-Keplerian rotation from the disk mass contribution:

$$\Omega^2 = \frac{GM_*}{r^3} + \frac{1}{r} \frac{\partial \Phi_{\text{disk}}}{\partial r}. \quad (10)$$

By assuming that the disk is marginally unstable, and by maintaining the self-regulated state condition, the amplitude of the radial and azimuthal velocity perturbations $A_{v_r}(r)$ and $A_{v_\phi}(r)$ are determined²³:

$$A_{v_r}(r) = 2jm\chi\beta^{-1/2} \left( \frac{M_{\text{disk}}(r)}{M_*} \right)^2 v_{\text{kep}}(r), \quad (11)$$

$$A_{v_\phi}(r) = -\frac{1}{2}j\chi\beta^{-1/2} \left( \frac{M_{\text{disk}}(r)}{M_*} \right)^2 v_{\text{kep}}(r), \quad (12)$$

in which $M_{\text{disk}}(r)$ is the disk mass enclosed within radius $r$. With a surface density profile $\Sigma_0(r) \propto r^{-p}$, then $M_{\text{disk}}(r) \propto r^{p+2}$ and the amplitude of the radial perturbation is described by $A_{v_r}(r) \propto r^{-2p+7/2}$. For $p < 7/4$, $A_{v_r}(r)$ is an increasing function of radius. The factor of imaginary number $j$ in equation (11) has important physical consequences: when the real component of $A_{v_r}(r)$ is taken (equation (8)), the radial-velocity
perturbation is $\pi/2$ out of phase with the spiral density perturbation (equation (5)) and convergent at the locations at which $\Sigma_{\text{spir}}$ takes a maximum. Explicitly,

$$v_r(r, \frac{\pi}{2}) \propto -\sin(m\frac{\pi}{2} + \psi(r)), \qquad (13)$$

$$\Sigma_{\text{spir}}(r, \frac{\pi}{2}) \propto \cos(m\frac{\pi}{2} + \psi(r)). \qquad (14)$$

For qualitative visual comparison with the AB Aur moment 1 map in Fig. 2a, we compute the projected line-of-sight velocity field of a gravitationally unstable disk with $\beta = 10$ and $M_{\text{disk}}/M_* = 0.3$ in Fig. 2b. We set $m = 3$ and $\alpha_{\text{pitch}} = 15^\circ$ to approximately match the $^{13}\text{CO}$ spirals in the AB Aur disk (Fig. 1) and assume $p = 1.0$ and $\chi = 1.0$ (ref. 13). The dominant azimuthal wavenumber is expected to be inversely related to the disk-to-star mass ratio $q$, roughly obeying $m \approx 1/q$ (refs. 12,13,16), so our choice of $m = 3$ is consistent with $M_{\text{disk}}/M_* \approx 0.3$.

## Revealing global spiral structure

We obtain the residual moment maps shown in Fig. 1 using a variation on the conventional high-pass-filtering (also known as unsharp masking) technique. The conventional method is to convolve the image with a Gaussian kernel and subtract the blurred image from the original. It is a common technique to increase the visual contrast of variations in an image and has been used successfully to reveal spiral structure disks (for example, refs. 28,34,35,80–83). Here we perform the convolution with a radially expanding kernel (https://github.com/jjspeedie/ expanding_kernel)—that is, with a Gaussian kernel whose FWHM, $w$, increases with radial distance from the image centre with a simple power-law dependence:

$$w(r) = w_0 \times (r/r_0)^\gamma, \qquad (15)$$

in which $w_0$ is the kernel width at $r_0 = 1''$. A radially expanding kernel provides a way to highlight variations more evenly throughout the disk, given the spatial scales of the variations (which are expected to track with the local scale height and increase with radius) and the dynamical range of the variations, which fall with radius. After experimentation, we use $w_0 = 0.3''$ and $\gamma = 0.25$, although we emphasize that this is a qualitative choice and the key spiral features, such as their locations, are robust against a variety of choices in kernel parameters. The high-pass-filter technique is also flexible to the disk emission surface morphology and can capture global-scale deviations from Keplerian rotation in the background disk. Extended Data Fig. 4 compares the residual moment 1 maps in $^{13}\text{CO}$ and $C^{18}\text{O}$ obtained after subtracting the axisymmetric geometrically thin Keplerian model (equation (1)) versus after subtracting a blurred version of the moment 1 map made with the expanding kernel filter. The Keplerian residuals (panels c and h) show signs of global-scale deviation from Keplerian: the east (west) side is generally blueshifted (redshifted), hinting at super-Keplerian rotation, signatures of disk mass contributing to the total mass of the system. Although spiral structure is indeed also visible in the Keplerian residuals, the expanding kernel residuals (panels e and j) reveal the underlying spiral structure in a spatially even manner, indicating that the expanding kernel background model (panels d and i) more successfully captures the quasi-local background disk velocity. We note that this background model is non-axisymmetric; it shows excess blueshifted velocity in the southeast quadrant of the disk such that the contour of $v_{\text{los}} = v_{\text{sys}}$ diverges westward from the minor axis south of the star, possibly indicative of a global disk warp. This is what necessitates a detrending of the line centres to isolate the sinusoidal component of the southern minor-axis PV wiggle in Fig. 4a (see the section 'Measuring the magnitude of the minor-axis PV wiggle'). Filtered moment maps for the synthetic ALMA observations of the simulated SPH GI disk are shown in Extended Data Fig. 2.

## Global kinematics of self-gravitating spiral arms

Radially convergent motion (as in Fig. 1b–d insets) serves as a kinematic signature for the location of self-gravitating spiral arms at disk azimuths at which the radial-velocity perturbation contributes sufficiently strongly to the observed velocity field and thus cannot be a fully unambiguous locator at disk azimuths away from the minor axis. Extended Data Fig. 5c,g provides maps of velocity residuals from Keplerian for the 2D analytic GI disk model and the SPH GI disk simulation. The convergent motion towards the spiral spines is visible for a range of azimuths around the minor axis but becomes progressively less clear moving towards the major axis as the azimuthal velocity (super-Keplerian rotation) contributes progressively more to the line of sight. However, high-pass filtering (panel h) captures and removes the background super-Keplerian rotation, leaving a residual map that resembles the isolated radial component (panel d). Extended Data Fig. 5i–l overlays the locations of $^{13}\text{CO}$ spirals in the AB Aur disk (from filtered moment 0/2; Fig. 1c,d) onto the filtered moment 1 maps, to illustrate where convergent motion does or does not serve as a locator throughout the disk. Ambiguity occurs around the major axis, which is a location of transition in the sign of $v_r \sin i \sin\phi$ (first term of equation (2)), and when two spirals are not well separated and their motions superimpose. Three of the seven spiral structures in VLT/SPHERE scattered light seem to be spatially associable with those in $^{13}\text{CO}$ (S1, S5 and S7; panel l inset). Offsets in the southeast quadrant of the disk (S2, S3 and S4) may be further indication of a disk warp (Extended Data Fig. 4d,i) or other non-trivial phenomena (for example, vertical density and temperature gradients, or projection effects$^{84}$).

The kinematic signatures observed in the present ALMA dataset—probing disk scales from about 100 to 1,000 AU—are recognizably different from what is expected for planet-driven perturbations. Planetary wakes are dampened and become nearly circular as they propagate away from the planet$^{85–87}$, whereas GI-driven spirals maintain their modest pitch angles with radius and the amplitude of the induced velocity perturbations depends on the enclosed disk mass (equations (11) and (12)). In the planetary case, the density and radial-velocity perturbations are in phase (their peaks spatially coincide) and the pattern of motion within an arm along a radial cross-section is divergent$^{88,89}$. Overall, the essential characteristic of GI-induced spirals is that they occur globally$^{22,23}$ (see Fig. 1 and Extended Data Figs. 2, 3 and 5). In previous datasets investigating smaller spatial scales (within the central cavity of the AB Aur disk) planetary candidates P1/f1 (refs. 35,45), P2/b (refs. 45,46,90–92) and f2 (ref. 35) are known to be associated with—or driving—spiral arms, as observed in VLT/SPHERE scattered light and/or ALMA $^{12}\text{CO}$ emission. As shown in Extended Data Fig. 9, owing to their small separations ($\le 0.7''$), kinematic signatures from these candidates are inaccessible to our ALMA observations. Clump-like signals ‘c’ and ‘d’ seen by HST/STIS (ref. 46) at wide separations (approximately 2.75” and approximately 3.72”, respectively) are in locations tentatively suggestive of constituting spiral arm fragments and may warrant further investigation.

## PV analysis

We use the robust 1.5 image cubes for our PV analysis to maximize the recovery of emission at large disk radii. Owing to the clear association with a self-gravitating spiral arm (Fig. 1b–d insets), we target the wiggle on the southern minor axis. A clear spiral arm in moment 0/2 crossing the northern minor axis is also observed, but at the outer edge of the recovered $^{13}\text{CO}$ and $C^{18}\text{O}$ emission (about 3’’; see Extended Data Fig. 5k,l). We obtain the PV diagrams shown in Fig. 3 using eddy$^{67}$ to extract spectra from pixels within a 0.5°-wide wedge-shaped mask oriented 90° clockwise of the redshifted major axis (shown in Fig. 3 insets). Our quantitative analysis of the minor-axis PV wiggles is performed with maps of the line centres made using the quadratic method of bettermoments$^{41,66}$, which fits a quadratic curve to the spectrum in each pixel of the cube:
$$I(\nu) = a_0 + a_1(\nu - \nu_{\text{peak}}) + a_2(\nu - \nu_{\text{peak}})^2, \quad (16)$$

in which $\nu_{\text{peak}}$ is the channel of peak intensity in the spectrum. We select this approach over the traditional intensity-weighted mean velocity (moment 1) method specifically for its ability to provide well-characterized, statistically meaningful uncertainties on the line centre, $\sigma_{v\text{los}}$ (ref. 41). The statistical uncertainty on each line centre is computed as:

$$\sigma_{v\text{los}} = \sqrt{\frac{\sigma_I^2}{8} \left( \frac{3}{a_2^2} + \frac{a_1^2}{a_2^4} \right)}, \quad (17)$$

in which $\sigma_I$ is the r.m.s. noise of the intensities (see ref. 41 for a derivation). The quadratic method also has the advantage of being unaffected by sigma clipping and of automatically distinguishing the front side of the disk from the back side$^{41}$. Before the quadratic fitting, we spectrally smooth the data with a Savitzky-Golay filter of polynomial order 1 and filter window length of ten channels (420 m s$^{-1}$) in the case of $^{13}$CO and three channels (252 m s$^{-1}$) in the case of $^{18}$O. The former was also applied to the two synthetic ALMA $^{13}$CO image cubes generated from the SPH simulations. We extract the values from the resulting line-centre and line-uncertainty maps within the same wedge mask described above. The extracted line-centre values are shown as yellow points in Fig. 3 and the uncertainties are shown as yellow-shaded regions in Fig. 4a.

## Measuring the magnitude of the minor-axis PV wiggle

Following ref. 23, we measure the 'magnitude' of a minor-axis PV wiggle as the standard deviation of the line-centre values over a radial range. Bounded by the inner central cavity and the outer edge of $^{18}$O emission, we use a radial range of 1" to 5". We estimate the uncertainty on the magnitude measurement using a resampling procedure: we take 10,000 draws from Gaussian distributions centred on the observed line centres with standard deviation $\sigma_{v\text{los}}$ (equation (17)) to create 10,000 instances of the minor-axis PV wiggle, compute their magnitudes and then report the uncertainty as the standard deviation of those 10,000 magnitude estimates.

As well as the wiggle, the $^{13}$CO and $^{18}$O emissions on the southern disk minor axis also exhibit an underlying monotonic blueward trend with disk radius, seen in Fig. 3a,b as a subtle downward bend with radius of the line centres or equivalently in Fig. 2a as a westward or clockwise shift in the contour of $\nu_{v\text{los}} = \nu_{\text{sys}}$. We earmark this feature as a possible disk warp (Extended Data Fig. 4d,i) and use a least-squares-fitting approach to isolate the sinusoidal component of the PV wiggle. This approach yields the background trend line that minimizes the standard deviation of the residuals, thus providing the most conservative estimate for the magnitude of the detrended PV wiggle. We fit a quadratic trend line (Extended Data Fig. 6a) as it more closely resembles the high-pass-filter background curve than a linear one (Extended Data Fig. 6b,c). We show the quadratically detrended PV wiggles in Fig. 4a and report their magnitudes in Fig. 4b. We find very similar magnitudes for both the $^{13}$CO and $^{18}$O wiggles, despite $^{18}$O probably tracing lower optical depths in the AB Aur disk. This empirically substantiates comparisons with the 2D analytic model (next section).

Performing the same procedure outlined above on the synthetic $^{13}$CO minor-axis PV wiggle of the GI disk in the SPH simulation, we find a wiggle magnitude of 39.1 ± 1.9 m s$^{-1}$ (Extended Data Fig. 7).

## Constraining disk mass with quantitative comparisons with analytic models

We perform quantitative comparisons between the observed $^{13}$CO and $^{18}$O minor-axis PV wiggles and the projected radial-velocity component in our analytic model, $v_s\text{,sini}$ (ref. 23). From equations (8) and (11), the projected radial velocity on the minor axis ($\phi = \pi/2$) is:

$$v_r(r, \frac{\pi}{2}) \sin i = -2m\chi\beta^{-1/2} \left( \frac{M_{\text{disk}}(r)}{M_*} \right)^2 \nu_{\text{Kep}}(r) \sin \left( m\frac{\pi}{2} + \psi(r) \right) \sin i. \quad (18)$$

This curve reflects the disk mass enclosed within the inner and outer radii of the model, which we set to span the same projected radial range as the observed PV wiggles (1" to 5"). We compute 3,600 of these curves for a 60 × 60 grid of models with (total enclosed) $M_{\text{disk}}/M_*$ linearly spaced $\in [0, 0.4]$ and $\beta$ logarithmically spaced $\in [10^{-2}, 10^2]$. Again, we set $m = 3$ and $\alpha_{\text{pitch}} = 15^\circ$ to match the AB Aur disk and assume $p = 1.0$ and $\chi = 1.0$ (ref. 13). For qualitative comparison, we plot an example analytic minor-axis PV wiggle behind the data in Fig. 4a; the model has $\beta = 10$ and $M_{\text{disk}}/M_* = 0.3$. We show in Extended Data Fig. 8 that $m = 3$ reproduces the observed wiggles better than other choices and that $p = 1.5$ could also provide a satisfying match, whereas $p = 2.0$ is too steep. Because the wiggle amplitude is independent of $\alpha_{\text{pitch}}$ (equation (11)), the magnitude is constant with $\alpha_{\text{pitch}}$ when sampled over the same range in phase (not shown).

We measure the minor-axis PV wiggle magnitude of the 3,600 models and present the resulting magnitude map in Fig. 4c. By drawing contours in the Fig. 4c map at the magnitude values measured for AB Aur (37.4 ± 2.9 m s$^{-1}$ in $^{13}$CO and 44.2 ± 1.3 m s$^{-1}$ in $^{18}$O), we find every combination of $M_{\text{disk}}/M_*$ and $\beta$ that satisfies the observations. Repeating this procedure with our synthetic ALMA observations of the SPH GI disk simulation shown in Fig. 3c, we find that this technique successfully recovers the disk mass set in the underlying SPH simulation (Extended Data Fig. 7).

For independent physical estimates of plausible $\beta$ values between 1" and 5" (155 to 780 AU), we rely on radiative cooling prescriptions$^{93,94}$. From equation (39) in ref. 94, $\beta$ is a function of $r$ and depends on $M_{\text{disk}}$ through the surface density $\Sigma$. We assume $T = \left( \frac{\phi L_*}{8\pi r^2 \sigma_{\text{SB}}} \right)^{1/4}$, in which $\sigma_{\text{SB}}$ is the Stefan-Boltzmann constant, $L_* = 59 L_\odot$ is the stellar luminosity of AB Aur$^{46}$ and $\phi = 0.02$ represents the flaring angle$^{95}$. We use the DSHARP Rosseland mean opacity$^{96}$ $\kappa_R = \kappa_R(T, a_{\text{max}})$ for a power-law grain-size distribution truncated at $a_{\text{max}}$. We set $a_{\text{max}}$ to 0.1 mm and the dust-to-gas mass ratio to $f = 0.1\%$, based on radial drift arguments and lack of (sub-)millimetre emission at these large radii. We compute a $\beta(r)$ profile for each $M_{\text{disk}}/M_* \in [0, 0.4]$ and extract the values at 1" and 5". We overlay the resulting $\beta(M_{\text{disk}}/M_*)$ ranges as white-shaded regions in Extended Data Fig. 8 (in which the dependence on $p$ arises from the dependence on $\Sigma$) and in Fig. 4 as white horizontal bars at a selection of $M_{\text{disk}}/M_*$ values. For example, for $M_{\text{disk}}/M_* = 0.2$ and $p = 1.0$, we find $\beta(1') = 5.3$ and $\beta(5') = 3.6 \times 10^{-2}$. Although knowledge of cooling in disks is very limited, these estimates help to emphasize that not all values of $\beta$ are equally likely.

## Data availability

All observational data products presented in this work are available through the CANFAR Data Publication Service at https://doi.org/10.11570/24.0087. All simulated data products are available at https://doi.org/10.5281/zenodo.11668694. The raw ALMA data are publicly available at the ALMA archive (https://almascience.nrao.edu/aq/) under project ID 2021.1.00690.S. The raw VLT/SPHERE data are publicly available from the ESO Science Archive Facility (https://archive.eso.org/eso/eso_archive_main.html) under programme 0104.C-0157(B). Source data are provided with this paper.

## Code availability

ALMA data-reduction and imaging scripts are available at https://jjspeedie.github.io/guide.2021.1.00690.S. The following Python packages were used in this work: bettermoments (https://github.com/richteague/bettermoments), eddy (https://github.com/richteague/eddy),
51. Salyk, C. et al. Measuring protoplanetary disk accretion with H I Pfund β. Astrophys. J. 769, 21 (2013).

52. Rice, W. K. M. & Armitage, P. J. Time-dependent models of the structure and stability of self-gravitating protoplanetary discs. Mon. Not. R. Astron. Soc. 396, 2228-2236 (2009).

53. Hartmann, L., Calvet, N., Gullbring, E. & D'Alessio, P. Accretion and the evolution of T Tauri disks. Astrophys. J. 495, 385-400 (1998).

54. Dong, R., Najita, J. R. & Brittain, S. Spiral arms in disks: planets or gravitational instability? Astrophys. J. 862, 103 (2018).

55. Sicilia-Aguilar, A., Henning, T. & Hartmann, L. W. Accretion in evolved and transitional disks in CEP OB2: looking for the origin of the inner holes. Astrophys. J. 710, 597-612 (2010).

56. Tang, Y. W. et al. The circumstellar disk of AB Aurigae: evidence for envelope accretion at late stages of star formation? Astron. Astrophys. 547, A84 (2012).

57. Nakajima, T. & Golimowski, D. A. Coronagraphic imaging of pre-main-sequence stars: remnant envelopes of star formation seen in reflection. Astron. J. 109, 1181-1198 (1995).

58. Grady, C. A. et al. Hubble Space Telescope space telescope imaging spectrograph coronagraphic imaging of the Herbig AE star AB Aurigae. Astrophys. J. Lett. 523, L151-L154 (1999).

59. Rivière-Marichalar, P. et al. AB Aur, a Rosetta stone for studies of planet formation. I. Chemical study of a planet-forming disk. Astron. Astrophys. 642, A32 (2020).

60. Ediss, G. A. et al. in Proc. 15th International Symposium on Space Terahertz Technology (ed. Narayanan, G.) 181-188 (ISST, 2004).

61. Cornwell, T. J. Multiscale CLEAN deconvolution of radio synthesis images. IEEE J. Sel. Top. Signal Process. 2, 793-801 (2008).

62. Kepley, A. A. et al. Auto-multithresh: a general purpose automasking algorithm. Publ. Astron. Soc. Pac. 132, 024505 (2020).

63. Leroy, A. K. et al. PHANGS-ALMA data processing and pipeline. Astrophys. J. Suppl. Ser. 255, 19 (2021).

64. Jorsater, S. & van Moorsel, G. A. High resolution neutral hydrogen observations of the barred spiral galaxy NGC 1365. Astron. J. 110, 2037 (1995).

65. Czekala, I. A. et al. Molecules with ALMA at Planet-forming Scales (MAPS). II. CLEAN strategies for synthesizing images of molecular line emission in protoplanetary disks. Astrophys. J. Suppl. Ser. 257, 2 (2021).

66. Teague, R. & Foreman-Mackey, D. bettermoments: a robust method to measure line centroids. Zenodo https://doi.org/10.5281/zenodo.1419753 (2018).

67. Teague, R. eddy: extracting protoplanetary disk dynamics with Python. J. Open Source Softw. 4, 1220 (2019).

68. Gaia Collaboration. The Gaia mission. Astron. Astrophys. 595, A1 (2016).

69. Piétu, V., Guilloteau, S. & Dutrey, A. Sub-arcsec imaging of the AB Aur molecular disk and envelope at millimeter wavelengths: a non Keplerian disk. Astron. Astrophys. 443, 945-954 (2005).

70. Price, D. J. et al. Phantom: a smoothed particle hydrodynamics and magnetohydrodynamics code for astrophysics. Publ. Astron. Soc. Aust. 35, e031 (2018).

71. Bate, M. R., Bonnell, I. A. & Price, N. M. Modelling accretion in protobinary systems. Mon. Not. R. Astron. Soc. 277, 362-376 (1995).

72. Cullen, L. & Dehnen, W. Inviscid smoothed particle hydrodynamics. Mon. Not. R. Astron. Soc. 408, 669-683 (2010).

73. Pinte, C., Ménard, F., Duchêne, G. & Bastien, P. Monte Carlo radiative transfer in protoplanetary disks. Astron. Astrophys. 459, 797-804 (2006).

74. Pinte, C. et al. Benchmark problems for continuum radiative transfer. High optical depths, anisotropic scattering, and polarisation. Astron. Astrophys. 498, 967-980 (2009).

75. Pinte, C. et al. Kinematic evidence for an embedded protoplanet in a circumstellar disk. Astrophys. J. Lett. 860, L13 (2018).

76. Li, D. et al. An ordered magnetic field in the protoplanetary disk of AB Aur revealed by mid-infrared polarimetry. Astrophys. J. 832, 18 (2016).

77. Hillenbrand, L. A., Strom, S. E., Vrbka, F. J. & Keene, J. Herbig Ae/Be stars: intermediate-mass stars surrounded by massive circumstellar accretion disks. Astrophys. J. 397, 613-643 (1992).

78. Natta, A. et al. A reconsideration of disk properties in Herbig Ae stars. Astron. Astrophys. 371, 186-197 (2001).

79. Lodato, G. Classical disc physics. New Astron. Rev. 52, 21-41 (2008).

80. Rosotti, G. P. et al. Spiral arms in the protoplanetary disc HD100453 detected with ALMA: evidence for binary–disc interaction and a vertical temperature gradient. Mon. Not. R. Astron. Soc. 491, 1335-1347 (2020).

81. Meru, F. et al. On the origin of the spiral morphology in the Elias 2–27 circumstellar disk. Astrophys. J. Lett. 839, L24 (2017).

82. Zhang, Y. et al. Disk evolution study through imaging of nearby young stars (DESTINYS): diverse outcomes of binary–disk interactions. Astron. Astrophys. 672, A145 (2023).

83. Norfolk, B. J. et al. The origin of the Doppler flip in HD 100546: a large-scale spiral arm generated by an inner binary companion. Astrophys. J. Lett. 936, L4 (2022).

84. Ginski, C. et al. Direct detection of scattered light gaps in the transitional disk around HD 97048 with VLT/SPHERE. Astron. Astrophys. 595, A112 (2016).

85. Goodman, J. & Rafikov, R. R. Planetary torques as the viscosity of protoplanetary disks. Astrophys. J. 552, 793-802 (2001).

86. Rafikov, R. R. Nonlinear propagation of planet-generated tidal waves. Astrophys. J. 569, 997-1008 (2002).

87. Ogilvie, G. I. & Lubow, S. H. On the wake generated by a planet in a disc. Mon. Not. R. Astron. Soc. 330, 950-954 (2002).

88. Bollati, F., Lodato, G., Price, D. J. & Pinte, C. The theory of kinks – I. A semi-analytic model of velocity perturbations due to planet–disc interaction. Mon. Not. R. Astron. Soc. 504, 5444–5454 (2021).

89. Hilder, T., Fasano, D., Bollati, F. & Vandenberg, J. Wakeflow: a Python package for semi-analytic models of planetary wakes. J. Open Source Softw. 8, 4863 (2023).

90. Zhou, Y. et al. UV-optical emission of AB Aur b is consistent with scattered stellar light. Astron. J. 166, 220 (2023).

91. Biddle, L. I., Bowler, B. P., Zhou, Y., Franson, K. & Zhang, Z. Deep Paß imaging of the candidate accreting protoplanet AB Aur b. Astron. J. 167, 172 (2024).

92. Currie, T. Direct imaging detection of the protoplanet AB Aur b at wavelengths covering Paß. Res. Notes AAS 8, 146 (2024).

93. Zhu, Z., Dong, R., Stone, J. M. & Rafikov, R. R. The structure of spiral shocks excited by planetary-mass companions. Astrophys. J. 813, 88 (2015).

94. Zhang, S. & Zhu, Z. The effects of disc self-gravity and radiative cooling on the formation of gaps and spirals by young planets. Mon. Not. R. Astron. Soc. 493, 2287-2305 (2020).

95. Dullemond, C. P. et al. The Disk Substructures at High Angular Resolution Project (DSHARP). VI. Dust trapping in thin-ringed protoplanetary disks. Astrophys. J. Lett. 869, L46 (2018).

96. Birnstiel, T. et al. The Disk Substructures at High Angular Resolution Project (DSHARP). V. Interpreting ALMA maps of protoplanetary disks in terms of a dust model. Astrophys. J. Lett. 869, L45 (2018).

**Acknowledgements** We thank our referees for their careful and insightful comments that improved the manuscript. We thank K. Kratter for enlightening discussions and valuable suggestions. J.S. thanks R. Loomis, S. Wood and T. Ashton at the North American ALMA Science Center (NAASC) for providing science support and technical guidance on the ALMA data as part of a data reduction visit to the NAASC, which was financed by the NAASC. The reduction and imaging of the ALMA data were performed on NAASC computing facilities. J.S. thanks C. Pinte, D. Price and J. Calcino for support with MCFOST, L. Keyte and F. Zagaria for discussions on self-calibrating ALMA data and C. White for sharing perceptually uniform colour maps. J.S. acknowledges financial support from the Natural Sciences and Engineering Research Council of Canada (NSERC) through the Canada Graduate Scholarships Doctoral (CGSD) programme. R.D. acknowledges financial support provided by the NSERC through a Discovery Grant, as well as the Alfred P. Sloan Foundation through a Sloan Research Fellowship. C.L. and G.L. acknowledge funding from the European Union's Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement no. 823823 (RISE DUSTBUSTERS project). C.L. acknowledges funding from the UK Science and Technology Facilities Council (STFC) through the consolidated grant ST/W0009971.1. B.V. acknowledges funding from the ERC CoG project PODCAST no. 864965. Y.-W.T. acknowledges support through National Science and Technology Council grant nos. 111-2112-M-001-064- and 112-2112-M-001-066-. J.-H. was supported by Japan Society for the Promotion of Science (JSPS) KAKENHI grant nos. 21H00059, 22H01274 and 23K03463. This paper makes use of the following ALMA data: ADS/JAO. ALMA#2021.1.00690.S. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), MOST and ASIAA (Taiwan) and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC; https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. Based on data products created from observations collected at the European Organisation for Astronomical Research in the Southern Hemisphere under ESO programme 0104.C-0157(B). This work has made use of the SPHERE Data Centre, jointly operated by OSUG/IPAG (Grenoble), PYTHEAS/LAM/CESAM (Marseille), OCA/Lagrange (Nice), Observatoire de Paris/LESIA (Paris) and Observatoire de Lyon. This research used the Canadian Advanced Network for Astronomical Research (CANFAR) operated in partnership by the Canadian Astronomy Data Centre and the Digital Research Alliance of Canada, with support from the National Research Council of Canada, the Canadian Space Agency, CANARIE and the Canada Foundation for Innovation.

**Author contributions** R.D. led the ALMA proposal. J.S. processed the ALMA data. J.H. processed the VLT/SPHERE data. C.H. performed the SPH simulations. J.S. performed the radiative-transfer calculations. C.L. and G.L. developed the analytic model. J.S. performed all presented analyses. J.S. and R.D. wrote the manuscript. All co-authors provided input to the ALMA proposal and/or the manuscript.

**Competing interests** The authors declare no competing interests.

**Additional information**

**Supplementary information** The online version contains supplementary material available at https://doi.org/10.1038/s41586-024-07877-0.

**Correspondence and requests for materials** should be addressed to Jessica Speedie or Ruobing Dong.

**Peer review information** Nature thanks Jonathan Williams and the other, anonymous, reviewer(s) for their contribution to the peer review of this work. Peer reviewer reports are available.

**Reprints and permissions information** is available at http://www.nature.com/reprints.
<figure><img src="image_7.png" /></figure>

**Extended Data Fig. 1** | Moment maps: AB Aur observations and GI disk simulations. **a–c**, Integrated intensity (moment 0), intensity-weighted mean velocity (moment 1) and intensity-weighted linewidth (moment 2) maps for the ALMA $^{13}$CO observations towards AB Aur. Panel **b** appears in the main text as

Fig. 2a. **d–f**, Moment 0, 1 and 2 maps for the synthetic ALMA $^{13}$CO observations of the SPH GI disk simulation. Like the AB Aur observations, the simulated GI disk shows a prominent GI wiggle along the southern minor axis (indicated by white arrows).
<figure><img src="image_8.png" /><figcaption>**Extended Data Fig. 2** | Filtered moment maps: AB Aur observations and Gl disk simulations. Expanding kernel filter residuals of the maps shown in Extended Data Fig. 1, highlighting global spirals and velocity disturbances</figcaption></figure>

generated by Gl. Panels a–c appear in the main text as Fig. 1b–d. The minor-axis Gl wiggle indicated by arrows in Extended Data Fig. 1b,e is shown here as an isovelocity contour at $v_{\text{los}} = v_{\text{sys}} \pm v_{\text{chan}}$ in all panels.
<figure><img src="image_9.png" /></figure>

**Extended Data Fig. 3** | Global GI wiggles in analytic models, SPH simulations and the AB Aur disk. Isovelocity contours in line-of-sight velocity maps at the velocity values indicated by the colour bar. **a**, $v_{\text{los}}$ map of the 2D analytic GI disk model (shown in Fig. 2b). **b**, $v_{\text{los}}$ map of the 2D analytic Keplerian disk model (shown in Fig. 2b inset). **c**, Synthetic ALMA $^{13}$CO moment 1 map for the 3D SPH GI

disk simulation (shown in Fig. 2c). **d**, Synthetic ALMA $^{13}$CO moment 1 map for the 3D SPH Keplerian disk simulation (shown in Fig. 2c inset). **e**, Observed ALMA $^{13}$CO moment 1 map for the AB Aur disk, imaged with robust 0.5 (shown in Fig. 2a). **f**, Like **e** but imaged with robust 1.5.
<figure><img src="image_10.png" /><figcaption>**Extended Data Fig. 4** | Obtaining velocity residuals in the AB Aur disk.
a, ALMA $^{13}$CO moment 1 map, imaged with robust 0.5, as shown in Fig. 2a.
b, Background model made with a Keplerian rotation profile, assuming a geometrically thin axisymmetric disk (equation (1)). c, Velocity residuals after subtracting the model in panel b. Global spiral substructure is visible but</figcaption></figure>

unevenly so. The model does not capture the non-axisymmetric emission surface morphology and/or super-Keplerian rotation. d, Background model made with the expanding kernel filter (equation (15)). e, Velocity residuals after subtracting the model in panel d, as shown in Fig. 1b. f–j, Like a–e but with the ALMA C$^{18}$O moment 1 map, imaged with robust 1.5.
<figure><img src="image_11.png" /></figure>

**Extended Data Fig. 5** | Kinematics of GI-driven spiral arms. **a–d**, 2D analytic modelling<sup>23</sup>. **e–h**, Synthetic ALMA <sup>13</sup>CO observations of the 3D SPH GI disk simulation. **i–l**, ALMA observations of the AB Aur disk. **a**, Disk surface density (equations (4) and (5)). **b**, Line-of-sight velocity (equation (2)), as in Fig. 2b. **c**, Velocity residuals from Keplerian (that is, subtracting Fig. 2b inset). **d**, Line-of-sight component of the radial velocity (first term of equation (2)). **e**, Filtered

moment 0. **f**, Moment 1. **g**, Moment 1 residuals from Keplerian. **h**, Filtered moment 1. **i**, ALMA <sup>13</sup>CO filtered moment 0. **j**, ALMA <sup>13</sup>CO filtered moment 2. **k**, ALMA <sup>13</sup>CO filtered moment 1. **l**, ALMA <sup>18</sup>O filtered moment 1 (robust 1.5). Panel l inset overlays the VLT/SPHERE H-band scattered-light spirals S1–S7 (refs. 37,38) in red and <sup>13</sup>CO spirals S1–S9 we identify in black.
<figure><img src="image_12.png" /></figure>

**Extended Data Fig. 6** | Methods for isolating the sinusoidal component of the southern minor-axis PV wiggle in the AB Aur disk. a, Detrending the ALMA $^{13}$CO line centres from Fig. 3a with linear and quadratic trend lines found by a least-squares fit. b, Detrending with the expanding kernel high-pass filter, varying the kernel width parameter $w_0$ and keeping the kernel radial power-law index fixed to $\gamma = 0.25$ (equation (15)). We find the background trend lines by

extracting the velocity values from the high-pass-filter background map (for example, Extended Data Fig. 4d) within the same 0.5°-wide wedge-shaped mask as we do for the line centres, positioned along the southern disk minor axis. c, Like b but varying $\gamma$ and keeping $w_0$ fixed to $w_0 = 0.30''$. The high-pass-filter detrending approach converges to the same measured PV wiggle magnitude as the quadratic fit approach.
<figure><img src="image_13.png" /></figure>

**Extended Data Fig. 7** | PV wiggle morphology, magnitude and disk mass recovery in the SPH GI disk simulation. Like Fig. 4 but for the synthetic ALMA observations of the SPH GI disk simulation. **a**, The synthetic ALMA $^{13}$CO line centres along the southern minor axis from Fig. 3c after quadratic detrending. Uncertainties on the line centres are shown by yellow-shaded regions. The magnitude of this PV wiggle is measured to be 39.1 ± 1.9 m s$^{-1}$. The analytic model shown in the background for qualitative comparison has the same parameters as the underlying SPH simulation ($M_{\text{disk}}/M_{\star} = 0.29$ and $\beta = 10$) and its PV wiggle magnitude is 39.0 m s$^{-1}$. **b**, As in Fig. 4c, a map of the minor-axis PV wiggle magnitude of 60 × 60 analytic models on a grid of disk-to-star mass ratios and cooling timescales. A contour is drawn at the measured magnitude of the synthetic $^{13}$CO PV wiggle in panel **a** and dashed lines represent the quoted uncertainties. The technique successfully recovers the disk mass set in the SPH simulation.
<figure><img src="image_14.png" /></figure>

**Extended Data Fig. 8 | Comparisons with further sets of analytic models.**

Like Fig. 4 but varying the azimuthal wavenumber *m* and surface density power-law index *p* in the comparison grid of analytic GI model disks. Each upper subpanel shows the quadratically detrended ¹³CO and C¹⁸O line centres (yellow) behind a demonstrative analytic PV wiggle (black) computed with the combination of *m* and *p* indicated by the row and column labels (keeping *M*<sub>disk</sub>/*M*<sub>*</sub> = 0.3 and *β* = 10 fixed). Each lower subpanel shows the corresponding

map of PV wiggle magnitude computed for a 60 × 60 grid of analytic models in *M*<sub>disk</sub>/*M*<sub>*</sub> and *β*, again with the combination of *m* and *p* indicated by the row and column labels. The two yellow contours are drawn at the magnitude values measured for the observed AB Aur ¹³CO and C¹⁸O southern minor-axis PV wiggles. The white-shaded region between two white curves represents plausible *β* ranges from *r* = 1–5". The combination shown in Fig. 4c is *m* = 3, *p* = 1.0.
<figure><img src="image_15.png" /></figure>

**Extended Data Fig. 9** | **Candidate sites of planet formation.** Coloured crosses mark the locations of candidate protoplanets reported in the literature³⁵,⁴⁵,⁴⁶. A table providing the coordinates of the candidates on the sky, estimated masses and the reporting references is available as source data. a, Filtered ALMA ¹³CO moment 0 map, as in Fig. 1c. b, VLT/SPHERE H-band scattered-light image (ref. 35), as in Fig. 1a. The inset zooms into the central 2" × 2" region to

show the spiral structures in different tracers at spatial scales unresolved by the present ALMA observations. The H-band scattered-light image is shown after high-pass filtering and orange contours show the two spirals identified in ALMA ¹²CO J = 2–1 moment 0 (ref. 45) at levels from 25 to 50 mJy per beam km s⁻¹ in increments of 5 mJy per beam km s⁻¹.
Extended Data Table 1 | Details of the ALMA Band 6 observations

<table><thead><tr><th>UTC Date</th><th>Time on source <sup>a</sup><br>(min)</th><th>N<sub>ant</sub></th><th>Baselines<br>(m)</th><th>pwv<br>(mm)</th><th>Bandpass</th><th>Calibrators<br>Flux</th><th>Phase</th></tr></thead><tbody><tr><td>2022-04-19</td><td>41</td><td>41</td><td>14-500</td><td>0.74</td><td>J0510+1800</td><td>J0510+1800</td><td>J0438+3004</td></tr><tr><td>2022-05-15</td><td>32.4</td><td>45</td><td>15-680</td><td>0.69</td><td>J0510+1800</td><td>J0510+1800</td><td>J0438+3004</td></tr><tr><td>2022-07-17 <sup>b</sup></td><td>47.4</td><td>41</td><td>15-2617</td><td>0.24</td><td>J0510+1800</td><td>J0510+1800</td><td>J0438+3004</td></tr><tr><td>2022-07-19</td><td>47.4</td><td>42</td><td>15-2617</td><td>1.49</td><td>J0510+1800</td><td>J0510+1800</td><td>J0438+3004</td></tr><tr><td>2022-07-20</td><td>34.8</td><td>44</td><td>15-2617</td><td>2.85</td><td>J0510+1800</td><td>J0510+1800</td><td>J0438+3004</td></tr><tr><td>2022-07-21</td><td>47.4</td><td>42</td><td>15-2617</td><td>2.37</td><td>J0510+1800</td><td>J0510+1800</td><td>J0438+3004</td></tr><tr><td>2022-07-22</td><td>47.3</td><td>44</td><td>15-2617</td><td>0.84</td><td>J0510+1800</td><td>J0510+1800</td><td>J0438+3004</td></tr><tr><td>2022-07-22</td><td>47.4</td><td>41</td><td>15-2617</td><td>0.63</td><td>J0510+1800</td><td>J0510+1800</td><td>J0438+3004</td></tr></tbody></table>