The canonical theory for planet formation in circumstellar disks proposes that planets are grown from initially much smaller seedsSUP15SUP. The longconsidered alternative theory proposes that giant protoplanets can be formed directly from collapsing fragments of vast spiral armsSUP611SUP induced by gravitational instabilitySUP1214SUPif the disk is gravitationally unstable. For this to be possible the disk must be massive compared with the central star a disktostar mass ratio of 110 is widely held as the rough threshold for triggering gravitational instability inciting substantial nonKeplerian dynamics and generating prominent spiral armsSUP1518SUP. Although estimating disk masses has historically been challengingSUP1921SUP the motion of the gas can reveal the presence of gravitational instability through its effect on the diskvelocity structureSUP2224SUP. Here we present kinematic evidence of gravitational instability in the disk around AB Aurigae using deep observations of SUP13SUPCO and CSUP18SUPO line emission with the Atacama Large Millimetersubmillimeter 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.
['Astrophysics - Earth and Planetary Astrophysics']
Gravitational instability in a planetforming disk
{'Gravitational instability in a planet-forming disk': "J. Speedie 1 , R. Dong 1 , 2 , C. Hall 3 , 4 , C. Longarini 5 , 6 , B. Veronesi 7 , T. Paneque-Carreño 8 , 9 , G. Lodato 5 , Y. Tang 10 , R. Teague 11 , J. Hashimoto 12 , 13 , 14  \n- 1 Department of Physics & Astronomy, University of Victoria, Victoria, BC, V8P 5C2, Canada\n- 2 Kavli Institute for Astronomy and Astrophysics, Peking University, Beijing 100871, People's Republic of China  \n- Department of Physics and Astronomy, The University of Georgia, Athens, GA 30602, USA  \n- 7 Univ Lyon, Univ Lyon1, Ens de Lyon, CNRS, Centre de Recherche Astrophysique de Lyon UMR5574, F-69230, Saint-Genis,-Laval, France\n- 8\n- Leiden Observatory, Leiden University, P.O. Box 9513, NL-2300 RA Leiden, the Netherlands 9 European Southern Observatory, Karl-Schwarzschild-Str 2, 85748 Garching, Germany\n- 10 Academia Sinica, Institute of Astronomy and Astrophysics, 11F of AS/NTU AstronomyMathematics Building, No.1, Sec. 4, Roosevelt Rd., Taipei, Taiwan  \n- 12 Astrobiology Center, National Institutes of Natural Sciences, 2-21-1 Osawa, Mitaka, Tokyo 1818588, Japan\n- 13 Subaru Telescope, National Astronomical Observatory of Japan, Mitaka, Tokyo 181-8588, Japan 14 Department of Astronomy, School of Science, Graduate University for Advanced Studies (SOKENDAI), Mitaka, Tokyo 181-8588, Japan  \nThe 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 (GI) 12-14 -if the disk is gravitationally unstable. For this to be possible, the disk must be massive compared to the central star: a disk-to-star mass ratio of 1 / 10 is widely held as the rough threshold for triggering GI, inciting significant non-Keplerian dynamics and generating prominent spiral arms 15-18 . While estimating disk masses has historically been challenging 19-21 , the motion of the gas can reveal the presence of GI 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 C 18 O line emission with the Atacama Large Millimeter/submillimeter Array (ALMA). The observed kinematic signals strongly resemble predictions from simulations and analytic modeling. From quantitative comparisons, we infer a disk mass of up to 1 / 3 the stellar mass enclosed within 1 '' to 5 '' on the sky.  \nWe targeted the disk around AB Aurigae (AB Aur), a 2 . 5 -4 . 4 Myr old 25-28 Herbig Ae 29 star of intermediate mass ( M ⋆ = 2 . 4 M ⊙ ) 26,27,30 at a distance of 155 . 9 ± 0 . 9 pc 31 . AB Aur is at a relatively late stage of protostellar evolution, classified as a Class II Young Stellar Object 32,33 (YSO). To probe the velocity structure of the disk, we obtained deep ALMA Band 6 observations of molecular emission lines 13 CO ( J = 2 -1 ) and C 18 O ( J = 2 -1 ) with high velocity resolution (channel widths of v chan = 42 m/s and 84 m/s respectively). The observations were taken in two array configurations with baselines ranging 14 to 2 , 216 m, reaching a total on-source integration time of 5.75 hours. 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 line width (moment 2). This collection is shown in Extended Data Figure 1.  \nTo reveal the spiral arms in the disk, we apply a high-pass filter 34 (see Methods) to the ALMA 13 CO moment maps (Figure 1bcd). In the filtered line-of-sight velocity (moment 1) map, we observe spiral-shaped disturbances in the gas velocity field throughout the disk (Figure 1b). With the filtered velocity-integrated intensity (moment 0) and line width (moment 2) maps, we visually highlight regions of peak density and temperature (Figure 1cd). Compression and shock-heating are expected to lead to temperature enhancements (and thus localized line broadening) within GI-induced density spirals in self-regulating disks 13,23 . The VLT/SPHERE H -band scattered light image of AB Aur originally presented in Boccaletti et al. (2020) 35 is shown for comparison (Figure 1a). Scattered light comes from the disk surface, probing the distribution of (sub-)micron-sized dust usually well-coupled with the gas. Previous simulations have shown that GIinduced 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 , though not all occupy the same radial region and some may be branches of adjacent arms 38 . The disk rotates counter-clockwise (the spiral arms are trailing), and the south side is the near side, tilted toward us 38-40 .  \nTo provide a qualitative comparison to the ALMA observations, we run 3D smoothed-particle hydrodynamic (SPH) simulations of a gravitationally unstable disk (see Methods). The simulations were postprocessed 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 ⋆ ≳ 0 . 1 ), we set the total gas mass to 0 . 3 × the mass of the star. For sustained spiral arms, we set the cooling timescale to 10 × 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 Figures 1 and 2). Overall, the AB Aur disk hosts a global architecture of spiral arms at 100 to 1,000 au scales across all azimuths in multi-wavelength observations tracing different disk components and quantities, strongly indicating ongoing gravitational instability.  \nOne characteristic kinematic feature in the AB Aur disk can be found in the isovelocity curve at the systemic velocity v sys in the moment 1 map - Figure 2a shows a sinusoidal pattern at v los = v sys (along the minor axis; white color), more prominent towards the south. This signature, known as a 'minor axis GI wiggle' 22 , has been predicted in hydrodynamic simulations 22,24 and analytic theory 23 as a clear  \nkinematic signature of gravitational instability (Figure 2bc). It is one of a global set of GI wiggles in isovelocity curves we observe throughout the AB Aur disk (Extended Data Figure 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 (Figure 2c), completely distinct from the linear pattern found in a disk undergoing Keplerian rotation with no radial motions (Figure 2bc insets).  \nAmong all GI wiggles, the minor axis GI wiggle has been known and targeted in past studies for its convenience in quantitative analysis 23,24 . Due to projection effects, only the radial and vertical components of the disk velocity field ( v r or v z ) contribute to v 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 22 . 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 sys where the spiral crosses the minor axis (c.f. Extended Data Figure 5). The filtered moment 1 map in Figure 1b displays redand blue-shift patterns corresponding to convergent flows toward spiral S5 (visible in both scattered light and 13 COmoment 0 and 2; Figure 1acd), supporting the interpretation that the GI wiggle along the southern minor axis in Figure 2a is generated by a self-gravitating spiral arm.  \nHaving identified evidence of gravitational instability in disk kinematics and in the detections of spirals across multiple tracers and moment maps, we now quantitatively analyze the GI wiggle along the southern minor axis to constrain the disk mass. We extract the 13 CO and C 18 O emission spectra along the southern disk minor axis (Figure 3ab) 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 Figure 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 centers 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 center and yields statistically meaningful and robust uncertainties 42 . We find remarkably similar sinusoidal morphology between the PV wiggles in 13 CO and C 18 O emission (Figure 4a).  \nTheoretical 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 disk /M ⋆ ) 2 and β -1 / 2 (Eqns. 11 & 18 in 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 Longarini et al. (2021) 23 , we employ a statistical metric to quantify the 'magnitude' of the minor axis PV wiggle, defined as the standard deviation of the line center 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 ± 2 . 9 m/s for the southern minor axis PV wiggle in 13 CO and 44 . 2 ± 1 . 3 m/s in C 18 O (Figure 4b). For comparison, the gravitationally unstable disk in the SPH simulation has a southern minor axis PV wiggle in 13 CO emission with quantitatively similar amplitude and sinusoidal morphology (Figure 3c), and a magnitude of 39 . 1 ± 1 . 8 m/s (Extended Data Figure 7a).  \nQuantifying the minor axis PV wiggle magnitude as above, we perform comparisons against analytic models to identify the combinations of disk mass ( M disk /M ⋆ ) and cooling timescale ( β ) that satisfy the AB Aur observations. A proof of concept of this technique with the SPH simulation is shown in Extended Data Figure 7b. Using the analytic modeling code giggle 1 of Longarini et al. (2021) 23 (Methods), we calculate the minor axis PV wiggle magnitude in gravitationally unstable disk models for 60 × 60 combinations of M disk /M ⋆ and β , letting each vary within the ranges 0 . 0 ≤ M disk /M ⋆ ≤ 0 . 4 and 10 -2 ≤ β ≤ 10 2 . A demonstrative analytic curve for the minor axis PV wiggle from the same model shown in Figure 2b is underlaid in Figure 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 ≲ M disk /M ⋆ ≲ 0 . 3 for a cooling timescale of 0 . 1 < β < 10 . This result is robust to plausible variations in the analytic model parameter choices (Extended Data Figure Figure 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, Figure 4c also displays ranges of β derived from independent radiative cooling prescriptions (see Methods).  \nThe detection of GI in the disk around AB Aur, a Class II YSO 32,33 , demonstrates that gravitational instability can take place during later evolutionary stages. This result, together with previous reports of multiple protoplanet candidates in and amongst spiral arms in the system 35,44-46 (Extended Data Figure 9), provides a direct observational connection between gravitational instability and planet formation. Looking forward, the AB Aur system can be an ideal testbed for understanding how planet formation is facilitated by GI-induced spiral arms - whether by fragmentation into gas clumps enabled by rapid cooling 7-10 ( β ≲ 3 ), or by dust collapse of solids concentrated within spiral arms sustained by slow cooling 47-50 ( β ≳ 5 ).  \n2  \nFigure 1: Global spirals in the AB Aur disk. (a) VLT/SPHERE H -band scattered light image of the AB Aur disk (ref. 35 ) tracing spiral structure in (sub-)micron-sized dust grains. The labelled spirals S1-S7 are taken from previous works (ref. 37,38 ). (b) Filtered ALMA 13 CO intensity-weighted mean velocity (moment 1) map, revealing residual gas motion within the bulk flow. The synthesised 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 13 CO integrated intensity (moment 0) map, highlighting peaks in the gas density and/or temperature. (d) Filtered ALMA 13 CO emission line width (moment 2) map, showing localized line broadening within the spiral arms. Insets in panels (c,d) zoom into the same region as panel (b) inset, showing enhanced gas density/temperature caused by the radially converging flows around S5. 5  \n<!-- image -->  \nFigure 2: Detection of the GI wiggle in the AB Aur disk. (a) ALMA 13 CO intensity-weighted mean velocity (moment 1) map showing line-of-sight velocity ( v los ) of gas in the AB Aur disk. The observations display the 'GI wiggle' along the minor axis (arrow) predicted by Hall et al. (2020) 22 as a clear kinematic signature of gravitational instability. (b) v los map of a gravitationally unstable disk at the inclination and position angle of the AB Aur disk, computed with 2D analytic modeling 23 . Self-gravitating spiral arms crossing the minor axis induce radial motion that appears as a wiggle (arrow). (c) Synthetic ALMA 13 CO moment 1 map of the 3D SPH GI disk simulation, revealing the same GI wiggle signature (arrow). The insets in panels (b) and (c) show corresponding images for Keplerian disks with no radial gas motion, where the isovelocity curve at the systematic velocity appears as a straight line along the minor axis.  \n<!-- image -->  \nFigure 3: The GI wiggle in position-velocity space, or 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 centers are shown as yellow points. The insets at the bottom right of each panel show the corresponding line center map, with black 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 centers are extracted. In all PV diagram panels, the grey box in the bottom left corner has horizontal width equal to the beam major axis and vertical height equal to the channel width. (a) ALMA observations of the AB Aur disk in 13 CO and (b) in C 18 O. (c,d) Synthetic ALMA 13 COobservations 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 β = 10 . In (d) , the simulated disk has its velocity structure artificially post-processed to be Keplerian.  \n<!-- image -->  \nFigure 4: PVwiggle mophology, magnitude, and constraints on the AB Aur disk mass. (a) The ALMA 13 CO and C 18 Oline centers along the southern minor axis from Figure 3ab, after quadratic detrending (Methods). Uncertainties on the line centers are shown by yellow shaded regions. For qualitative comparison, the minor axis PV wiggle of the analytic GI model disk from Figure 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 in 13 CO and 44 . 2 ± 1 . 3 m/s in C 18 O. (c) A map of the minor axis PV wiggle magnitude of 3600 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 σ , 2 σ and 3 σ departures from a Keplerian signal in 13 CO (c.f. Fig. 3d). White horizontal bars indicate independently derived β ranges at a selection of M disk /M ⋆ values (Methods).  \n<!-- image -->", 'References': "- 1. Salyk, C. et al. Measuring Protoplanetary Disk Accretion with H I Pfund β . Astrophys. J. 769 , 21 (2013).\n- 2. 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).\n- 3. van den Ancker, M. E. et al. HIPPARCOS data on Herbig Ae/Be stars: an evolutionary scenario. Astron. Astrophys. 324 , L33-L36 (1997).\n- 4. 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).\n- 5. 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).\n- 6. Garufi, A. et al. The SPHERE view of the Taurus star-forming region. arXiv e-prints arXiv:2403.02158 (2024).\n- 7. Hartmann, L., Calvet, N., Gullbring, E. & D'Alessio, P. Accretion and the Evolution of T Tauri Disks. Astrophys. J. 495 , 385-400 (1998).\n- 8. Dong, R., Najita, J. R. & Brittain, S. Spiral Arms in Disks: Planets or Gravitational Instability? Astrophys. J. 862 , 103 (2018).\n- 9. 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).\n- 10. 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).\n- 11. Nakajima, T. & Golimowski, D. A. Coronagraphic Imaging of Pre-Main-Sequence Stars: Remnant Envelopes of Star Formation Seen in Reflection. Astron. J. 109 , 1181 (1995).\n- 12. 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).\n- 13. 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).\n- 14. Fukagawa, M. et al. Spiral Structure in the Circumstellar Disk around AB Aurigae. Astrophys. J. Lett. 605 , L53-L56 (2004).\n- 15. Hall, C. et al. The Temporal Requirements of Directly Observing Self-gravitating Spiral Waves in Protoplanetary Disks with ALMA. Astrophys. J. 871 , 228 (2019).  \n- 16. Ediss, G. A. et al. Narayanan, G. (ed.) ALMA Band 6 Cartridge: Design and Performance . (ed.Narayanan, G.) Fifteenth International Symposium on Space Terahertz Technology , 181-188 (2004).\n- 17. Cornwell, T. J. Multiscale CLEAN Deconvolution of Radio Synthesis Images. IEEE Journal of Selected Topics in Signal Processing 2 , 793-801 (2008).\n- 18. Kepley, A. A. et al. Auto-multithresh: A General Purpose Automasking Algorithm. PASP 132 , 024505 (2020).\n- 19. Leroy, A. K. et al. PHANGS-ALMA Data Processing and Pipeline. Astrophys. J. Suppl. Ser. 255 , 19 (2021).\n- 20. Jorsater, S. & van Moorsel, G. A. High Resolution Neutral Hydrogen Observations of the Barred Spiral Galaxy NGC 1365. Astron. J. 110 , 2037 (1995).\n- 21. Czekala, I. 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).\n- 22. Teague, R. & Foreman-Mackey, D. bettermoments: A robust method to measure line centroids. Zenodo (2018).\n- 23. Teague, R. & Foreman-Mackey, D. A Robust Method to Measure Centroids of Spectral Lines. Research Notes of the American Astronomical Society 2 , 173 (2018).\n- 24. Teague, R. eddy: Extracting Protoplanetary Disk Dynamics with Python. The Journal of Open Source Software 4 , 1220 (2019).\n- 25. 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).\n- 26. Gaia Collaboration et al. The Gaia mission. Astron. Astrophys. 595 , A1 (2016).\n- 27. Gaia Collaboration et al. Gaia Data Release 3. Summary of the content and survey properties. Astron. Astrophys. 674 , A1 (2023).\n- 28. 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).\n- 29. Price, D. J. et al. Phantom: A Smoothed Particle Hydrodynamics and Magnetohydrodynamics Code for Astrophysics. Publ. Astron. Soc. Aust. 35 , e031 (2018).\n- 30. Bate, M. R., Bonnell, I. A. & Price, N. M. Modelling accretion in protobinary systems. Mon. Not. R. Astron. Soc. 277 , 362-376 (1995).\n- 31. Cullen, L. & Dehnen, W. Inviscid smoothed particle hydrodynamics. Mon. Not. R. Astron. Soc. 408 , 669-683 (2010).  \n- 32. Gammie, C. F. Nonlinear Outcome of Gravitational Instability in Cooling, Gaseous Disks. Astrophys. J. 553 , 174-183 (2001).\n- 33. Hall, C. et al. Predicting the Kinematic Evidence of Gravitational Instability. Astrophys. J. 904 , 148 (2020).\n- 34. Terry, J. P. et al. Constraining protoplanetary disc mass using the GI wiggle. Mon. Not. R. Astron. Soc. 510 , 1671-1679 (2022).\n- 35. Longarini, C. et al. Investigating Protoplanetary Disk Cooling through Kinematics: Analytical GI Wiggle. Astrophys. J. Lett. 920 , L41 (2021).\n- 36. 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).\n- 37. Pinte, C., Ménard, F., Duchêne, G. & Bastien, P. Monte Carlo radiative transfer in protoplanetary disks. Astron. Astrophys. 459 , 797-804 (2006).\n- 38. Pinte, C. et al. Benchmark problems for continuum radiative transfer. High optical depths, anisotropic scattering, and polarisation. Astron. Astrophys. 498 , 967-980 (2009).\n- 39. Pinte, C. et al. Kinematic Evidence for an Embedded Protoplanet in a Circumstellar Disk. Astrophys. J. Lett. 860 , L13 (2018).\n- 40. 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).\n- 41. Hillenbrand, L. A., Strom, S. E., Vrba, F. J. & Keene, J. Herbig Ae/Be Stars: Intermediate-Mass Stars Surrounded by Massive Circumstellar Accretion Disks. Astrophys. J. 397 , 613-643 (1992).\n- 42. Natta, A. et al. A reconsideration of disk properties in Herbig Ae stars. Astron. Astrophys. 371 , 186197 (2001).\n- 43. Currie, T. et al. Images of embedded Jovian planet formation at a wide separation around AB Aurigae. Nature Astronomy 6 , 751-759 (2022).\n- 44. Lodato, G. Classical disc physics. New Astronomy Reviews 52 , 21-41 (2008).\n- 45. Cossins, P., Lodato, G. & Clarke, C. J. Characterizing the gravitational instability in cooling accretion discs. Mon. Not. R. Astron. Soc. 393 , 1157-1173 (2009).\n- 46. 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).\n- 47. 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).  \n- 48. 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).\n- 49. 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).\n- 50. Pérez, L. M. et al. Spiral density waves in a young protoplanetary disk. Science 353 , 1519-1521 (2016).\n- 51. Meru, F. et al. On the Origin of the Spiral Morphology in the Elias 2-27 Circumstellar Disk. Astrophys. J. Lett. 839 , L24 (2017).\n- 52. 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).\n- 53. 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).\n- 54. 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).\n- 55. Goodman, J. & Rafikov, R. R. Planetary Torques as the Viscosity of Protoplanetary Disks. Astrophys. J. 552 , 793-802 (2001).\n- 56. Rafikov, R. R. Nonlinear Propagation of Planet-generated Tidal Waves. Astrophys. J. 569 , 997-1008 (2002).\n- 57. 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).\n- 58. 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).\n- 59. Hilder, T., Fasano, D., Bollati, F. & Vandenberg, J. Wakeflow: A Python package for semi-analytic models of planetary wakes. The Journal of Open Source Software 8 , 4863 (2023).\n- 60. Zhou, Y. et al. UV-Optical Emission of AB Aur b is Consistent with Scattered Stellar Light. arXiv e-prints arXiv:2308.16223 (2023).\n- 61. 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).\n- 62. Currie, T. Direct Imaging Detection of the Protoplanet AB Aur b at Wavelengths Covering Pa β . Research Notes of the American Astronomical Society 8 , 146 (2024).\n- 63. Zhu, Z., Dong, R., Stone, J. M. & Rafikov, R. R. The Structure of Spiral Shocks Excited by Planetarymass Companions. Astrophys. J. 813 , 88 (2015).  \n- 64. 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).\n- 65. 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).\n- 66. 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).\n- 67. 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).", 'Methods': "Additional information on the source. AB Aur is accreting from the disk at a rate ˙ M ∼ 10 -7 M ⊙ yr -1 (ref. 1 ), within the range expected for modest GI-driven accretion ( 10 -7 -10 -6 M ⊙ yr -1 ) 2 . This accretion rate, taken together with the current age t 0 = 2 . 5 -4 . 4 Myr 3-6 , implies a high 'latent disk mass': M latent disk = ˙ M ( t 0 ) × t 0 = 0 . 25 -0 . 44 M ⊙ , or M latent disk /M star ∼ 0 . 1 -0 . 2 . M latent disk provides an accretion rate-based assessment of disk mass, assuming a constant stellar accretion rate ˙ M and we are observing the system mid-way through the disk's lifetime 7, 8 . This is a conservative estimate as the accretion rate at earlier epochs is likely higher 9 . In millimeter continuum observations, the disk shows a dust ring at ∼ 1 '' and a cavity inside 10 , likely caused by the trapping of millimeter-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 10-13 is likely encouraging GI by providing a source of mass to maintain a high M disk /M ⋆ value 14,15 .  \nALMAobservations. We observed AB Aur with ALMA in April, May and September 2022 under ALMA program ID 2021.1.00690.S (PI: R. Dong). Measurements were taken with the Band 6 receivers 16 in array configurations C-3 (2 execution blocks) and C-6 (6 execution blocks). In total, the 8 execution blocks reached an on-source integration time of 5.75 hours, making this the longest fine-kinematics ( v chan < 100 m/s) program toward a single protoplanetary disk to date. Extended Data Table 1 provides details of the observations. We centered one spectral window (SPW) at the 13 CO J = 2 -1 molecular emission line transition rest frequency (220.3986 GHz), covering a bandwidth of 58.594 MHz with 1920 channels, resulting in the highest achievable spectral resolution of 41.510 m/s after default spectral averaging with N = 2 by Hanning smoothing within the correlator data processor. A second SPW was centered at the C 18 O J = 2 -1 rest frequency (219.5603 GHz) covering the same bandwidth with half as many channels (960 channels; due to sharing a baseband with another SPW), achieving a 83.336 m/s spectral resolution. To enable self-calibration, our correlator setup sampled the continuum in another SPW centered 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 center in the uv -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 center realignments and calibration gain tables (that we generated with the continuum data) to the line data. We performed continuum subtraction in the uv -plane using the uvcontsub task.  \nAll imaging was performed with the CASA tclean task. We used the multiscale deconvolution algorithm 17 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 18 , we adopted an imaging strategy similar to PHANGS-ALMA 19 , in which we clean conservatively, with a broad mask ( usemask='pb'  \nand pbmask=0.2 ), forcing frequent major cycles 2 . 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.02 times the peak residual. We adopted 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 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 FOV out to the primary beam 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 for 13 CO and 84 m/s for C 18 O respectively (nearly native channel spacing). The CLEAN threshold was set to 5 × the rms noise measured in 20 line-free channels of the dirty image cube. We applied JvM correction 20,21 and primary beam correction. The rms noise in the resulting 13 CO cubes imaged with robust 0 . 5 and robust 1 . 5 is 2 . 0 mJy/beam and 1 . 2 mJy/beam respectively, and 0 . 6 mJy/beam in the C 18 Ocube imaged with robust 1 . 5 .  \nWe used the robust 0 . 5 image cubes for our position-position analysis (moment maps; Figures 1 & 2) and the robust 1 . 5 cubes for our position-velocity analysis (PV diagrams and line centers; Figures 3 & 4). We made the moment 0, 1 and 2 maps using the bettermoments 22,23 methods 'collapse\\_zeroth, 'collapse\\_first', and 'collapse\\_percentiles', respectively. We note that we calculate our 'moment 2' maps as the average of the red- and blue-shifted line widths about the intensity-weighted median line center (i.e., as the average of the wpdVr and wpdVb maps). Mathematically this is a different approach to find the line width than the classic moment 2 approach, though in our case we find the two yield nearly identical outcomes. We applied sigma-clipping at 5 × the rms noise and performed no spectral smoothing.  \nGeometric properties. We used the Python package eddy 24 to infer geometric properties of the disk, namely to constrain the disk center x 0 , y 0 , the disk inclination i , the position angle PA , the systemic velocity v sys , and the dynamical stellar mass M ⋆ . We performed an MCMC to fit the C 18 Omoment 1 map (Extended Data Figure 4) with a geometrically thin Keplerian disk rotation profile:  \nv 0 = √ GM ⋆ r · sin i · cos ϕ + v sys , (1)  \nwhere r is the disk radius, ϕ 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 · (ref. 10,25 ), and the distance to 155.9 pc (Gaia DR3 26,27 ). We assumed flat priors for all values and spatially downsampled the rotation map to the beam FWHM prior to 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 additional 10,000 steps. The resulting posterior distributions were x 0 = -5 ± 7 mas, y 0 = -17 ± 7 mas, PA = 236 . 7 ± 0 . 3 · , M ⋆ = 2 . 23 ± 0 . 02 M ⊙ , and v sys = 5858 ± 5 m/s, where we report the uncertainties represented by the 16th and 84th percentiles about the median value. The latter three values are consistent with constraints from previous observations 4, 5, 10, 25, 28 .  \nHydrodynamic simulations and synthetic ALMA observations. We performed 3D global smoothedparticle hydrodynamic (SPH) simulations with the PHANTOM code 29 using 1 million SPH particles. We assumed a central star mass of 2 . 4 M ⊙ (ref. 4, 5 ), represented by a sink particle 30 with accretion radius set to 60 au. The initial inner and outer disk radii were set to r in , SPH = 80 au and r out , SPH = 500 au, respectively. We set the initial gas mass to 0 . 7 M ⊙ , corresponding to M disk /M ⋆ = 0 . 29 . The surface density profile follows Σ ∝ r -p (where the power-law index p = 1 . 0 ), and the sound speed profile follows c s ∝ r -q (where q = 0 . 25 ). The initial disk aspect ratio was set to H/r = 0 . 05 at 80 au. We set α SPH such that α min ≤ α SPH ≤ α max , with α min = 0 . 001 and α max = 1 . 0 , with the value of α SPH set by the Cullen & Dehnen (2010) 31 switch that increases viscosity only in the case of converging flows. This results in a Shakura-Sunyaev viscosity of α SS ≈ 0 . 01 throughout the disk.  \nWe assumed an adiabatic equation of state, with heating from compressional P d V work and shock heating. The disk cools by Gammie cooling 32 (a.k.a. β -cooling) where the cooling timescale is proportional to the local dynamical time by the factor β , such that t cool ( r ) = β Ω -1 ( r ) , where Ω( r ) = ( GM ⋆ /r 3 ) 1 / 2 is the Keplerian frequency. We set β = 10 , a typical value used or found in simulations 33-36 . 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 in , SPH to 1.1 r out , SPH .  \nWe computed the disk thermal structure and 13 CO ( J = 2 -1 ) model line cubes using the Monte Carlo radiative transfer code MCFOST 37,38 . We assumed the 13 CO molecule is in local thermodynamic equilibrium (LTE) with its surroundings and that the dust is in thermal equilibrium with the gas ( T gas = T dust ). We set the 13 CO/H 2 abundance to 7 × 10 -7 (ref. 33,39 ) and we used ≈ 10 7 photon packets to calculate T dust . Voronoi tesselation was performed 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 ⊙ radiating isotropically at an effective temperature T eff = 9770 K, set to match AB Aur 40-43 . The disk was given an inclination of 23 . 2 · , a position angle of 236 . 7 · (where PA is measured east of north to the red-shifted major axis), and placed at a distance of 155 . 9 pc, all consistent with the AB Aur system.  \nWe used the same PHANTOM simulation to create both the GI and Keplerian model line cubes shown in Figure 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 13 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 2048 × 2048 pixels of size 0 . 02 '' . We assumed a turbulent velocity of 0 . 05 km/s.  \nWe generated synthetic ALMA image cubes from the 13 CO model line cubes using syndisk 3 to match the properties of the observed AB Aur 13 CO image cubes (robust 0 . 5 and 1 . 5 ). In the latter case  \nthe model line cube was convolved with a beam of size 0 . 390 '' × 0 . 274 '' and PA -1 . 4 · . Correlated noise was added with an rms of 1 . 2 mJy/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 uv -plane. The synthetic cubes were collapsed into moment maps following the same procedure as the AB Aur data (Extended Data Figure 1).  \nAnalytic modeling. We analytically compute the velocity fields of gravitationally unstable disks using the giggle 4 package developed by Longarini et al. (2021) 35 . Working in 2D polar coordinates ( r , ϕ ), giggle considers a geometrically thin disk with surface density profile Σ 0 ∝ r -p and inclination i , centered on a star of mass M ⋆ . It computes the projected line-of-sight velocity field as  \nv los = ( v r sin ϕ + v ϕ cos ϕ ) sin i + v sys , (2)  \nwhere v r and v ϕ are the radial and azimuthal components of the disk velocity field. The basic state of the disk (i.e., considering only the gravitational potential contribution from the central star) is assumed to be Keplerian: v r = 0 and v ϕ = v Kep . The scheme of the model is to determine the perturbations in v r and v ϕ generated by gravitational instability by taking into account the additional 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 disc 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 Σ spir / Σ 0 saturated to a finite value proportional to the cooling timescale β : 44,45  \nA Σ spir Σ 0 = χβ -1 / 2 , (3)  \nwhere the proportionality factor χ is of order unity 35,45 . The imprinted spiral density perturbation is assumed to be small relative to the background surface density, so that all the relevant quantities (density Σ , gravitational potential Φ , velocities v r and v ϕ , and enthalpy h ) can be written as a linear sum of the basic state and the perturbation:  \nX ( r, ϕ ) = X 0 ( r ) + X spir ( r, ϕ ) . (4)  \nThe spiral perturbation in density is given the form  \nΣ spir ( r, ϕ ) = ℜ [ A Σ spir e j ( mϕ + ψ ( r )) ] , (5)  \nwhere j = √ -1 (as we are using i to represent the disk inclination), and m is the azimuthal wavenumber. The 'shape function' ψ ( r ) is described by m and the spiral pitch angle α pitch as:  \nψ ( r ) = m tan α pitch log r , (6)  \nwhich is related to the radial wavenumber k by dψ/dr = k . The spiral density perturbation necessarily introduces a corresponding perturbation to the gravitational potential:  \nΦ spir ( r, ϕ ) = -2 πG | k | Σ spir ( r, ϕ ) . (7)  \nThe negative proportionality Φ spir ∝ -Σ spir is the definition of self-gravitating spiral arms. As a result, corresponding perturbations in the azimuthal and radial velocities are driven:  \nv r ( r, ϕ ) = ℜ [ A v r ( r ) · e j ( mϕ + ψ ( r )) ] , (8)  \nv ϕ ( r, ϕ ) = ℜ [ A v ϕ ( r ) · e j ( mϕ + ψ ( r )) ] + r Ω , (9)  \nwhere we note r Ω = v Kep because the angular frequency Ω includes super-Keplerian rotation from the disk mass contribution:  \nΩ 2 = GM ⋆ r 3 + 1 r ∂ Φ disk ∂r . (10)  \nBy assuming 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 ϕ ( r ) are determined: 35  \nA v r ( r ) = 2 jmχβ -1 / 2 ( M disk ( r ) M ⋆ ) 2 v Kep ( r ) , (11)  \nA v ϕ ( r ) = -1 2 jχβ -1 / 2 ( M disk ( r ) M ⋆ ) v Kep ( r ) , (12)  \nwhere M disk ( r ) is the disk mass enclosed within radius r . With a surface density profile Σ 0 ( r ) ∝ r -p , then M disk ( r ) ∝ r -p +2 , and the amplitude of the radial perturbation is described by A v r ( r ) ∝ r -2 p +7 / 2 . For p < 7 / 4 , A v r ( r ) is an increasing function of radius. The factor of imaginary number j in Eqn. 11 has important physical consequences: when the real component of A v r ( r ) is taken (Eqn. 8), the radial velocity perturbation is π/ 2 out of phase with the spiral density perturbation (Eqn. 5), and convergent at the locations where Σ spir takes a maximum. Explicitly,  \nv r ( r, ϕ ) ∣ ∣ ∣ ∣ ϕ = π/ 2 ∝ -sin ( m π 2 + ψ ( r ) ) , (13)  \nΣ spir ( r, ϕ ) ∣ ∣ ∣ ∣ ϕ = π/ 2 ∝ cos ( m π 2 + ψ ( r ) ) . (14)  \nFor qualitative visual comparison with the AB Aur moment 1 map in Figure 2a, we compute the projected line-of-sight velocity field of a gravitationally unstable disk with β = 10 and M disk /M ⋆ = 0 . 3 in Figure 2b. We set m = 3 and α pitch = 15 · to approximately match the 13 CO spirals in the AB Aur disk (Figure 1), and assume p = 1 . 0 and χ = 1 . 0 (ref. 45 ). The dominant azimuthal wavenumber is expected to be inversely related to the disk-to-star mass ratio q , roughly obeying m ∼ 1 /q (ref. 45-47 ), so our choice of m = 3 is consistent with M disk /M ⋆ ≈ 0 . 3 .  \n̸  \nRevealing global spiral structure. Weobtain the residual moment maps shown in Figure 1 using a variation on the conventional high-pass filtering (a.k.a. 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 (e.g. 6, 48-53 ). Here, we perform the convolution with a radially expanding kernel 5 - that is, with a Gaussian kernel whose FWHM, w , increases with radial distance from the image center (i.e., with disk radius) with a simple power-law dependence:  \nw ( r ) = w 0 · ( r/r 0 ) γ , (15)  \nwhere 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 adopt w 0 = 0 . 3 '' and γ = 0 . 25 , though we emphasize 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 Figure 4 compares the residual moment 1 maps in 13 CO and C 18 O obtained after subtracting the axisymmetric geometrically thin Keplerian model (Eqn. 1) vs. 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 blue-shifted (red-shifted), hinting at super-Keplerian rotation, signatures of disk mass contributing to the total mass of the system. While 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 displays excess blueshifted velocity in the southeast quadrant of the disk such that the contour of v los = v 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 centers to isolate the sinusoidal component of the southern minor axis PV wiggle in Figure 4a (see section 'Measuring the magnitude of AB Aur's minor axis PV wiggle'). Filtered moment maps for the synthetic ALMA observations of the simulated SPH GI disk are shown in Extended Data Figure 2.  \nGlobal kinematics of self-gravitating spiral arms. Radially convergent motion (as in Figure 1bcd insets) serves as a kinematic signature for the location of self-gravitating spiral arms at disk azimuths where 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 Figure 5c and g provide maps of velocity residuals from Keplerian for the 2D analytic GI disk model and the SPH GI disk simulation. The convergent motion toward the spiral spines is visible for a range of azimuths around the minor axis, but becomes progressively less clear moving toward 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  \nresembles the isolated radial component (panel d). Extended Data Figure 5i-l overlays the locations of 13 CO spirals in the AB Aur disk (from filtered moment 0/2; Figure 1cd) onto the filtered moment 1 maps, in order 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 ϕ (first term of Eqn. 2), and when two spirals are not well separated and their motions superimpose. Three of the seven spiral structures in VLT/SPHERE scattered light appear to be spatially associable with those in 13 CO (S1, S5, S7; panel l inset). Offsets in the southeast quadrant of the disk (S2, S3, S4) may be further indication of a disk warp (Extended Data Figure 4di), or other non-trivial phenomena (e.g., vertical density and temperature gradients, projection effects 54 ).  \nThe kinematic signatures observed in the present ALMA dataset -probing disk scales ∼ 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 55-57 , 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 (Eqns. 11 & 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 58,59 . Overall, the essential characteristic of GI-induced spirals is that they occur globally 33,35 (c.f. Figure 1, Extended Data Figures 2, 3 & 5). In previous datasets probing smaller spatial scales -within the AB Aur disk's central cavity- planetary candidates P1/f1 (ref. 25,48 ), P2/b (ref. 25,43,60-62 ), and f2 (ref. 48 ) are known to be associated with -or driving- spiral arms, as observed in VLT/SPHERE scattered light and/or ALMA 12 CO emission. As shown in Extended Data Figure 9, due to their small separations ( ≲ 0 . 7 '' ), kinematic signatures from these candidates are inaccessible to our ALMA observations. Clump-like signals 'c' and 'd' seen by HST/STIS (ref. 43 ) at wide separations ( ∼ 2 . 75 '' and ∼ 3 . 72 '' respectively) are in locations tentatively suggestive of constituting spiral arm fragments and may warrant further investigation.  \nPosition-velocity analysis. Weuse the robust 1 . 5 image cubes for our position velocity analysis to maximize the recovery of emission at large disk radii. Owing to the clear association with a self-gravitating spiral arm (Figure 1bcd 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 CO and C 18 O emission ( ∼ 3 '' ; c.f. Extended Data Figure 5kl). We obtain the position-velocity diagrams shown in Figure 3 using eddy 24 to extract spectra from pixels within a 0 . 5 · -wide wedge-shaped mask oriented 90 · clockwise of the red-shifted major axis (shown in Figure 3 insets). Our quantitative analysis of the minor axis PV wiggles is performed with maps of the line centers made using the quadratic method of bettermoments 22,23 , which fits a quadratic curve to the spectrum in each pixel of the cube:  \nI ( v ) = a 0 + a 1 ( v -v peak ) + a 2 ( v -v peak ) 2 , (16)  \nwhere v 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 center, σ v los (ref. 23 ). The statistical uncertainty on  \neach line center is computed as:  \nσ v los = √ σ 2 I 8 ( 3 a 2 2 + a 1 2 a 2 4 ) , (17)  \nwhere σ I is the rms noise of the intensities (see ref. 23 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 23 . Prior to the quadratic fitting we spectrally smooth the data with a Savitzky-Golay filter of polynomial order 1 and filter window length of 10 channels ( 420 m/s) in the case of 13 CO and 3 channels ( 252 m/s) in the case of C 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 center and line uncertainty maps within the same wedge mask described above. The extracted line center values are shown as yellow points in Figure 3 and the uncertainties are shown as yellow shaded regions in Figure 4a.  \nMeasuring the magnitude of the minor axis PV wiggle. Following Longarini et al. (2021) 35 , we measure the 'magnitude' of a minor axis PV wiggle as the standard deviation of the line center values over a radial range. Bounded by the inner central cavity and the outer edge of C 18 O emission, we adopt a radial range of 1 . 0 '' to 5 . 0 '' . We estimate the uncertainty on the magnitude measurement using a resampling procedure: we take 10,000 draws from Gaussian distributions centered on the observed line centers with standard deviation σ v los (Eqn. 17) to create 10,000 instances of the minor axis PV wiggle; we compute their magnitudes; and then report the uncertainty as the standard deviation of those 10,000 magnitude estimates.  \nIn addition to the wiggle, the 13 CO and C 18 O emission on the southern disk minor axis also exhibit an underlying monotonic blueward trend with disk radius, seen in Figure 3ab as a subtle downward bend with radius of the line centers, or equivalently in Figure 2a as a westward or clockwise shift in the contour of v los = v sys . We earmark this feature as a possible disk warp (Extended Data Figure 4di), and adopt a least-squares fitting approach to isolate the sinusoidal component of the PV wiggle. This approach yields the background trendline 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 trendline (Extended Data Figure 6a) as it more closely resembles the high-pass filter background curve than a linear one (Extended Data Figure 6bc). We show the quadratically-detrended PV wiggles in Figure 4a and report their magnitudes in Figure 4b. We find very similar magnitudes for both the 13 CO and C 18 O wiggles, despite C 18 O likely tracing lower optical depths in the AB Aur disk. This empirically substantiates comparisons with the 2D analytic model (next section).  \nPerforming 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 (Extended Data Figure 7).  \nConstraining disk mass with quantitative comparisons to analytic models. We perform quantitative comparisons between the observed 13 CO and C 18 Ominor axis PV wiggles and the projected radial velocity component in our analytic model, v r sin i (ref. 35 ). From Eqns. 8 and 11, the projected radial velocity on the  \nminor axis ( ϕ = π/ 2 ) is:  \nv r ( r, ϕ ) ∣ ∣ ∣ ∣ ϕ = π/ 2 · sin i = -2 mχβ -1 / 2 ( M disk ( r ) M ⋆ ) 2 v Kep ( r ) sin ( m π 2 + ψ ( r ) ) · sin i . (18)  \nThis 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 3600 of these curves for a 60 × 60 grid of models with (total enclosed) M disk /M ⋆ linearly spaced ∈ [0 . 0 , 0 . 4] and β logarithmically spaced ∈ [10 -2 , 10 2 ] . Again we set m = 3 and α pitch = 15 · to match the AB Aur disk, and assume p = 1 . 0 and χ = 1 . 0 (ref. 45 ). For qualitative comparison, we plot an example analytic minor axis PV wiggle behind the data in Figure 4a; the model has β = 10 and M disk /M ⋆ = 0 . 3 . We show in Extended Data Figure Figure 8 that m = 3 reproduces the observed wiggles better than other choices, and that p = 1 . 5 could also provide a satisfying match, while p = 2 . 0 is too steep. Since the wiggle amplitude is independent of α pitch (Eqn. 11), the magnitude is constant with α pitch when sampled over the same range in phase (not shown).  \nWe measure the minor axis PV wiggle magnitude of the 3600 models and present the resulting magnitude map in Figure 4c. By drawing contours in the Figure 4c map at the magnitude values measured for AB Aur ( 37 . 4 ± 2 . 9 m/s in 13 CO and 44 . 2 ± 1 . 3 m/s in C 18 O), we find every combination of M disk /M ⋆ and β that satisfy the observations. Repeating this procedure with our synthetic ALMA observations of the SPH GI disk simulation shown in Figure 3c, we find that this technique successfully recovers the disk mass set in the underlying SPH simulation (Extended Data Figure 7).  \nFor independent physical estimates of plausible β values between 1 '' to 5 '' ( 155 to 780 au), we rely on radiative cooling prescriptions 63,64 . From Equation 39 of Zhang & Zhu (2020) 64 , β is a function of r and depends on M disk through the surface density Σ . We assume T = ( ϕL ⋆ 8 π r 2 σ SB ) 1 / 4 , where σ SB is the StefanBoltzmann constant, L ⋆ = 59 L ⊙ is the stellar luminosity of AB Aur 43 , and ϕ = 0 . 02 represents the flaring angle 65 . We use the DSHARP Rosseland mean opacity 66 κ R = κ R ( T, a max ) for a power-law grain size distribution truncated at a max . We set a 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-)mm emission at these large radii. We compute a β ( r ) profile for each M disk /M ⋆ ∈ [0 . 0 , 0 . 4] and extract the values at 1 '' and 5 '' . We overlay the resulting β ( M disk /M ⋆ ) ranges as white shaded regions in Extended Data Figure 8 (where the dependence on p arises from the dependence on Σ ), and in Figure 4 as white horizontal bars at a selection of M disk /M ⋆ values. For example, for M disk /M ⋆ = 0 . 2 and p = 1 . 0 , we find β (1 '' ) = 5 . 3 and β (5 '' ) = 3 . 6 × 10 -2 . While knowledge of cooling in disks is very limited, these estimates help to emphasize that not all values of β are equally likely.  \nData 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 . This includes final reduced and calibrated ALMA measurement sets, image cubes and moment maps, and processed SPHERE data. All simulated data products including hydrodynamic simulations and synthetic ALMA data are available at https://doi.org/10.5281/zenodo.11668694 . The raw ALMA data are publicly available via the ALMA archive https://almascience.nrao.edu/aq/ under project ID 2021.1.00690.S. The raw VLT/SPHERE data are publicly available via the ESO Science Archive Facility https://archive.  \neso.org/eso/eso\\_archive\\_main.html under programme 0104.C-0157(B).  \nCode availability ALMA data reduction and imaging scripts are available at https://jjspeedie. github.io/guide.2021.1.00690.S . The Python packages used in this work are available: bettermoments ( https://github.com/richteague/bettermoments ), eddy ( https://github.com/richteague/ eddy ), giggle v0 ( http://doi.org/10.5281/zenodo.10205110 ), PHANTOM ( https:// github.com/danieljprice/phantom ), MCFOST ( https://github.com/cpinte/mcfost ).  \nAcknowledgements We thank our referees for their careful and insightful comments that improved the manuscript. We thank Kaitlin Kratter for enlightening discussions and valuable suggestions. J.S. thanks Ryan Loomis, Sarah Wood and Tristan 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 funded by the NAASC. The reduction and imaging of the ALMA data was performed on NAASC computing facilities. J.S. thanks Christophe Pinte, Daniel Price and Josh Calcino for support with MCFOST, Luke Keyte and Francesco Zagaria for discussions on self-calibrating ALMA data, and Chris White for sharing perceptually uniform colormaps. J.S. acknowledges financial support from the Natural Sciences and Engineering Research Council of Canada (NSERC) through the Canada Graduate Scholarships Doctoral (CGS D) program. R.D. acknowledges financial support provided by the Natural Sciences and Engineering Research Council of Canada 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 Sklodowska-Curie grant agreement # 823823 (RISE DUSTBUSTERS project). C.L. acknowledges funding from UK Science and Technology research Council (STFC) via the consolidated grant ST/W000997/1. B.V. acknowledges funding from the ERC CoG project PODCAST No 864965. Y.W.T. acknowledges support through NSTC grant 111-2112-M-001-064- and 112-2112-M-001-066-. J.H. was supported by JSPS KAKENHI Grant Numbers 21H00059, 22H01274, 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 Astronomy 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 Canadian Foundation for Innovation.  \nAuthor 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.  \nCompeting Interests The authors declare that they have no competing financial interests.  \nCorrespondence Correspondence and requests for materials should be addressed to: J.S. (email: jspeedie@uvic.ca), R.D. (email: rbdong@uvic.ca).", 'Extended Data': "Extended Data Figure 1: Moment maps: AB Aur observations and GI disk simulations. (ac) Integrated intensity (moment 0), intensity-weighted mean velocity (moment 1), and intensityweighted line width (moment 2) maps for the ALMA 13 CO observations toward AB Aur. Panel (b) appears in the main text as Figure 2a. (d-f) Moment 0, 1, and 2 maps for the synthetic ALMA 13 COobservations of the SPH GI disk simulation. Like the AB Aur observations, the simulated GI disk displays a prominent GI wiggle along the southern minor axis (indicated by white arrows).  \n<!-- image -->  \nExtended Data Figure 2: Filtered moment maps: AB Aur observations and GI disk simulations. Expanding kernel filter residuals of the maps shown in Extended Data Figure 1, highlighting global spirals and velocity disturbances generated by GI. Panels (a-c) appear in the main text as Figure 1b-d. The minor axis GI wiggle indicated by white arrows in Extended Data Figure 1b and e is shown here as an isovelocity contour at v los = v sys ± v chan in all panels.  \n<!-- image -->  \nExtended Data Table 1: Details of the ALMA Band 6 observations.  \nExtended Data Figure 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 los map of the 2D analytic GI disk model (shown in Figure 2b). (b) v los map of the 2D analytic Keplerian disk model (shown in Figure 2b inset). (c) Synthetic ALMA 13 CO moment 1 map for the 3D SPH GI disk simulation (shown in Figure 2c). (d) Synthetic ALMA 13 CO moment 1 map for the 3D SPH Keplerian disk simulation (shown in Figure 2c inset). (e) Observed ALMA 13 CO moment 1 map for the AB Aur disk, imaged with robust 0 . 5 (shown in Figure 2a). (f) Like (e), but imaged with robust 1 . 5 .  \n<!-- image -->  \nExtended Data Figure 4: Obtaining velocity residuals in the AB Aur disk. (a) ALMA 13 CO moment 1 map, imaged with robust 0 . 5 , as shown in Figure 2a. (b) Background model made with a Keplerian rotation profile, assuming a geometrically thin axisymmetric disk (Eqn. 1). (c) Velocity residuals after subtracting the model in panel (b). Global spiral substructure is visible, but unevenly so. The model does not capture the non-axisymmetric emission surface morphology and/or superKeplerian rotation. (d) Background model made with the expanding kernel filter (Eqn. 15). (e) Velocity residuals after subtracting the model in panel (d), as shown in Figure 1b. (f-j) Like (a-e) but with the ALMA C 18 O moment 1 map, imaged with robust 1 . 5 .  \n<!-- image -->  \nExtended Data Figure 5: Kinematics of GI-driven spiral arms. (a-d) 2D analytic modeling (Longarini et al. 2021) 35 . (e-h) Synthetic ALMA 13 CO observations of the 3D SPH GI disk simulation. (i-l) ALMA observations of the AB Aur disk. (a) Disk surface density (Eqn. 4 & 5). (b) Line-of-sight velocity (Eqn. 2), as in Figure 2b. (c) Velocity residuals from Keplerian (i.e., subtracting Figure 2b inset). (d) Line-of-sight component of the radial velocity (first term of Eqn. 2). (e) Filtered moment 0. (f) Moment 1. (g) Moment 1 residuals from Keplerian. (h) Filtered moment 1. (i) ALMA 13 CO filtered moment 0. (j) ALMA 13 CO filtered moment 2. (k) ALMA 13 CO filtered moment 1. (l) ALMA C 18 O filtered moment 1 (robust 1 . 5 ). Panel (l) inset overlays the VLT/SPHERE H -band scattered light spirals S1-S7 (ref. 14,67 ) in red, and 13 CO spirals S1-S9 we identify in black. 29  \n<!-- image -->  \nExtended Data Figure 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 centers from Figure 3a with linear and quadratic trendlines 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 γ = 0 . 25 (Equation 15). We find the background trendlines by extracting the velocity values from the high-pass filter background map (e.g. Extended Data Figure 4d) within the same 0 . 5 · -wide wedge-shaped mask as we do for the line centers, positioned along the southern disk minor axis. (c) Like (b), but varying γ 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.  \n<!-- image -->  \n<!-- image -->  \nExtended Data Figure 7: PV wiggle morphology, magnitude, and disk mass recovery in the SPH GI disk simulation. Like Figure 4, but for the synthetic ALMA observations of the SPH GI disk simulation. (a) The synthetic ALMA 13 CO line centers along the southern minor axis from Figure 3c, after quadratic detrending. Uncertainties on the line centers are shown by yellow shaded regions. The magnitude of this PV wiggle is measured to be 39 . 1 ± 1 . 9 m/s. The analytic model shown in the background for qualitative comparison has the same parameters as the underlying SPH simulation ( M disk /M ⋆ = 0 . 29 and β = 10 ) and its PV wiggle magnitude is 39 . 0 m/s. (b) As in Figure 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.  \n<!-- image -->  \nExtended Data Figure 8: Comparisons to additional sets of analytic models. Like Figure 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 13 CO and C 18 O line centers (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 disk /M ⋆ = 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 disk /M ⋆ 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 13 CO and C 18 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 Figure 4c is m = 3 , p = 1 . 0 . 32  \nExtended Data Figure 9: Candidate sites of planet formation. Coloured × 's mark the locations of candidate protoplanets reported in the literature 25,43,48 . A table providing the candidates' coordinates on the sky, estimated masses and the reporting references is available as source data . (a) Filtered ALMA 13 CO moment 0 map, as in Figure 1c. (b) VLT/SPHERE H -band scattered light image (ref. 48 ), as in Figure 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 12 CO J = 2 -1 moment 0 (ref. 25 ) at levels from 25 to 50 mJy/beam km/s in increments of 5 mJy/beam km/s.  \n<!-- image -->"}