# Molecules with ALMA at Planet-forming Scales (MAPS). XIX. Spiral Arms, a Tail, and Diffuse Structures Traced by CO around the GM Aur Disk

Jane Huang, Edwin A. Bergin, Karin I. Öberg, Sean M. Andrews, Richard Teague, Charles J. Law, Paul Kalas, Yuri Aikawa, Jaehan Bae, Jennifer B. Bergner, et al.

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

► To cite this version:

Jane Huang, Edwin A. Bergin, Karin I. Öberg, Sean M. Andrews, Richard Teague, et al. Molecules with ALMA at Planet-forming Scales (MAPS). XIX. Spiral Arms, a Tail, and Diffuse Structures Traced by CO around the GM Aur Disk. The Astrophysical Journal Supplement Series, 2021, 257, 10.3847/1538-4365/ac143e. insu-03672362

HAL Id: insu-03672362

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

Submitted on 19 Sep 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_2.png" /></figure>

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

# Molecules with ALMA at Planet-forming Scales (MAPS). XIX. Spiral Arms, a Tail, and Diffuse Structures Traced by CO around the GM Aur Disk

Jane Huang<sup>1,2,25</sup> (D), Edwin A. Bergin<sup>1</sup> (D), Karin I. Öberg<sup>2</sup> (D), Sean M. Andrews<sup>2</sup> (D), Richard Teague<sup>2</sup> (D), Charles J. Law<sup>2</sup> (D), Paul Kalas<sup>3,4,5</sup> (D), Yuri Aikawa<sup>6</sup> (D), Jaehan Bae<sup>7,8,25</sup> (D), Jennifer B. Bergner<sup>9,25</sup> (D), Alice S. Booth<sup>10,11</sup> (D), Arthur D. Bosman<sup>1</sup> (D), Jenny K. Calahan<sup>1</sup> (D), Gianni Cataldi<sup>6,12</sup> (D), L. Ilsedere Cleeves<sup>13</sup> (D), Ian Czekala<sup>3,14,15,16,17,25</sup> (D), John D. Ilee<sup>11</sup> (D), Romane Le Gal<sup>2,18,19,20</sup> (D), Viviana V. Guzmán<sup>21</sup> (D), Feng Long<sup>2</sup> (D), Ryan A. Loomis<sup>22</sup> (D), François Ménard<sup>19</sup> (D), Hideko Nomura<sup>12</sup> (D), Chunhua Qi<sup>2</sup> (D), Kamber R. Schwarz<sup>23,25</sup> (D), Takashi Tsukagoshi<sup>12</sup> (D), Merel L. R. van 't Hoff<sup>1</sup> (D), Catherine Walsh<sup>11</sup> (D), David J. Wilner<sup>2</sup> (D), Yoshihide Yamato<sup>6</sup> (D), and Ke Zhang<sup>1,24,26</sup> (D)

## Abstract

The concentric gaps and rings commonly observed in protoplanetary disks in millimeter continuum emission have lent the impression that planet formation generally proceeds within orderly, isolated systems. While deep observations of spatially resolved molecular emission have been comparatively limited, they are increasingly suggesting that some disks interact with their surroundings while planet formation is underway. We present an analysis of complex features identified around GM Aur in $^{12}$CO J = 2 − 1 images at a spatial resolution of ~40 au. In addition to a Keplerian disk extending to a radius of ~550 au, the CO emission traces flocculent spiral arms out to radii of ~1200 au, a tail extending ~1800 au southwest of GM Aur, and diffuse structures extending from the north side of the disk up to radii of ~1900 au. The diffuse structures coincide with a “dust ribbon” previously identified in scattered light. The large-scale asymmetric gas features present a striking contrast with the mostly axisymmetric, multi-ringed millimeter continuum tracing the pebble disk. We hypothesize that GM Aur’s complex gas structures result from late infall of remnant envelope or cloud material onto the disk. The morphological similarities to the SU Aur and AB Aur systems, which are also located in the L1517 cloud, provide additional support to a scenario in which interactions with the environment are playing a role in regulating the distribution and transport of material in all three of these Class II disk systems. This paper is part of the MAPS special issue of the Astrophysical Journal Supplement.

*Unified Astronomy Thesaurus concepts: Protoplanetary disks (1300); Planet formation (1241); Radio interferometry (1346); Astrochemistry (75)*

## 1. Introduction

Observing the physical and chemical conditions in protoplanet-
tary disks is key for elucidating how planets form. Carbon
monoxide (CO) has long been one of the most important probes
of disk structure and evolution because of its high abundance in

disks and the existence of readily detectable rotational transitions
at millimeter wavelengths. Early interferometric CO detections
toward pre-main-sequence stars established the presence of gas-
rich, Keplerian disks (e.g., Sargent & Beckwith 1987; Koerner
et al. 1993). Disks with detections of bright, extended Keplerian
CO emission naturally became popular targets for chemical
studies at millimeter wavelengths, initially with single-dish
telescopes (e.g., Dutrey et al. 1997; Thi et al. 2004) and then
with interferometers (e.g., Qi et al. 2003; Öberg et al. 2010).

<footnote n="1">¹Department of Astronomy, University of Michigan, 233 West Hall, 1085 S. University Ave., Ann Arbor, MI 48109, USA; jnhuang@umich.edu</footnote>

<footnote n="2">²Center for Astrophysics | Harvard & Smithsonian, 60 Garden St., Cambridge, MA 02138, USA</footnote>

<footnote n="3">³Department of Astronomy, 501 Campbell Hall, University of California, Berkeley, CA 94720-3411, USA</footnote>

<footnote n="4">⁴SETI Institute, Carl Sagan Center, 189 Bernardo Ave., Mountain View, CA 94043, USA</footnote>

<footnote n="5">⁵Institute of Astrophysics, FORTH, GR-71110 Heraklion, Greece</footnote>

<footnote n="6">⁶Department of Astronomy, Graduate School of Science, The University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-0033, Japan</footnote>

<footnote n="7">⁷Earth and Planets Laboratory, Carnegie Institution for Science, 5241 Broad Branch Road NW, Washington, DC 20015, USA</footnote>

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

<footnote n="9">⁹Department of Geophysical Sciences, University of Chicago, Chicago, IL 60637, USA</footnote>

<footnote n="10">¹⁰Leiden Observatory, Leiden University, 2300 RA Leiden, The Netherlands</footnote>

<footnote n="11">¹¹School of Physics & Astronomy, University of Leeds, Leeds LS2 9JT, UK</footnote>

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

<footnote n="13">¹³University of Virginia, 530 McCormick Rd, Charlottesville, VA 22904, USA</footnote>

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

<footnote n="15">¹⁵Center for Exoplanets and Habitable Worlds, 525 Davey Laboratory, The Pennsylvania State University, University Park, PA 16802, USA</footnote>

<footnote n="16">¹⁶Center for Astrostatistics, 525 Davey Laboratory, The Pennsylvania State University, University Park, PA 16802, USA</footnote>

<footnote n="17">¹⁷Institute for Computational & Data Sciences, The Pennsylvania State University, University Park, PA 16802, USA</footnote>

<footnote n="18">¹⁸IRAP, Université de Toulouse, CNRS, CNES, UT3, Toulouse, France</footnote>

<footnote n="19">¹⁹Univ. Grenoble Alpes, CNRS, IPAG, F-38000 Grenoble, France</footnote>

<footnote n="20">²⁰IRAM, 300 rue de la piscine, F-38406 Saint-Martin d'Hères, France</footnote>

<footnote n="21">²¹Instituto de Astrofísica, Pontificia Universidad Católica de Chile, Av. Vicuña Mackenna 4860, 7820436 Macul, Santiago, Chile</footnote>

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

<footnote n="23">²³Lunar and Planetary Laboratory, University of Arizona, 1629 E. University Blvd, Tucson, AZ 85721, USA</footnote>

<footnote n="24">²⁴Department of Astronomy, University of Wisconsin-Madison, 475 N. Charter St, Madison, WI 53706, USA</footnote>

<footnote>Received 2021 February 19; revised 2021 June 10; accepted 2021 June 29; published 2021 November 3</footnote>

<footnote n="25">25 NASA Hubble Fellowship Program Sagan Fellow.</footnote>

<footnote n="26">²⁶ NASA Hubble Fellow.</footnote>
Five such disks (IM Lup, GM Aur, AS 209, HD 163296, and MWC 480) constitute the sample for Molecules with ALMA at Planet-forming Scales (MAPS), an ALMA Cycle 6 Large Program that conducted the first high-angular-resolution millimeter wavelength survey of protoplanetary disk chemistry (Öberg et al. 2021). These sources were selected because of the rich and varied substructures that have previously been imaged with the Atacama Large Millimeter/submillimeter Array (ALMA) in a variety of molecules (e.g., Huang et al. 2017; Bergner et al. 2019; Pegues et al. 2020) as well as in millimeter continuum emission (e.g., Andrews et al. 2018; Huang et al. 2018; Long et al. 2018; Huang et al. 2020a).

The distribution and transportation of gas are key inputs for models of disk evolution and planet formation. Simulations of the formation of dust gaps and rings via planet-disk interactions typically adopt simple, axisymmetric gas surface density profiles, with no exchange of material with their environments (e.g., Dipierro et al. 2015; Bae et al. 2017; Zhang et al. 2018). While deep molecular line observations support this treatment for some widely modeled systems, including most of the MAPS targets, recent observations of systems such as HL Tau, RU Lup, and BHB1 have shown that concentric gap and ring dust structures can coexist with large-scale, non-Keplerian stream or spiral structures associated with the disk (e.g., ALMA Partnership et al. 2015; Yen et al. 2017; Andrews et al. 2018; Yen et al. 2019; Alves et al. 2020; Huang et al. 2020b). However, since most ALMA CO disk observations have been perfunctory, with integration times on the order of a couple minutes (e.g., Barenfeld et al. 2016; Ansdell et al. 2018), the range and prevalence of possible gas structures within and around disks are not yet well-characterized.

In this article, we analyze complex CO features identified in MAPS observations of the disk around GM Aur (International Celestial Reference System 04:55:10.981, +30:21:59.376), a K5.5 T Tauri star at a distance of 159 pc in the Taurus–Auriga star-forming region (e.g., Espaillat et al. 2010; Bailer-Jones et al. 2018; Gaia Collaboration et al. 2018). As one of the first disks to be identified as “transitional” (i.e., with a central dust cavity) and to be spatially resolved in CO emission, it has been subject to extensive observations from centimeter to X-ray wavelengths in order to probe the disk structure, chemistry, accretion behavior, and signatures of planet-disk interactions (e.g., Strom et al. 1989; Koerner et al. 1993; Bergin et al. 2004; Calvet et al. 2005; Hughes et al. 2009; Espaillat et al. 2019). Much like most other disks that have now been imaged at moderate-to-high resolution, the GM Aur disk’s millimeter continuum displays a series of nearly axisymmetric gaps and rings, (e.g., Huang et al. 2020a). Meanwhile, low-resolution (~2") $^{12}$CO observations have traced large-scale asymmetric structures deviating from the expected Keplerian velocity field of the disk (Dutrey et al. 2008; Hughes et al. 2009). The MAPS program improves upon previously published GM Aur $^{12}$CO observations by an order of magnitude in angular resolution, revealing a tremendously intricate system with flocculent spiral arms, a tail, and diffuse structures. This large-scale complexity sets GM Aur apart from the other MAPS disks, which exhibit nearly axisymmetric CO emission (Law et al. 2021a). The GM Aur observations and data reduction are presented in Section 2, the molecular features are described in Section 3, their possible origins are discussed in Section 4, and our conclusions are summarized in Section 5.

## 2. Observations and Data Reduction

### 2.1. ALMA observations

A full description of the observational setup and self-calibration procedure for the MAPS Large Program is provided in Öberg et al. (2021).<sup>27</sup> We briefly summarize the key aspects. GM Aur was observed in the spectral setting containing $^{12}$CO $J=2-1$ over the course of eight execution blocks shared with MWC 480. Three were taken in the C43-4 configuration (baselines ranging from 14 to 1397 m) in 2018 October and November for a combined on-source integration time of one hour. The other five were taken in the C43-7 configuration (baselines ranging from 40 to 3638 m) in 2019 August, for a combined on-source integration time of two hours. Following self-calibration, a set of fiducial images was produced at a spatial resolution of 0.′′15 and spectral resolution of 0.2 km s$^{-1}$ for all Band 6 lines covered in the survey (Czekala et al. 2021). The fiducial imaging procedure aimed to facilitate high-resolution studies of the radial and vertical distributions of molecular emission. Analysis of the structure and chemistry of GM Aur’s Keplerian disk region as traced by CO isotopologues is presented in Law et al. (2021a, 2021b), Zhang et al. (2021), Bosman et al. (2021a), and Schwarz et al. (2021).

However, GM Aur’s $^{12}$CO emission is visible in some channels up to ~12" from the phase center. The ALMA Technical Handbook<sup>28</sup> suggests two definitions for maximum recoverable scale (MRS): one calculated using the shortest baseline $L_{min}$ ($\theta_{MRS} \approx 0.6\lambda_{obs}/L_{min}$), and the more conservative one using the fifth percentile of the baseline lengths ($\theta_{MRS} \approx 0.983\lambda_{obs}/L_5$). Given these definitions, the estimated MRS of GM Aur’s CO emission ranges from ~2" to 11". To improve sensitivity to large-scale emission, we reimaged $^{12}$CO $J=2-1$ with a Gaussian uv taper and a robust value of 1.0, resulting in a synthesized beam of 0.′′28 × 0.′′23 (−14°4). Unless otherwise specified, figures in this article are based on this version of the $^{12}$CO image, which we will subsequently refer to as the MAPS XIX version. The reimaging was performed with CASA 5.6.1 (McMullin et al. 2007). We applied the multiscale CLEAN algorithm (Cornwell 2008) with scales of [0, 0.′′4, 1", 2"). Since $^{12}$CO was also reobserved at higher spectral resolution than the other lines, we reimaged $^{12}$CO with channel widths of 0.1 km s$^{-1}$ for better recovery of the kinematic details. Due to the irregular emission morphology, we used CASA’s auto-multithresh algorithm (Kepley et al. 2020) to draw the CLEAN mask. The auto-multithresh algorithm searches the cube for significant emission, beginning with a relatively conservative mask and then expanding to encompass more emission during subsequent major cycles. The mask was initialized with full coverage of the primary beam from 5.2–6.4 km s$^{-1}$, where the emission is the broadest, because auto-multithresh algorithm does not readily mask diffuse emission. After some experimentation, the following auto-multithresh parameters were selected: sidelobethreshold=3.0, noisethreshold=4.0, ldownoisethreshold=1.5, and minbeamfrac=0.3. The CLEAN threshold was set to 5 mJy, corresponding to ~3× the rms of line-free channels in the dirty image.

<footnote n="27"><sup>27</sup> Self-calibration and imaging scripts, as well as imaging products, are available at http://alma-maps.info.</footnote>

<footnote n="28"><sup>28</sup> https://almascience.nrao.edu/documents-and-tools/cycle7/alma-technical-handbook/view</footnote>
As noted in Jorsater & van Moorsel (1995), the default procedure at the end of CLEAN is to add the residuals (in units of janskys per dirty beam) to the CLEAN model convolved with the synthesized beam (in units of janskys per clean beam) under the assumption that these units are roughly equivalent. However, this assumption does not hold when the dirty beam deviates significantly from the Gaussian shape of the synthesized beam, as is the case for the MAPS data set. Throughout the MAPS paper series, this unit mismatch is referred to as the “JvM effect.” To counteract this issue, the MAPS Large Program applied a “JvM correction” to all image products, where the residuals are scaled by a factor of $\epsilon$ = clean beam volume/dirty beam volume before being added to the convolved CLEAN model (Czekala et al. 2021). As noted in Czekala et al. (2021), lowering the CLEAN threshold does not obviate the need for the JvM correction, because CLEAN would then tend erroneously to incorporate more noise into the CLEAN model—and also because the units of the residuals would still be inconsistent with those of the convolved CLEAN model. We refer the reader to Czekala et al. (2021) for full details on the MAPS imaging workflow and the impact of the JvM effect and correction on the MAPS data sets. For the reimaged GM Aur $^{12}$CO observations presented in this work, $\epsilon$ = 0.46. A primary beam correction was applied to the resulting image using *impcor* in CASA. The resulting rms of the image cube, measured within a circle with a radius of 10" in line-free channels, is 1 mJy beam$^{-1}$.

The fiducial images for $^{13}$CO $J=2-1$, which was observed simultaneously with $^{12}$CO $J=2-1$, exhibited hints of residual structure outside its Keplerian CLEAN mask. While we used the auto-multithresh algorithm for $^{12}$CO because its emission is bright and irregular, we chose to use a simple circular mask with a diameter of 11" to reimage $^{13}$CO in order to check that noise fluctuations were not being amplified by selective masking. The CLEAN threshold was set to 4 mJy, corresponding to ~3× the rms of line-free channels in the dirty image. To increase sensitivity, a Gaussian uv taper and robust value of 1.0 were applied, resulting in a synthesized beam of $0.39 \times 0.35$ (−4°.4). Since the $^{13}$CO spectral windows had coarser spectral resolution than $^{12}$CO, the $^{13}$CO image cube was produced with a channel spacing of 0.2 km s$^{-1}$. In other respects, the $^{13}$CO reimaging followed the same procedure as the $^{12}$CO reimaging. The resulting rms of the image cube, following residual scaling corrections with $\epsilon$ = 0.60, is 1 mJy beam$^{-1}$.

## 2.2. Archival Hubble Space Telescope Observations

For comparison with the ALMA data, Wide Field Planetary Camera 2 (WFPC2) observations of GM Aur from program HST-GTO/WF2-6223 (PI: Trauger) were retrieved from the Hubble Space Telescope (HST) archive (now the Mikulski Archive for Space Telescopes). Two exposures lasting 200 seconds each (data sets U2RD0406T and U2RD0409T) were taken on 1995 July 29 in the F675W filter (central wavelength: 6717 Å) with the PC1 camera, which has a pixel scale of 0.0455. The data available on the archive were processed with the standard OPUS pipeline (Swam & Swade 1999). PSF subtraction was performed on GM Aur using the reference star HD 283572, which was observed as part of the same program. An independent reduction of these data was previously presented in Stapelfeldt et al. (1997). While HST images of GM Aur have been published at other wavelengths, we use the WFPC2 F675W image as the basis for comparison because the WFPC2 images presented in Hornbeck et al. (2016)

had shorter exposure times, and the field of view of the NICMOS coronagraphic images presented in Schneider et al. (2003) is more limited.

# 3. Results

## 3.1. Complex Structures Traced by $^{12}$CO

Figure 1 shows integrated intensity, intensity-weighted velocity, and peak brightness temperature maps of $^{12}$CO $J=2-1$ emission toward GM Aur, as well as the high-resolution 1.1 mm continuum image from Huang et al. (2020a) for comparison. The integrated intensity map includes channels with LSRK velocities from −1.6 to 12.5 km s$^{-1}$, based on the velocity range where disk emission is detected above the $3\sigma$ level. Because emission at any given position spans a much narrower velocity range than the overall velocity range encompassed by the integrated intensity map (see the $^{12}$CO $J=2-1$ channel maps in Appendix A), only pixels in individual channels above the $3\sigma$ level are included in the integrated intensity map, in order to reduce noise contributions from signal-free regions. The fiducial integrated intensity maps presented in Law et al. (2021a) use modified Keplerian masks, but this approach is not suitable for GM Aur’s $^{12}$CO emission, due to the extended, non-Keplerian structures.

To produce a peak brightness temperature map, we imaged the $^{12}$CO line from visibilities that have not been continuum-subtracted, in order to avoid artificially reducing the intensities in regions where continuum emission is absorbed by optically thick $^{12}$CO (e.g., Weaver et al. 2018). While the large-scale features of interest in this work lie outside the millimeter continuum and are therefore unaffected by continuum subtraction, we elect to make peak brightness temperature maps without continuum subtraction, in order to obtain results consistent with Law et al. (2021b). We otherwise follow the imaging steps described in Section 2. Similar image cubes were also produced with the channel gridding shifted by +0.03 and −0.03 km s$^{-1}$. A map is generated for each cube by calculating the peak intensity as a function of velocity at each pixel. Emission below the $5\sigma$ level is masked. The three maps were then median-stacked to mitigate channelization artifacts (e.g., Christiaens et al. 2014). The median peak intensity map was then converted to brightness temperatures using the Planck equation (i.e., not the Rayleigh-Jeans approximation). The intensity-weighted velocity map was created with the same velocity range as the other two maps, with pixels below the $5\sigma$ level masked. Higher masking thresholds were used for the peak intensity and intensity-weighted velocity maps compared to the integrated intensity map because the former are more sensitive to noisy outliers.

The $^{12}$CO emission is broad and asymmetric, extending up to a distance of ~12" (~1900) au from the star. Several arm-like structures emerge from the edge of the disk, more prominently in the peak brightness temperature map compared to the integrated intensity map. A long “tail” extends from the southwest side of the disk. Diffuse emission surrounds the northern side of the disk. Spatial filtering afflicts the appearance of the largest-scale $^{12}$CO emission by creating negative CLEAN “bowls” in channels near the systemic velocity (5.4–6.0 km s$^{-1}$). However, the appearance of the extended structures in the central channels is broadly consistent with the emission morphology observed with the Submillimeter Array’s more compact antenna configurations (Hughes et al. 2009). The irregularity and large extent of the $^{12}$CO emission presents a stark contrast with the continuum
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Moment maps of $^{12}$CO J = 2 − 1 toward GM Aur and a comparison with the millimeter continuum. Top left: An integrated intensity map of $^{12}$CO J = 2 − 1 emission toward GM Aur. A logarithmic color stretch is used to highlight faint structures. The scale of the synthesized beam, 0.′′28 × 0.′′23 (−14.′′4), is represented by the white ellipse in the lower left corner. Top right: A peak brightness temperature map of $^{12}$CO J = 2 − 1. A logarithmic color stretch is used to highlight faint structures. Bottom left: The $^{12}$CO intensity-weighted velocity map. Bottom right: The high-resolution 1.1 mm GM Aur continuum image from Huang et al. (2020a). An arcsinh color stretch is used to highlight faint structures. The scale of the synthesized beam, 45 mas × 25 mas (2.′.2), is represented by the white ellipse in the lower left corner. All images are shown on the same scale, to emphasize the different spatial extents of the molecular and dust emission. North is up and east is to the left.</figcaption></figure>

emission, which is comparatively compact (extending up to a radius of ~250 au) and nearly axisymmetric.

The $^{12}$CO velocity map traces the Keplerian disk out to a radius of ~550 au, with the northeast side redshifted and the southwest side blueshifted relative to the systemic velocity of 5.61 km s⁻¹ (Huang et al. 2020a). The extended emission exhibits more complex kinematic behavior—the southwest tail and diffuse northwest emission are slightly redshifted relative to the systemic velocity, the diffuse emission northeast of the disk becomes more blueshifted from south to north, and the southeastern edge of the $^{12}$CO emission is blueshifted relative to the Keplerian disk. The following sections examine the complex CO emission structures in more detail.

### 3.1.1. *The Spiral Arms*

Four spiral arms (S1, S2, S4, and S5) are identified based on visual inspection of GM Aur's $^{12}$CO J = 2 − 1 channel maps, along with an additional tentative arm (S3). Figure 2 shows the channel maps with logarithmic spiral curves overlaid (i.e., of the form $R(\varphi) = R_0 \exp(b\varphi)$, where R is the distance of the spiral from GM Aur, $\varphi$ is the azimuthal angle in radians, and $R_0$ and b are free parameters). The coordinate system is defined to be in the plane of a geometrically thin disk with a position angle of 57.°17 east of north and inclination of 53.°21 (Huang et al. 2020a). The

northeast side of the disk major axis corresponds to $\varphi = \frac{\pi}{2}$, and $\varphi$ increases clockwise (see Figure 10 from Huang et al. (2020a) for a schematic). Table 1 lists the parameter values of the overlaid curves, which were chosen to produce a reasonable visual match to the channel maps. A parameterization with constant pitch angle (i.e., a logarithmic spiral) is adopted for simplicity, given the relatively low signal-to-noise ratio of the spiral arm emission (typically ranging from 5–10$\sigma$ in a given location) and the small angular range over which the arms are detected (usually less than half a winding).

In general, a feature is considered to be detected if its emission exceeds the 4$\sigma$ level in more than one channel. Nevertheless, because of the patchiness of some of the emission and blending of signal between the arms and with the Keplerian disk in some channels, reasonable alternative spiral identifications may be possible. In particular, S3 is the least certain because of the large discontinuities between detected emission features. While the tapering performed for the MAPS XIX version of the $^{12}$CO image generally makes the faint emission from the spiral arms more readily apparent in the channel maps, the lower resolution obscures some of the emission from S4 and S5 at smaller disk radii. The interior regions of these arms are more readily discernible in the higher-resolution fiducial $^{12}$CO channel maps (Figure 3).
<figure><img src="image_5.png" /><figcaption>**Figure 2.** $^{12}$CO channel maps toward GM Aur with pink curves overlaid on the gas spiral structures. The dotted curve for S3 denotes less certainty in its identification. The light blue contours show the 5σ, 15σ, and 25σ emission levels, where σ = 1 mJy beam$^{-1}$. The isovelocity contours of the Keplerian disk are drawn in light gray, with the solid curves tracing the front of the disk and the dotted curves tracing the back side. The upper right corner of each panel is labeled with the LSRK velocity (km s$^{-1}$). The scale of the synthesized beam, 0.′′28 × 0.′′23 (−14.74), is represented by the white ellipse in the lower left corner of each panel. North is up and east is to the left.</figcaption></figure>

In addition to the spiral arms identified in Table 1, the $^{12}$CO
peak brightness temperature map exhibits features resembling
short arms on the northeast side of the disk. A comparison
between the peak brightness temperature map and the $^{12}$CO
channel maps, though, suggests that projection effects from the
back side of the disk may be responsible for these arm-like
emission features (Figure 4).

Assuming that the spiral arms are in the plane of the disk, they
are detected in $^{12}$CO at radii from ~340 to ~1200 au from the

star, or ~1900 au if the tentative S3 arm is included.
Nevertheless, it is plausible that the spiral arms do not lie in
the same plane as the disk or with one another. The analysis of
the $^{12}$CO emitting height of the Keplerian portion of the disk in
Law et al. (2021b) indicates that the $J=2-1$ emission generally
originates from well above the midplane ($z/r \sim 0.4$), although
the emitting height trends downward at larger radii and becomes
highly uncertain at the projected distances of the spiral arms.
Furthermore, the spiral arms deviate from Keplerian motion, as
<figure><img src="image_6.png" /><figcaption>Figure 2. (Continued.)</figcaption></figure>

<table><caption>**Table 1**
Parameters of Proposed Spirals</caption><thead><tr><th>ID<br>(1)</th><th>R<sub>0</sub><br>(au)<br>(2)</th><th>b<br>(3)</th><th>Pitch Angle<sup>a</sup><br>(deg)<br>(4)</th><th>&phi; Range<sup>b</sup><br>(radians)<br>(5)</th></tr></thead><tbody><tr><td>S1</td><td>1860</td><td>0.45</td><td>24</td><td>(-0.9&pi;, -0.3&pi;)</td></tr><tr><td>S2</td><td>1235</td><td>0.4</td><td>22</td><td>(-0.7&pi;, -0.2&pi;)</td></tr><tr><td>S3?</td><td>1200</td><td>0.8</td><td>39</td><td>(-0.2&pi;, 0.18&pi;)</td></tr><tr><td>S4</td><td>825</td><td>0.4</td><td>22</td><td>(-0.5&pi;, 0.25&pi;)</td></tr><tr><td>S5</td><td>500</td><td>0.5</td><td>27</td><td>(-0.25&pi;, 0.38&pi;)</td></tr></tbody></table>

**Notes.**

shown in Figure 2, where isovelocity contours for the Keplerian disk are drawn in gray for visual guidance. Details of how the isovelocity contours are computed, with adjustments for the

pressure gradient, are provided in Appendix B. The non-Keplerian motions suggest that the arms could be oriented out of the plane of the disk. As demonstrated in Dong et al. (2016) and Uyama et al. (2020), projection effects can make even symmetric spiral arms appear to have dramatically different pitch angles and extents. Thus, while spiral arm parameters are provided for guidance, they should be interpreted with caution, e.g., variations in the listed pitch angles for different arms could in large part reflect uncertainties associated with the viewing angle.

The spiral arms' deviations from Keplerian motion are further highlighted in Figure 5 (this and all subsequent figures are based on the MAPS XIX image versions). In the intensity-weighted velocity map, the portion of S4 and S5 on the southeast side of the disk stand out as being blueshifted relative to the disk's Keplerian motion. To examine velocity deviations in more detail, a velocity residual map was created by computing a model Keplerian velocity field for the front side of the disk using the stellar properties, disk orientation, and <sup>12</sup>CO emitting height described in Appendix B, then subtracting

<footnote><sup>a</sup> The pitch angle of a logarithmic spiral is arctan *b*.</footnote>

<footnote><sup>b</sup> Approximate azimuthal range over which the spiral arm is visible in the <sup>12</sup>CO channel maps.</footnote>
<figure><img src="image_7.png" /><figcaption>**Figure 3.** Channel maps from the fiducial MAPS ¹²CO image, with curves for S4 and S5 overlaid to highlight the emission from these spirals at smaller disk radii. An arcsinh color stretch is used to help make the faint spiral features more visible. The spiral emission manifests as protrusions from the Keplerian disk. The isovelocity contours of the Keplerian disk are drawn in gray, with the solid curves tracing the front of the disk and the dotted curves tracing the back side. The upper right corner of each panel is labeled with the LSRK velocity (km s⁻¹). The synthesized beam (0″·15 × 0″·15) is shown in the lower left corner of each panel. North is up and east is to the left.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 4.** A comparison of a ¹²CO channel map with the peak brightness temperature map suggests that the apparent arm-like structures to the northeast of the disk are due to emission from the disk's back side. Left: ¹²CO emission at 6.6 km s⁻¹, with the front and back sides of the Keplerian disk emission marked with solid and dashed contours, respectively. Right: ¹²CO peak brightness temperature map with the isovelocity contours corresponding to v<sub>lsr</sub> = 6.6 km s⁻¹ overlaid. The back side of the disk overlaps with the apparent arm-like structures.</figcaption></figure>

it from the intensity-weighted velocity map. Unlike the isovelocity contour calculations used for the channel maps, the pressure gradient is not included in the velocity residual calculations, because Equation (B2) yields imaginary solutions for $v_\phi$ at large distances from the star (i.e., in the vicinity of the diffuse structures north of the disk) and would result in a discontinuous residual map. Overall, the residual map is quite complex, not only because of the spiral arms, but also because of the southwest tail and diffuse northern emission. The non-Keplerian motion of the spiral arms can be traced to smaller radii in the velocity residual map compared to the intensity-weighted velocity map. In general, S4 and S5 are blueshifted relative to the expected Keplerian velocities, while S1 and S2 are redshifted.

Figure 6 shows the ¹²CO velocity dispersion map, which was created with the CASA *immoments* task. Pixels below the 4σ level were clipped in order to avoid skewing from signal-free channels. For the most part, the spiral arms do not stand out prominently in their velocity dispersion, which is often computed to search for regions of increased turbulence.

At smaller radii, S4 coincides with protrusions in the contour levels, which may hint at enhanced velocity dispersion in its vicinity. The southwest side of the disk, which overlaps with S1 and S2, exhibits more clearly elevated velocity dispersion, which can be attributed to the blending of multiple velocity components, including the two arms, the Keplerian disk, and the southwest tail (discussed further in Section 3.1.2).

### 3.1.2. The Southwest Tail

The southwest tail is visible in the ¹²CO channel maps from ~5.3–6.0 and 6.3–6.8 km s⁻¹. In between these velocity ranges, only broad emission is observed, implying either that the diffuse northern emission component is in front of the tail at these velocities or that they are blended due to similar brightnesses. To characterize the geometry of the tail, we produced an integrated intensity map from ~5.3–6.0 km s⁻¹, without applying any intensity cutoffs (Figure 7). The higher velocity channels are excluded in order to avoid confusion with the diffuse northern emission, but this exclusion does not substantially affect characterization of the tail’s geometry, because the tail emission at higher velocities is compact. The rms, measured in a signal-free region of this map, is 0.6 mJy beam⁻¹ km s⁻¹. To measure the orientation of the tail, we took a series of vertical cuts across the image spaced 0″·28 apart (i.e., comparable to the synthesized beam) in the east–west direction and identified local intensity maxima along these cuts. We consider the starting point of the tail to be 0″·64 west of GM Aur, which is the closest point to the star at which the tail can be distinguished from the Keplerian disk emission via separate local maxima in a vertical slice across the image. Vertical cuts were taken until the local maximum fell below 3σ, which occurs 10″·2 west of GM Aur. The uncertainty in the position of each local maximum is taken to be the standard deviation of the Gaussian synthesized beam. Since the tail does not exhibit substantial curvature, we fit a line, $y = m_{\text{tail}}x + b_{\text{tail}}$, to the positions of all the local maxima using least-squares minimization. The variables x and y denote the offsets (in arcseconds) from the center of the GM Aur disk (which is aligned with the phase center) in the eastern and northern
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Left: An integrated intensity map of ¹²CO emission from 5.3 to 6.0 km s⁻¹, highlighting emission from the southwest tail. The pink ellipse at the center of the image shows the orientation of the GM Aur disk. The red points show the positions of local maxima measured in vertical slices of the image, along with 1σ error bars. The dark blue line segment shows the best fit to the red points. The thin, light blue lines show the positions at which intensities are measured across the width of the tail. Middle: The integrated intensity of ¹²CO measured along the best-fit line to the tail, starting from the northeast end. The lavender shading shows the 1σ noise level. Right: The intensities measured along positions A and B, respectively. The measurements and 1σ error bars are plotted in black. The intensity along position A is measured within a narrower interval in order to exclude emission from other disk structures. The offset from tail center is defined with respect to the best-fit line segment drawn in the left-hand panel, with negative offsets southeast of the line and positive offsets northwest. The light blue curves show the best-fit Gaussians to the measured intensity profiles. The gray Gaussian curves in the upper right-hand side of each plot show the width of the synthesized beam.</figcaption></figure>

### 3.1.3. Diffuse Northern Emission

Diffuse emission north of GM Aur is present in the ¹²CO channel maps near the disk systemic velocity (5.61 km s⁻¹). Five filamentary components (F1, F2, F3, F4, and F5) are tentatively identified within this diffuse emission and highlighted with elliptical arcs in Figure 8. The elliptical arcs are parameterized as $x = A \cos t + x_0$, $y = B \sin t + y_0$, where x and y are the offsets (in arcseconds) east and north of GM Aur, respectively, and t is an angular variable in radians. This parameterization is selected because it reasonably approximates the morphology of the filamentary components, but other parameterizations may be similarly suitable. Table 2 lists the parameter values for each filamentary component, which are determined via visual inspection. As with the spiral structures, the patchiness of the diffuse northern emission creates some ambiguity as to which structure individual emission components should be assigned. Furthermore, a possible concern in the interpretation of the CO observations is that spatial filtering may be inducing the appearance of filamentary substructure within the diffuse northern emission. However, the identification of some of these filamentary components is supported by their alignment with structures observed in scattered light (see Section 3.2).

Figure 9 shows position–velocity diagrams extracted along the arcs drawn in Figure 8. In general, F1 and F2 exhibit decreasing offsets from the systemic velocity with increasing projected distance from GM Aur. In other words, at larger projected distances, the gas in F1 and F2 is moving more slowly along the line of sight relative to GM Aur. Meanwhile, F3, F4, and F5 do not exhibit much variation in velocity with projected distance from the star. F3 and F4 are slightly blueshifted with respect to $v_{\text{sys}}$, while F5 is slightly redshifted. The Keplerian disk is responsible for the bright second emission component visible at a projected distance of ~500 au and 6 km s⁻¹ in the diagrams for

F3 and F4. At projected distances less than ~700 au, emission from F1 and F2 is spread over a wide velocity range. The Keplerian disk contributes to some of the bright emission visible at smaller projected distances and (relatively) large velocities in the diagrams for F1 and F2, but as can be seen in Figure 8, much of the emission visible at the smaller projected distances originates from clumps at the inner tips of F1 and F2. (We refer to these as clumps based on their appearance in the channel maps, but note that the clumpiness may be due to the relatively low signal-to-noise ratio or imaging artifacts rather than the intrinsic physical properties).

The diffuse emission northwest of the disk is mostly redshifted relative to the systemic velocity, and does not exhibit as much substructure as the emission northeast of the disk. The brightest emission within the northwest component traces a broad arc (which is subsequently referred to as the “northwest ridge”), as shown in the peak brightness map created from channels between 6.0 and 7.1 km s⁻¹ and with a 5σ clip (Figure 10). The shape of the ridge was measured in a manner similar to the southwest tail, as described in Section 3.1.2. We created an integrated intensity map from 6.0 to 7.1 km s⁻¹ (with no pixels clipped) and took a series of vertical cuts across the image spaced 0.′′28 apart (i.e., roughly one synthesized beam) in the east–west direction. Local intensity maxima were identified along these cuts between 0.′′40 to 6.′′28 west of the star, and their positions relative to GM Aur were fit with a quadratic, $y = a_{\text{NW Ridge}}x^2 + b_{\text{NW Ridge}}x + c_{\text{NW Ridge}}$, where the positive x and y directions are east and north of the star, respectively. We use this parameterization because a second-order polynomial is the simplest curve that approximates the shape of the northwest ridge. Using least-squares minimization, we find that the best fit is $a_{\text{NW Ridge}} = -0.109 \pm 0.009 \text{ arcsec}^{-1}$, $b_{\text{NW Ridge}} = -0.73 \pm 0.06$, and $c_{\text{NW Ridge}} = 0.53 \pm 0.09 \text{ arcsec}^{-1}$. Figure 10 shows the best-fit curve overlaid on the integrated intensity map.
<figure><img src="image_10.png" /><figcaption>**Figure 8.** $^{12}$CO channel maps toward GM Aur with pink curves overlaid to mark tentative identifications of filamentary components within the diffuse emission north of the disk. The light blue contours show the $4\sigma$ and $8\sigma$ emission levels, where $\sigma = 1$ mJy beam$^{-1}$. The isovelocity contours of the Keplerian disk are drawn in light gray, with the solid curves tracing the front of the disk and the dotted curves tracing the back side. The upper right corner of each panel is labeled with the LSRK velocity (km s$^{-1}$). Synthesized beams are shown in the lower left corner of each panel.</figcaption></figure>

A position–velocity diagram extracted along the best-fit curve to the northwest ridge is shown at the bottom of Figure 10. At smaller velocities and projected distances within a few hundred au from GM Aur, some of the emission along the curve comes from the Keplerian disk (see Figure 20). At velocities higher than $6.3 \text{ km s}^{-1}$, though, the diagram traces only emission associated

with the ridge. As the velocity increases, the emission extent of the ridge in individual channels becomes more compact. Although the diffuse northern emission had previously been ascribed to cloud contamination by Hughes et al. (2009), the tendency for the diffuse structures to have higher line-of-sight velocities at smaller projected distances from GM Aur suggests
<table><caption>**Table 2**
Parameters of Elliptical Arcs Used to Make PV Diagrams of Tentative Filament Structures</caption><thead><tr><th>ID</th><th>A<br>(arcsec)</th><th>B<br>(arcsec)</th><th>x<sub>0</sub><br>(arcsec)</th><th>y<sub>0</sub><br>(arcsec)</th><th>t Range<sup>a</sup><br>(radians)</th></tr></thead><tbody><tr><td>F1</td><td>6.5</td><td>4</td><td>9.3</td><td>1.1</td><td>(5π/9, π)</td></tr><tr><td>F2</td><td>7</td><td>4</td><td>8.25</td><td>1.7</td><td>(5π/9, π)</td></tr><tr><td>F3</td><td>7.9</td><td>5.3</td><td>7.9</td><td>1.5</td><td>(5π/12, 11π/12)</td></tr><tr><td>F4</td><td>8.9</td><td>11</td><td>8</td><td>-1.4</td><td>(5π/9, 31π/36)</td></tr><tr><td>F5</td><td>7.2</td><td>5.6</td><td>-6.8</td><td>6.4</td><td>(-π/6, 7π/18)</td></tr></tbody></table>

**Note.**

a $t = 0$ points east and $t = π/2$ points north.

that gas is moving faster closer to the star, which is consistent with the gas motion being influenced by GM Aur's gravity. Thus, while the extended emission may technically originate from cloud material, the kinematics imply close interaction with the GM Aur system.

At several of the channels redshifted relative to $v_{sys}$, the full extent of the northern portion of the Keplerian disk is not visible (Figure 11). To demonstrate the difference from the expected appearance of the Keplerian disk, Figure 11 also shows the corresponding blueshifted channels, which would otherwise be roughly mirror images across the disk minor axis (with some small deviations due to spiral structures in the outer disk). We interpret the truncated appearance of the northern side of the Keplerian disk in redshifted channels to be due to foreground obscuration by cooler, optically thick CO gas in the diffuse northwest emission structures. The relative positions and velocities of the diffuse structures and the disk thus imply that the large-scale northwest gas is moving toward the disk; i.e., it is infalling.

### 3.1.4. Summary of CO Features

Figure 12 illustrates the positions of all identified ¹²CO features in relation to one another. While portions of these features are readily visible in the integrated intensity and peak brightness maps, others required inspection of channel maps and velocity maps in order to ascertain their positions. In projection, these features exhibit a wide range of orientations, morphologies, and extents. Although features are labeled differently depending on their appearance (e.g., the northwest ridge is broader than the northeast filamentary structures and the southwest tail does not exhibit the curvature of the other structures), viewing angle may be responsible for some of this morphological variation. Thus, different labels do not necessarily imply different physical origins. In addition, the relative positions of some features identified as separate structures could instead be part of a single larger structure, with the emission from the Keplerian disk or lack of sensitivity creating the appearance of separate features. In particular, the orientations of the southwest tail, F1, F2, and F3 are suggestive because they all lie roughly along the disk major axis. Finally, several of the structures appear to extend close to the FWHM of ALMA's primary beam (~27" at 230 GHz). Given the drop-off in sensitivity with distance from the phase center, it is likely that our data do not recover the full extent of these structures. In the future, mosaicked observations and better uv coverage will be key to fully characterizing the complexity of GM Aur's immediate environment.

<figure><img src="image_11.png" /><figcaption>**Figure 9.** Position-velocity diagrams of the tentatively identified filaments in the diffuse northern emission, extracted along the elliptical arcs drawn in Figure 8. Contours are drawn at the 4 and 8 mJy beam⁻¹ levels (corresponding to 4 and 8σ in the channel maps). The white dashed line marks GM Aur's systemic velocity.</figcaption></figure>

### 3.2. Comparison to HST Scattered-light Observations

While none of the extended ¹²CO structures have counterparts in millimeter continuum emission, which nominally traces
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Top: A ¹²CO peak brightness map from channels between 6.0 and 7.1 km s⁻¹, showing a broad “ridge” of emission northwest of the disk. The collimated component to the south is part of the southwest tail. Contours are drawn in increments of 5 mJy beam⁻¹, starting from 10 and ending at 40 mJy beam⁻¹. The blue ellipse shows the position angle and inclination of the GM Aur disk. Middle: The integrated intensity map from 6.0 to 7.1 km s⁻¹, with the measured positions of the local maxima along the northwest ridge and the 1σ error bars plotted in red. The best-fit quadratic curve to the shape of the northwest ridge is drawn in light blue. The southwest tail is marked in pink. Bottom: Position–velocity diagram extracted along the best-fit curve to the northwest ridge. Contours are drawn at [4, 8, 16, 32] mJy beam⁻¹.</figcaption></figure>

millimeter-sized dust grains in the disk midplane, some of the structures have counterparts in optical and near-infrared scattered-light images, which trace submicron-sized dust grains

that are better coupled to the gas. Figure 13 compares the ¹²CO integrated intensity map to the WFPC2 F675W image of GM Aur. F1 and F3 follow the southern and northern edges, respectively, of the elongated scattered-light feature northeast of the disk, while F2 falls roughly along the middle of the feature. This scattered-light feature was previously detected in R-band with the University of Hawaii 2.2 m telescope (Kalas 1997) and at 1.1 and 1.6 µm with HST NICMOS (Schneider et al. 2003). Schneider et al. (2003) named this feature the “blue ribbon.” Substructure within the northeast dust ribbon (e.g., distinguishing F2 from F1 and F3) is not clearly visible in any of these scattered-light images.

A point source, marked with a blue circle, is visible slightly to the southeast of the dust ribbon in the F675W image shown in Figure 13. Kalas (1997) suggested that a close encounter between the point source and GM Aur was responsible for the dust ribbon, given their apparent proximity. A search of the Vizier database indicates that the point source coincides with USNO-B1.0 1203-0070338, located at J2000 04:55:11. 894 + 30:22:03.032 (Monet et al. 2003). This source does not appear in the Gaia database (Gaia Collaboration et al. 2016, 2018), so no distance estimate appears to be available. However, the proper motion reported in the USNO-B catalog (μα = −152 ± 18 mas yr⁻¹, μδ = −124 ± 129 mas yr⁻¹) is inconsistent with the median proper motion measured by Luhman (2018) for the L1517 cloud (μα = 4.7 ± 0.8 mas yr⁻¹, μδ = −24.5 ± 1.7 mas yr⁻¹), in which GM Aur is located. Thus, we consider it unlikely that USNO-B1.0 1203-0070338 is physically in close proximity to GM Aur.

Faint nebulosity is visible northwest of the disk in the WFPC2 image, overlapping with some of the diffuse ¹²CO emission. A bright tail-like structure appears to the southwest of the disk. This feature is also visible in a WFPC2 F555W image published in Hornbeck et al. (2016), although it was not specifically described as a tail. The scattered-light tail is offset from where the tail is brightest in ¹²CO emission. Instead, as shown in a comparison between the WFPC2 F675W image and several of the ¹²CO channels, the scattered-light tail flanks the southeastern edge of the CO emission tail (Figure 14). This visual offset may be the consequence of an offset in the τ = 1 surfaces, but radiative transfer modeling would be needed to confirm. The ¹²CO tail contours appear to bend slightly at the outer tip of the scattered-light tail, suggesting some change in the physical properties. Overall, the correspondence between scattered light and ¹²CO structures indicates that the general morphologies of the complex gas features are not significantly distorted by spatial filtering.

### 3.3. Tentative ¹³CO Counterparts to ¹²CO Structures

Figure 15 shows the integrated intensity and intensity-weighted velocity maps for ¹³CO J = 2 − 1 toward GM Aur. The ¹³CO integrated intensity map was made by summing over all channels from −0.6 to 11.8 km s⁻¹, with no mask applied. The intensity-weighted velocity map was computed over the same velocity range, with pixels below the 5σ level masked. In contrast to ¹²CO, ¹³CO does not exhibit strong deviations from axisymmetry or Keplerian motion in these maps. However, the channel maps exhibit some faint non-Keplerian structures that overlap with the inner portions of S1, S2, S4, S5, the northwest ridge, and the southwest tail (Figure 16). No clear ¹³CO counterparts are identified for the northeast filamentary structures.
<figure><img src="image_13.png" /><figcaption>**Figure 11.** Top row: In several channels redshifted from $v_{\text{sys}}$, the northern portion of the Keplerian disk is partially obscured by diffuse emission structures. The blue contour shows the 40σ emission level. The Keplerian isovelocity contours are drawn in pink, with the solid curves denoting the front of the disk and the dotted curves denoting the back side. The upper right corner of each panel is labeled with the LSRK velocity in km s⁻¹. Bottom row: The corresponding blueshifted channels are shown, demonstrating that the northern portion of the Keplerian disk extends further than what is visible in the redshifted channels.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 12.** Maps of the ¹²CO integrated intensity (left) and peak brightness (right), annotated with the locations of the spiral arms, diffuse northern emission components, and the southwest tail. The “vignetting” effect at the edges of the integrated intensity map arises from the decreasing sensitivity away from the phase center, due to the primary beam.</figcaption></figure>

None of the other lines observed as part of the MAPS Large Program exhibit obvious counterparts to the ¹²CO structures (Law et al. 2021a; Öberg et al. 2021). Most notably, other bright lines observed toward GM Aur include HCO⁺ J = 1–0, CS J = 2–1, ¹³CO J = 1–0, HCN J = 3–2, and H₂CO J$_{K_aK_c}$ = 3₀₃–2₀₂. The HCO⁺ J = 3–2 line was observed separately by Huang et al. (2020a) at comparable spatial resolution and sensitivity, and also did not reveal any clear non-axisymmetric structures. The absence of counterparts identified

in deep observations of these lines underlines the importance of ¹²CO as a probe of the complex surroundings of disks.

### 3.4. Constraints on the Mass of GM Aur's Extended Gas Structures

The ¹²CO and ¹³CO J = 2–1 observations together can be used to place rough constraints on the mass of GM Aur's extended gas structures. In the optically thin, LTE limit, the
<figure><img src="image_15.png" /><figcaption>Figure 13. A comparison of GM Aur's $^{12}$CO J = 2 − 1 integrated intensity map to the WFPC2 F675W optical image. The pink curves mark the $^{12}$CO substructures that occur in close proximity to scattered-light features. The blue circle in the WFPC2 image marks the location of USNO-B1.0 1203-0070338. The linear artifacts intersecting with the disk are diffraction spike residuals. North is up and east is to the left. Both images are shown with a logarithmic color stretch to highlight the faint extended features.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Figure 14. The WFPC2 F675W image of GM Aur (grayscale) overlaid with $^{12}$CO channel map contours (5$\sigma$, 15$\sigma$, 25$\sigma$, and 35$\sigma$). The tail-like structure to the southwest of GM Aur in scattered light flanks the southeastern edge of the tail traced by $^{12}$CO emission. The $^{12}$CO LSRK velocity is labeled in the lower right corner of each panel. North is up and east is to the left.</figcaption></figure>

molecular column density can be estimated as

$$N = \frac{Q(T)}{g} \exp\left(\frac{E_u}{T}\right) \frac{4\pi}{A_{ul}hc} \int I_u dv, \quad (1)$$

where Q is the partition function, g is the upper state degeneracy, $E_u$ is the upper state energy, T is the excitation temperature, $A_{ul}$ is the Einstein A coefficient, and $\int I_u dv$ is the integral of the intensity over the velocity axis (e.g., Goldsmith & Langer 1999).

Given the nondetection of most of the extended structures in lines other than $^{12}$CO J = 2 − 1, neither the optical depth

nor the gas temperature of the extended structures is well-constrained. As shown in Figure 1, the brightness temperatures of the extended structures range from ~10–20 K, but it is not clear to what extent these brightness temperature variations reflect column density versus gas temperature variations. While these would be relatively cold gas temperatures, they are not implausible. For the sake of estimating a lower bound in the gas mass from $^{12}$CO, we adopt the optically thin approximation.

To exclude $^{12}$CO emission associated with the Keplerian disk from the mass estimate, we applied a Keplerian mask made with the keplerian_mask package (Teague 2020) to
<figure><img src="image_17.png" /><figcaption>**Figure 15.** The integrated intensity (left) and intensity-weighted velocity (right) maps for $^{13}$CO J = 2 − 1 toward GM Aur. A logarithmic color stretch is used for the integrated intensity map, to make the emission in the outer disk more visible. The synthesized beam size, 0.′′39 × 0.′′35 (−4.°4), is represented with a white ellipse in the lower left corner of each panel.</figcaption></figure>

the image cube. Details of the mask generation are provided in Section C. An integrated intensity map was produced from the masked image cube, including velocities from 3.6 to 7.1 km s⁻¹ (i.e., channels where non-Keplerian structures are visible). Since spatial filtering led to some negative artifacts in the $^{12}$CO image cube, negative pixels were excluded from the integrated intensity map. The resulting map is shown in Figure 17.

The total flux of the extended structures was then extracted from an ellipse centered on GM Aur with a position angle of 57.2, semimajor axis of $11.′′5$ (1830 au), and semiminor axis of $7.′′5$" (1190 au), as drawn in Figure 17. Since the extent of the emission is comparable to the primary beam FWHM, it is possible that GM Aur's large-scale structures continue beyond ALMA's field of view. Even within ALMA's field of view, the flux is likely underestimated due to spatial filtering. Thus, the flux measured from the ALMA observations should be treated as a lower bound. We find a total flux of 15 Jy km s⁻¹ and a mean integrated intensity of 4 mJy beam⁻¹ km s⁻¹ within the elliptical extraction region.

To convert the mean integrated intensity to a column density lower limit, we adopt an average gas temperature of 20 K, somewhat colder than the temperature that Zhang et al. (2021) derived for the CO emitting layer in the outer regions of GM Aur's Keplerian disk. The molecular constants for the $^{12}$CO J = 2 − 1 transition, which are taken from the Cologne Database for Molecular Spectroscopy (CDMS), are $Q(20\,\text{K}) = 7.6$, $g = 5$, $E_u = 16.6$ K, and $A_{ul} = 6.91 \times 10^{-7}$ s⁻¹ (Müller et al. 2001, 2005). Equation (1) yields a $^{12}$CO column density lower limit of $7 \times 10^{14}$ cm⁻². Assuming an ISM-like $^{12}$CO:H₂ number density ratio of $10^{-4}$ and that the gas is predominantly molecular hydrogen, we integrated the gas surface density over the area of the ellipse used to extract the $^{12}$CO flux and estimate that the gas mass is at least $2 \times 10^{-5}$ M☉.

Although the bulk of the extended $^{12}$CO structures are not detected in $^{13}$CO, $^{13}$CO can be used to set an upper bound on their gas mass. $I_ν$ is approximated as a Gaussian line profile with a peak of 3 mJy beam⁻¹ (i.e., 3 × the rms) and an approximately thermal FWHM of 0.2 km s⁻¹. The molecular constants for the $^{13}$CO J = 2 − 1 transition, which are taken from CDMS, are $Q(20\,\text{K}) = 15.8$, $g = 10$, $E_u = 15.9$ K, and $A_{ul} = 6.08 \times 10^{-7}$ s⁻¹ (Müller et al. 2001, 2005). (Note that the upper state degeneracy differs from the LAMDA database value (Schöier et al. 2005)

listed in Öberg et al. (2021) by a factor of 2 because the CDMS database accounts for hyperfine splitting in the partition function). Equation (1) yields a $^{13}$CO column density upper limit of ~6 × $10^{13}$ cm⁻². Assuming that the $^{12}$C/$^{13}$C ratio has an ISM-like value of 69 (e.g., Wilson 1999), the $^{12}$CO:H₂ ratio is $10^{-4}$, and that the gas is predominantly molecular hydrogen, we again integrate over the area of the ellipse used to extract the flux of the $^{12}$CO extended structures and estimate a nominal gas mass upper bound of $10^{-4}$ M☉. This could still be an underestimate of the gas mass of the extended structures if the $^{12}$CO:H₂ ratio is substantially less than $10^{-4}$ due to freeze-out, photodissociation, and/or chemical depletion. Models comparing GM Aur's CO column densities to the disk gas mass (as constrained by HD) suggest that, even in the warm molecular layer, where neither freeze-out nor photodissociation should significantly affect CO abundances, CO is still depleted by 1–2 orders of magnitude relative to ISM levels. Nevertheless, even with the uncertainty in the CO:H₂ ratio, the total mass contained in the extended structures appears to be considerably less than the gas mass of the Keplerian disk, which has been estimated to be ~0.02–0.2 M☉ (e.g., McClure et al. 2016; Schwarz et al. 2021).

# 4. Discussion

## 4.1. Possible Origins of GM Aur's Complex Gas Structures

### 4.1.1. Accretion of Remnant Envelope or Cloud Material

In Section 3.1.3, we argued that GM Aur's extended gas kinematics are consistent with infalling motion. Morphologically, GM Aur's southwest tail and diffuse northern emission structures resemble scaled-down versions of the infalling gas streams that have been detected in molecular emission toward some (younger) Class 0 and I sources (e.g., Yen et al. 2014, 2017; Alves et al. 2020; Pineda et al. 2020). Hydrodynamical simulations also indicate that infalling material can induce the formation of spiral-like structures in or around disks (e.g., Stamatellos et al. 2012; Bae et al. 2015; Lesur et al. 2015; Kuffmeier et al. 2017). Thus, GM Aur's complex CO structures could result from ongoing accretion of material from the cloud or a remnant envelope.

Notably, GM Aur underwent a “burst” in 2018, in which the stellar accretion rate changed by a factor of 3.5 within a couple of weeks (Robinson & Espaillat 2019). The accretion rate change
<figure><img src="image_18.png" /><figcaption>**Figure 16.** $^{13}$CO J = 2 – 1 channel maps annotated with the locations of $^{12}$CO features for which counterparts are tentatively identified in $^{13}$CO. The blue contours mark the 4 and 8σ emission levels, where σ = 1 mJy beam$^{-1}$. The Keplerian isovelocity contours are drawn in gray, with solid curves marking the front side of the disk and dashed curves marking the back side. The synthesized beam size, 0.′39 × 0.′35 (−4.′4), is represented with a white ellipse in the lower left corner of each panel.</figcaption></figure>

coincided with a change in the disk’s mid-infrared emission,
leading Espaillat et al. (2019) to hypothesize that the increased
stellar accretion rate resulted from an increase in the inner disk’s
surface density. At earlier evolutionary stages, the high stellar
accretion rates and powerful outbursts observed in FU Ori-like
systems have often been attributed to instabilities triggered by a
build-up of material in the disk due to envelope infall (e.g.,
Hartmann & Kenyon 1996; Zhu et al. 2009, 2010; Bae et al.
2014). Perhaps in an analogous fashion, infall from the large-
scale gas structures around GM Aur may be responsible for its

hypothesized inner disk surface density build-up and “burst”-like
behavior. Assessing the plausibility of such a link in the GM Aur
system will require spatially resolved observations of other Class
II systems exhibiting similar burst-like behavior.

A number of pre-main-sequence stars have long been known
to be surrounded by nebulosities detected in optical and infrared
images (e.g., Joy 1945; Herbig 1960). The nature of the physical
relationship between these stars and the nearby nebulosities,
though, had not always been clear. However, increasingly
powerful millimeter interferometers and high-contrast optical/
<figure><img src="image_19.png" /><figcaption>**Figure 17.** An integrated intensity map of $^{12}$CO J = 2–1 emission, with the Keplerian disk masked. Edge artifacts from the mask are visible in the central region of the image. An arcsinh color stretch is used to make faint structures more visible. The pink star marks the location of GM Aur. The dotted white ellipse denotes the area over which the flux of the extended structures is measured to estimate their gas mass. The synthesized beam is shown as an ellipse in the lower left corner.</figcaption></figure>

near-infrared imagers have been able to resolve large-scale spirals, streams, and tails that can be traced down to disk scales (e.g., Grady et al. 1999; Liu et al. 2016; Garufi et al. 2020; Huang et al. 2020b). It is interesting to note that GM Aur, AB Aur, and SU Aur, the three most extensively observed Class II systems in the L1517 cloud, all exhibit spiral structures in CO and/or scattered light, strong deviations from Keplerian kinematics, and connections (at least in projection) to large-scale arc or tail-like nebulosities detected in scattered light (e.g., Grady et al. 1999; Schneider et al. 2003; Tang et al. 2012; de Leon et al. 2015; Akiyama et al. 2019; Boccaletti et al. 2020; Ginski et al. 2021). Figure 18 shows the locations of these systems in the L1517 cloud and images of the extended structures detected in association with their disks. Based on their kinematics and morphology, the AB Aur and SU Aur systems have often been cited as likely examples of disks undergoing late infall, where “late” generally refers to Class II disks more than 1 Myr old (e.g., Grady et al. 1999; Tang et al. 2012; Akiyama et al. 2019; Ginski et al. 2021). The similar properties of the extended gas structures of GM Aur, AB Aur, and SU Aur, in spite of their disparate stellar and millimeter continuum disk properties, lend further weight to the idea that the gas structures result from processes external to rather than within the disk. On a more speculative note, AB Aur and SU Aur's extended structures appear more prominent in CO and scattered light compared to those of GM Aur (see aforementioned references). In projection, AB Aur and SU Aur are positioned in higher-density regions of the L1517 cloud than GM Aur is (see Figure 18). Observations of additional

disks in L1517 could help to determine what relationship may exist between the large-scale distribution of material in the cloud and the local disk environment.

Late infall has been hypothesized to provide access to a significant mass reservoir for planet formation, to misalign disks and therefore the orbits of planets that subsequently form, and to enhance stellar accretion rates, potentially leading to outbursts that heat up the disk (e.g., Manara et al. 2018; Dullemond et al. 2019; Kuffmeier et al. 2020). Systematic surveys of disks in gas tracers will be fundamental for determining the prevalence of late infall and therefore the magnitude of its potential impact on planet demographics. While $^{12}$CO and scattered-light observations can be used to establish whether the kinematics and morphology of the large-scale gas structures are indicative of late infall, they are not well-suited by themselves for assessing whether the mass or angular momentum transfer is sufficient to yield some of the dramatic effects predicted by models. Thus, a critical step in assessing the role that late infall might play in planet formation is to search for suitable spatially resolved mass and temperature tracers via observations of other molecular lines. Potentially useful targets include CS or CN, since they have been shown sometimes to trace extended structures around embedded disks (e.g., Garufi et al. 2021). While no obvious large-scale structures are identified around GM Aur in CS $J=2-1$ emission by Le Gal et al. (2021) or in CN $N=1-0$ by Bergner et al. (2021), these structures might be recoverable in brighter transitions.

### 4.1.2. Gravitational Instability

Gravitational instabilities (GI), potentially precipitated by infalling material, may drive spiral arm formation as well as the formation of tail-like structures through gas clumps ejected when the disk fragments (e.g., Boss 1997; Gammie 2001; Vorobyov & Basu 2005; Kratter et al. 2008; Harsono et al. 2011; Vorobyov et al. 2020). It is ambiguous from Toomre Q parameter estimates whether GM Aur is in fact unstable. From modeling HD emission, McClure et al. (2016) report that $Q \sim 1.3$ at a radius of 300 au, which would make the disk marginally gravitationally unstable. From modeling the same HD observations as well as new spatially resolved CO isotopologue observations from MAPS, Schwarz et al. (2021) find that the disk may be gravitationally unstable between $\sim 70–100$ au (where $Q < 1.7$ is used as the metric), but stable at all other radii. This would suggest that GM Aur's spiral arms, which are located hundreds of astronomical units from the star, are not a consequence of GI. However, Schwarz et al. (2021) note that their surface densities beyond 100 au may be underestimated (and therefore the Toomre Q parameter overestimated) because HD emission comes from warmer gas. Using multifrequency millimeter continuum observations and assuming a gas-to-dust ratio of 100, Sierra et al. (2021) find that the disk is gravitationally stable between 25 and 125 au, with Q values that are systematically higher than those estimated by Schwarz et al. (2021). Because the millimeter continuum emission is compact compared to the gas, the Q parameter estimates from Sierra et al. (2021) are not directly relevant to assessing the disk stability at the location of the spiral arms. However, the discrepant values that have been estimated for GM Aur highlight the challenges of assessing disk stability.

In principle, disk kinematics offer an independent avenue to assess whether a disk is unstable. Hall et al. (2020) use hydrodynamical simulations to demonstrate that GI should
<figure><img src="image_20.png" /><figcaption>**Figure 18.** Left: The positions of AB Aur, SU Aur, and GM Aur marked on a Herschel SPIRE (Griffin et al. 2010; Pilbratt et al. 2010) 500 µm map tracing dust emission in the L1517 cloud. The version of the map retrieved from the Herschel Science Archive is the Level 2.5 extended data processing product, a combination of ObsID 1342204843 and 1342204844 from the Herschel Gould Belt survey (André et al. 2010). Top right: ¹²CO J = 3 – 2 image of SU Aur (Ginski et al. 2021). Middle right: R-band image of AB Aur (Kalas 1997; Grady et al. 1999). Bottom right: ¹²CO J = 2 – 1 image of GM Aur from this work.</figcaption></figure>

induce a characteristic non-Keplerian “wiggle” in CO channel maps. Given the overall kinematic and geometric complexity of GM Aur’s CO emission, it is not obvious whether such a “wiggle” is present. We note that, in any case, disk instabilities alone would not completely account for GM Aur’s morphology, particularly the diffuse northern structures.

### 4.1.3. Other Explanations

Simulations have shown that close stellar encounters can create spiral and tail structures in disks (e.g., Clarke & Pringle 1993; Pfalzner 2003; Cuello et al. 2019; Vorobyov et al. 2020). Gas structures reminiscent of the GM Aur system have also been detected in several systems with stars separated by projected distances of a few hundred au, including RW Aur A and B (Rodríguez et al. 2018), AS 205 N and S (Kurtovic et al. 2018), and UX Tau A and C (Zapata et al. 2020).

Unlike these other systems, though, there are no strong candidates for objects that have undergone close encounters with GM Aur. Simulations indicate that, following an encounter, spiral

arms should only persist in the disk for several thousand years (e.g., Cuello et al. 2019). The one-dimensional velocity dispersion of the Taurus molecular cloud complex has been estimated to be ~3 km s⁻¹ (Galli et al. 2018). Supposing that GM Aur experienced an encounter 5000 years ago with a star moving away at a relatively large velocity of 10 km s⁻¹, their present-day separation would be ~11,000 au. Of the Taurus members identified in Luhman (2018) and Esplin & Luhman (2019) from the Gaia mission (Gaia Collaboration et al. 2016, 2018) and follow-up observations, the closest source to GM Aur is 2MASS J04552333+3027366, a brown dwarf located 6.2 arcmin (59,300 au) northeast of GM Aur. Another source with a similar parallax, Gaia DR2 156913778302440704, is located closer in projection, 45″ (7200 au) northeast of GM Aur (Gaia Collaboration et al. 2016, 2018). However, Gaia DR2 156913778302440704 is not classified as a member of Taurus in Luhman (2018), because its Gaia-measured proper motion (μα = 10.5 ± 0.2 mas yr⁻¹, μδ = −32.8 ± 0.1 mas yr⁻¹) differs significantly from the median proper motion measured for the L1517 region (μα = 4.7 ± 0.8 mas yr⁻¹, μδ = −24.5 ± 1.7 mas yr⁻¹). Luhman (2018) estimated
that, in regions where $A_J < 1$, their Taurus catalog is complete for spectral types earlier than M6–M7, which corresponds to a magnitude of $G \sim 19$. This limit is valid in the vicinity of GM Aur, for which $A_J = 0.14$ (Luhman et al. 2017). Luhman (2018) estimated that, at an age of several million years (i.e., roughly that of the Taurus star-forming region), spectral types of M6–M7 correspond to stellar masses of roughly $0.1 M_{\odot}$. Meanwhile, Ziegler et al. (2018) determined that Gaia can recover nearly all secondary companions with $G < 21$ and separations greater than $3''$. In addition, the SEEDS high-contrast imaging survey, with a field of view of $\sim 400$ au, rules out companions more massive than $\sim 2.5 M_{\odot}$ beyond 50 au from GM Aur (Uyama et al. 2017). Together, these studies appear to rule out the presence of stellar mass objects in GM Aur's vicinity. A more qualitative argument against GM Aur having undergone a recent stellar encounter is that stellar encounters are expected to lead to tidal truncation of disks (e.g., Breslau et al. 2014), but the GM Aur disk has one of the largest known radial extents among Class II disks in both molecular and millimeter continuum emission (e.g., Koerner et al. 1993; Tripathi et al. 2017).

Stellar or planetary companions can also excite spiral density waves (e.g., Lin & Papaloizou 1986; Tanaka et al. 2002). Stellar companions appear unlikely to be the cause of GM Aur's spiral structures, for the same reasons that stellar encounters are not deemed likely to be responsible. Meanwhile, planetary-mass companions are expected to result in small and localized perturbations to the disk gas kinematics, rather than the large-scale non-Keplerian motion observed toward GM Aur (e.g., Pérez et al. 2018; Pinte et al. 2018; Teague et al. 2019). At any rate, a planetary-mass companion would not account for either the southwest tail or diffuse northern emission.

## 4.2. Exploring the Chemical Implications of GM Aur's Complex Gas Structures

Given that GM Aur was observed as part of the MAPS disk chemistry survey, it is worthwhile to consider how the processes that might be responsible for GM Aur's extended CO structures could also influence the disk chemistry, and how best to test these ideas. Chemically speaking, the outer disk of GM Aur probed at millimeter wavelengths is notable for both its faint C₂H and bright H₂CO emission (adjusted for distance) relative to similarly large, massive disks (Guilloteau et al. 2016; Bergner et al. 2019; Pegues et al. 2020; Guzmán et al. 2021; Öberg et al. 2021).

At first glance, GM Aur's properties suggest that it should be hospitable to C₂H production. The factors thought to be most critical for C₂H production in disks are high volatile C/O ratios and strong UV radiation (e.g., Bergin et al. 2016; Bosman et al. 2021b). While the C/O ratio has not been measured directly for GM Aur, its gas-phase CO has been estimated to be depleted up to two orders of magnitude (e.g., McClure et al. 2016; Schwarz et al. 2021; Zhang et al. 2021). High volatile CO depletion is often thought to signify elevated C/O ratios (e.g., Kama et al. 2016; Bergin et al. 2016; Miotello et al. 2019). Meanwhile, GM Aur's UV radiation field (as quantified by $L_{Lyα}/L_{UV}$ and $L_{FUV}/L_{UV}$) is comparable to other Class II disks known to have bright C₂H emission (e.g., Arulanantham et al. 2020). One possible explanation for GM Aur's more modest C₂H emission is that, if the disk is indeed still accreting cloud or envelope material, then enough oxygen from either CO or photodesorbing water ice might be introduced into the surface layers to reduce hydrocarbon production without

resetting the chemistry in the deeper layers containing the bulk of the disk mass. Astrochemical modeling of infalling material onto a Class II disk will be useful for quantifying whether this mechanism accounts for the observed hydrocarbon behavior.

Although a large disk mass could be a straightforward explanation for GM Aur's bright H₂CO emission, it is still notable insofar as some of GM Aur's other molecular lines are faint compared to other sources with comparable or lower masses (Pegues et al. 2020; Guzmán et al. 2021; Öberg et al. 2021). While the number of measurements is small, observations of Class I disks suggest that their H₂CO column densities tend to be higher than in Class II disks (van't Hoff et al. 2020; Garufi et al. 2021). If GM Aur is accreting primordial gas, some of its chemical characteristics may bear closer resemblance to a young, embedded disk than to other Class II disks. A possible reason for GM Aur's bright H₂CO emission could be shock heating induced by infall and/or disk instabilities, which in turn may lead to enhanced H₂CO desorption from grains as well as increased gas-phase production (e.g., Ilee et al. 2011; Evans et al. 2019). An excitation temperature measurement of H₂CO in this disk would be useful to ascertain whether shock heating is indeed occurring.

While GM Aur has been a fairly popular target for disk chemistry observations, interferometric studies of molecular emission have otherwise preferentially targeted sources exhibiting relatively clean signatures of Keplerian rotation, in order to avoid confusion with cloud material (e.g., Piétu et al. 2007; Qi et al. 2008; Öberg et al. 2010; Bergin et al. 2016). By definition, choosing targets based on their Keplerian signatures has limited the extent to which Class II disks undergoing infall have been characterized. A broader sample needs to be observed to investigate whether there are indeed systematic chemical differences between sources showing evidence of interacting with their surroundings versus those that are more isolated.

## 5. Summary and Conclusions

As part of the MAPS Large Program, we observed CO toward the GM Aur disk at high spatial resolution and sensitivity. Our major findings are as follows:

1. In addition to a Keplerian disk visible up to a radius of ~550 au, ¹²CO emission reveals spiral arms out to a radius of ~1200 au, a tail extending up to ~1800 au southwest of GM Aur, and diffuse northern structures up to ~1900 au from the star. These large-scale non-axisymmetric gas structures contrast sharply with the nearly axisymmetric, multi-ringed millimeter continuum, which is detected only up to ~250 au from the star.

2. Portions of the diffuse northern structures and southwest tail have counterparts in earlier HST scattered-light observations tracing the small-grain distribution.

3. Based on the kinematics and morphology of GM Aur's extended gas structures, we hypothesize that the GM Aur disk is experiencing late infall, either from a remnant protostellar envelope or from surrounding cloud material.

4. The extended ¹²CO structures are only tentatively detected in ¹³CO. The absence of strong extended emission in other molecular tracers underscores the utility of observing bright ¹²CO lines to probe interactions between disks and their surroundings.
5. Given GM Aur's relatively weak C₂H and strong H₂CO emission in comparison with other massive Class II disks that have been observed, we speculate that late infall could serve as a "fountain of youth" that partially resets disk chemistry. Molecular line observations of other systems exhibiting signatures of late infall will be necessary in order to establish if and how environmental interactions affect disk and protoplanet compositions.

While recent ALMA surveys have taken advantage of the array's unparalleled spatial resolution and sensitivity to probe the inner regions of protoplanetary disks, serendipitous detections of extended gas structures such as spirals and tails have shown that much is left to be learned about the structure and evolution of disks on larger scales as well. Ascertaining how disks interact with their surroundings through deep, systematic ¹²CO observations will be key for determining how protoplanets ultimately acquire their starting material.

This paper makes use of ALMA data ADS/JAO.ALMA#2018.1.01055.L, ADS/JAO.ALMA#2017.1.01151.S, and ADS/JAO.ALMA#2018.1.01230.S. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada) and NSC and ASIAA (Taiwan), 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 research is based on observations made with the NASA/ESA Hubble Space Telescope obtained from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-26555. These observations are associated with program(s) HST-GO-6223. Herschel is an ESA space observatory with science instruments provided by European-led Principal Investigator consortia and with important participation from NASA. SPIRE has been developed by a consortium of institutes led by Cardiff University (UK) and including Univ. Lethbridge (Canada); NAOC (China); CEA, LAM (France); IFSI, Univ. Padua (Italy); IAC (Spain); Stockholm Observatory (Sweden); Imperial College London, RAL, UCL-MSSL, UKATC, Univ. Sussex (UK); and Caltech, JPL, NHSC, Univ. Colorado (USA). This development has been supported by national funding agencies: CSA (Canada); NAOC (China); CEA, CNES, CNRS (France); ASI (Italy); MCINN (Spain); SNSB (Sweden); STFC, UKSA (UK); and NASA (USA). This research has made use of data from the Herschel Gould Belt survey (HGBS) project (http://goudbelt-herschel.cea.fr). The HGBS is a Herschel Key Programme jointly carried out by SPIRE Specialist Astronomy Group 3 (SAG 3), scientists of several institutes in the PACS Consortium (CEA Saclay, INAF-IFSI Rome and INAF-Arcetri, KU Leuven, MPIA Heidelberg), and scientists of the Herschel Science Center (HSC). This research has made use of NASA's Astrophysics Data System and the VizieR catalog access tool, CDS, Strasbourg, France.

Support for J.H. was provided by NASA through the NASA Hubble Fellowship grant No. HST-HF2-51460.001-A awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS5-26555. J.H. and S.M.A. acknowledge funding support from the National Aeronautics and Space Administration under grant No. 17-XRP17 2-0012

issued through the Exoplanets Research Program. E.A.B and A.D.B. acknowledge support from NSF AAG grant No. 1907653. K.I.Ö. acknowledges support from the Simons Foundation (SCOL #321183) and an NSF AAG grant (No. 1907653). R.T. and F.L. acknowledge support from the Smithsonian Institution as Submillimeter Array (SMA) Fellows. C.J.L. acknowledges funding from the National Science Foundation Graduate Research Fellowship under grant No. DGE1745303. P.K. was supported by NNX15AD95G, part of NASA's Nexus for Exoplanet System Science (NExSS) research coordination network sponsored by NASA's Science Mission Directorate. Y.A. acknowledges support by NAOJ ALMA Scientific Research grant code 2019-13B and Grant-in-Aid for Scientific Research (S) 18H05222. J.B. acknowledges support by NASA through the NASA Hubble Fellowship grant #HST-HF2-51427.001-A awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Incorporated, under NASA contract NAS5-26555. J.B.B. acknowledges support from NASA through the NASA Hubble Fellowship grant #HST-HF2-51429.001-A, awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS5-26555. A.S.B acknowledges the studentship funded by the Science and Technology Facilities Council of the United Kingdom (STFC). J.K.C. acknowledges support from the National Science Foundation Graduate Research Fellowship under grant No. DGE 1256260 and the National Aeronautics and Space Administration FINESST grant, under grant No. 80NSSC19K1534. G.C. is supported by NAOJ ALMA Scientific Research grant code 2019-13B. L.I.C. gratefully acknowledges support from the David and Lucile Packard Foundation and Johnson & Johnson's WiSTEM2D Program. I.C. was supported by NASA through the NASA Hubble Fellowship grant HST-HF2-51405.001-A awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS5-26555. J.D.I. acknowledges support from the STFC under ST/T000287/1. R.L.G. acknowledges support from a CNES fellowship grant. V.V.G. acknowledges support from FONDECYT Iniciación 11180904 and ANID project Basal AFB-170002. F.M. acknowledges support from ANR of France under contract ANR-16-CE31-0013 (Planet-Forming-Disks) and ANR-15-IDEX-02 (through CDP "Origins of Life"). H.N. acknowledges support by NAOJ ALMA Scientific Research grant code 2018-10B and Grant-in-Aid for Scientific Research 18H05441. K.R.S. acknowledges the support of NASA through the NASA Hubble Fellowship grant No. HST-HF2-51419.001 awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS5-26555. T.T. is supported by JSPS KAKENHI grant Nos. JP17K14244 and JP20K04017. M.L.R.H. acknowledges support from the Michigan Society of Fellows. Y.Y. is supported by IGPEES, WINGS Program, the University of Tokyo. C.W. acknowledges financial support from the University of Leeds, STFC, and UKRI (grant Nos. ST/R000549/1, ST/T000287/1, MR/T040726/1). K.Z. acknowledges the support of the Office of the Vice Chancellor for Research and Graduate Education at the University of Wisconsin-Madison with funding from the Wisconsin Alumni Research Foundation, and the support of NASA through
Hubble Fellowship grant HST-HF2-51401.001. awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS5-26555. We thank Robert Gutermuth, Glenn Schneider, Jonathan Williams, Alvaro Hacar, and Megan Ansdell for useful discussions, as well as the anonymous referee for comments improving the manuscript.

Facilities: ALMA, HST (WFPC2), Herschel.

Software: `analysisUtils` (https://casaguides.nrao.edu/index.php/Analysis_Utilities), `AstroPy` (Astropy Collaboration et al. 2013), `CASA` (McMullin et al. 2007), `gofish` (Teague 2019),

`keplerian_mask` (Teague 2020), `VESPRAP` (Hunt 2007), `NumPy` (Harris et al. 2020), `numba` (Lam et al. 2015), `scikit-image` (van der Walt et al. 2014), `SciPy` (Virtanen et al. 2020).

## Appendix A

### $^{12}$CO Channel Maps

Figure 19 presents $^{12}$CO J = 2 – 1 channel maps toward GM Aur without annotations.

Figure 20 shows the best-fit curve to the northwest ridge overlaid on the $^{12}$CO channel maps.
<figure><img src="image_21.png" /><figcaption>**Figure 19.** Channel maps of the $^{12}$CO J = 2–1 emission toward the GM Aur protoplanetary disk. The synthesized beam is shown in the lower left corner and the LSRK velocity (km s$^{-1}$) is labeled in the upper right corner of every panel. Offsets from the disk center (arcseconds) are marked on the axes in the lower left corner of the figure. An arcsinh color stretch is used to highlight faint emission. North is up and east is to the left.</figcaption></figure>
<figure><img src="image_22.png" /><figcaption>Figure 19. (Continued.)</figcaption></figure>
<figure><img src="image_23.png" /><figcaption>**Figure 20.** Channel maps of the ¹²CO J = 2 – 1 emission toward GM Aur with the best-fit curve to the northwest ridge overlaid in pink. The light blue contours show the 5σ, 15σ, 25σ, and 35σ emission levels, where σ = 1 mJy beam⁻¹. The dark blue curves show the isovelocity contours of the Keplerian disk, with solid curves denoting the front of the disk and dashed curves denoting the back side.</figcaption></figure>

## Appendix B
Computing Isovelocity Contours for GM Aur's Keplerian Disk

To compute the isovelocity contours approximately corresponding to the portion of ¹²CO J = 2 – 1 emission tracing GM Aur's Keplerian disk, we adopt the following emitting height z (au) measured by Law et al. (2021b):

$$z(r) = 61.2 \text{ au} \left(\frac{r}{159 \text{ au}}\right)^{1.066} \exp\left[-\left(\frac{r}{599.0 \text{ au}}\right)^{4.988}\right], \quad (B1)$$

where r is the disk radius (au) in cylindrical coordinates. The azimuthal velocity of the gas $v_φ$ in the frame of the disk is
calculated with the expression

$$ \frac{v_{\varphi}^2}{r} = \frac{GM_*r}{(r^2 + z^2)^{\frac{3}{2}}} + \frac{1}{\rho(r, z)} \frac{\partial P(r, z)}{\partial r}, \quad (B2) $$

where $\rho$ is the gas density and $P$ is the gas pressure. The two terms on the right-hand side account for the effects of stellar gravity and the gas pressure gradient, respectively. A stellar mass of $M_* = 1.1 M_\odot$ is adopted (Macías et al. 2018). The gas pressure gradient term is often dropped in velocity calculations because purely Keplerian rotation can usually reproduce disk molecular emission well (e.g., Dutrey et al. 1994; Simon et al. 2000; Czekala et al. 2015). However, deviations from Keplerian motion due to the pressure gradient become more prominent at large disk radii and heights, so the effect of the pressure gradient should be considered in the analysis of well-resolved $^{12}$CO observations of large disks (e.g., Rosenfeld et al. 2013; Dullemond et al. 2020). For brevity, even though the disk velocities are slightly sub-Keplerian, we refer to the $^{12}$CO isovelocity contours calculated from Equation (B2) as the isovelocity contours of the Keplerian disk to make a distinction from the gas structures tracing the spiral arms, southwest tail, and diffuse northern structures. Equation (B2) neglects self-gravity, which may not be accurate for the GM Aur disk because some estimates place its mass at or close to the gravitationally unstable limit (e.g., McClure et al. 2016; Schwarz et al. 2021). While we can produce isovelocity contours that match reasonably well with GM Aur's emission without considering self-gravity, it should be kept in mind that the isovelocity contour calculations are intended to provide a visual aid rather than a model that captures the full physical complexity of the system.

The pressure is given by $P(r, z) = \rho(r, z)c_s^2$, where $c_s = \sqrt{\frac{k_B T(r)}{\mu m_H}}$ is the isothermal sound speed. We set $\mu = 2.37$. Based on the GM Aur $^{12}$CO $J=2-1$ brightness temperatures measured in Law et al. (2021b), the temperature structure is approximated as

$$ T(r) = 52 K \left( \frac{r}{100 \text{ au}} \right)^{-0.608}. \quad (B3) $$

The gas density profile is taken from the GM Aur model presented in Zhang et al. (2021). The disk gas surface density is

$$ \Sigma_g(r) = 9.4 \text{ g cm}^{-2} \left(\frac{r}{176 \text{ au}}\right)^{-1} \exp\left[-\left(\frac{r}{176 \text{ au}}\right)\right]. \quad (B4) $$

GM Aur's gas density is approximated as

$$ \rho(r, z) = \frac{\Sigma_g(r)}{\sqrt{2\pi} H(r)} \exp\left(-\frac{z^2}{2H(r)^2}\right), \quad (B5) $$

where

$$ H(r) = 7.5 \text{ au} \left( \frac{r}{100 \text{ au}} \right)^{1.35} \quad (B6) $$

is the gas pressure scale height.

The azimuthal velocity is related to the line-of-sight velocity $v_{los}$ by

$$ v_{los} = v_{\text{sys}} + v_{\varphi} \cos \varphi \sin i, \quad (B7) $$

where $i = 53^\circ 21$ is the disk inclination and $v_{sys} = 5.61 \text{ km s}^{-1}$ is the systemic velocity (Huang et al. 2020a). Using the gofish Python package (Teague 2019), the position of each pixel in the channel maps in the frame of the observer is mapped to the $r, \varphi, z$ coordinates of the upper and lower surfaces of the $^{12}$CO emission in the frame of the disk in order to compute the line of sight velocity at each position.

The $^{13}$CO $J=2-1$ isovelocity contours are computed similarly based on the emission height and brightness temperature measurements from Law et al. (2021b):

$$ z(r) = 18.0 \text{ au} \left(\frac{r}{159 \text{ au}}\right)^{4.539} \exp\left[-\left(\frac{r}{237.9 \text{ au}}\right)^{4.989}\right] \quad (B8) $$

and

$$ T(r) = 22 K \left( \frac{r}{100 \text{ au}} \right)^{-0.260}, \quad (B9) $$

respectively.

## Appendix C: Constructing a Keplerian Mask for GM Aur's $^{12}$CO Emission

Because GM Aur's extended $^{12}$CO structures overlap in projection with the Keplerian disk but have different line-of-sight velocities, we use a Keplerian mask to isolate emission originating from the extended structures. The keplerian_mask Python package (Teague 2020) can calculate Boolean masks based on the expected emitting region of a Keplerian disk for a given stellar mass, disk inclination and position angle, emitting surface height $z(r)$, and disk radius. We used the same stellar mass, disk orientation, and $z(r)$ used to draw the isovelocity contours described in Section B. The maximum radius was set to $r=550$ au based on visual inspection. The code also assumes that the intrinsic linewidth has the form

$$ \Delta V(r) = \Delta V_0 \times \left(\frac{r}{1''}\right)^{\Delta V_q} \quad (C1) $$

We used the code's default values of $\Delta V_0 = 300 \text{ m s}^{-1}$ and $\Delta V_q = -0.5$ because they yielded a mask that reasonably encompassed the emission from the Keplerian disk, but the adopted linewidth profile should not necessarily be taken as a close representation of the true disk conditions. Finally, the mask was convolved with a two-dimensional Gaussian with the same size as the synthesized beam in order to account for emission broadening
<figure><img src="image_24.png" /><figcaption>**Figure 21.** The ¹²CO Keplerian mask (outlined in pink) used to produce the integrated intensity map shown in Figure 17, overlaid onto the ¹²CO channel maps. The LSRK velocity (km s⁻¹) is noted in the upper right corner of each panel, and the synthesized beam is shown in the lower left corner.</figcaption></figure>

due to resolution limits. Figure 21 shows the mask overlaid on the ¹²CO channel maps.

## ORCID iDs

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

Edwin A. Bergin https://orcid.org/0000-0003-4179-6394

Karin I. Öberg https://orcid.org/0000-0001-8798-1347

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

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

Charles J. Law https://orcid.org/0000-0003-1413-1776

Paul Kalas https://orcid.org/0000-0002-6221-5360

Yuri Aikawa https://orcid.org/0000-0003-3283-6884

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

Jennifer B. Bergner https://orcid.org/0000-0002-8716-0482

Alice S. Booth https://orcid.org/0000-0003-2014-2121

Arthur D. Bosman https://orcid.org/0000-0003-4001-3589

Jenny K. Calahan https://orcid.org/0000-0002-0150-0125

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

L. Ilsedore Cleeves https://orcid.org/0000-0003-2076-8001

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

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

Romane Le Gal https://orcid.org/0000-0003-1837-3772

Viviana V. Guzmán https://orcid.org/0000-0003-4784-3040

Feng Long https://orcid.org/0000-0002-7607-719X

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

François Ménard https://orcid.org/0000-0002-1637-7393

Hideko Nomura https://orcid.org/0000-0002-7058-7682

Chunhua Qi https://orcid.org/0000-0001-8642-1786
Kamber R. Schwarz https://orcid.org/0000-0002-6429-9457

Takashi Tsukagoshi https://orcid.org/0000-0002-6034-2892

Merel L. R. van 't Hoff https://orcid.org/0000-0002-2555-9869

Catherine Walsh https://orcid.org/0000-0001-6078-786X

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

Yoshihide Yamato https://orcid.org/0000-0003-4099-6941

Ke Zhang https://orcid.org/0000-0002-0661-7517

## References

Akiyama, E., Vorobyov, E. I., Baobabu Liu, H., et al. 2019, AJ, 157, 165

ALMA Partnership, Brogan, C. L., Pérez, L. M., et al. 2015, ApJL, 808, L3

Alves, F. O., Cleeves, L. I., Girart, J. M., et al. 2020, ApJL, 904, L6

André, P., Men'shchikov, A., Bontemps, S., et al. 2010, A&A, 518, L102

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

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

Arulanantham, N., France, K., Cazzoletti, P., et al. 2020, AJ, 159, 168

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

Bae, J., Hartmann, L., & Zhu, Z. 2015, ApJ, 805, 15

Bae, J., Hartmann, L., Zhu, Z., & Nelson, R. P. 2014, ApJ, 795, 61

Bae, J., Zhu, Z., & Hartmann, L. 2017, ApJ, 850, 201

Bailer-Jones, C. A. L., Rybizki, J., Fouesneau, M., Mantelet, G., & Andrae, R. 2018, AJ, 156, 58

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

Bergin, E., Calvet, N., Sitko, M. L., et al. 2004, ApJL, 614, L133

Bergin, E. A., Du, F., Cleeves, L. I., et al. 2016, ApJ, 831, 101

Bergner, J. B., Öberg, K. I., Guzmán, V. V., et al. 2021, ApJS, 257, 11

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

Boccaletti, A., Di Folco, E., Pantin, E., et al. 2020, A&A, 637, L5

Bosman, A. D., Alarcón, F., Bergin, E. A., et al. 2021b, ApJS, 257, 7

Bosman, A., Bergin, E. A., Loomis, R. A., et al. 2021a, ApJS, 257, 15

Boss, A. P. 1997, Sci, 276, 1836

Breslau, A., Steinhausen, M., Vincke, K., & Pfalzner, S. 2014, A&A, 565, A130

Calvet, N., D'Alessio, P., Watson, D. M., et al. 2005, ApJ, 630, L185

Christiaens, V., Casassus, S., Perez, S., van der Plas, G., & Ménard, F. 2014, ApJL, 785, L12

Clarke, C. J., & Pringle, J. E. 1993, MNRAS, 261, 190

Cornwell, T. J. 2008, ISTSP, 2, 793

Cuello, N., Dipierro, G., Mentiplay, D., et al. 2019, MNRAS, 483, 4114

Czekala, I., Loomis, R. A., Teague, R., et al. 2021, ApJS, 257, 2

Czekala, I., Andrews, S. M., Jensen, E. L. N., et al. 2015, ApJ, 806, 154

de Leon, J., Takami, M., Karr, J. L., et al. 2015, ApJL, 806, L10

Dipierro, G., Price, D., Laibe, G., et al. 2015, MNRAS, 453, L73

Dong, R., Fung, J., & Chiang, E. 2016, ApJ, 826, 75

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

Dullemond, C. P., Küffmeier, M., Goicovic, F., et al. 2019, A&A, 628, A20

Dutrey, A., Guilloteau, S., & Guelin, M. 1997, A&A, 317, L55

Dutrey, A., Guilloteau, S., Piétu, V., et al. 2008, A&A, 490, L15

Dutrey, A., Guilloteau, S., & Simon, M. 1994, A&A, 286, 149

Espaillat, C., D'Alessio, P., Hernández, J., et al. 2010, ApJ, 717, 441

Espaillat, C. C., Macías, E., Hernández, J., & Robinson, C. 2019, ApJL, 877, L34

Esplin, T. L., & Luhman, K. L. 2019, AJ, 158, 54

Evans, M. G., Hartquist, T. W., Caselli, P., et al. 2019, MNRAS, 483, 1266

Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2016, A&A, 595, A2

Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2018, A&A, 616, A1

Galli, P. A. B., Loinard, L., Ortiz-Léon, G. N., et al. 2018, ApJ, 859, 33

Gammie, C. F. 2001, ApJ, 553, 174

Garufi, A., Avenhaus, H., Pérez, S., et al. 2020, A&A, 633, A82

Garufi, A., Podio, L., Codella, C., et al. 2021, A&A, 645, A145

Ginski, C., Facchini, S., Huang, J., et al. 2021, ApJL, 908, L25

Goldsmith, P. F., & Langer, W. D. 1999, ApJ, 517, 209

Grady, C. A., Woodgate, B., Schweiger, F. C., et al. 1999, ApJL, 523, L151

Griffin, M. J., Abergel, A., Abreu, A., et al. 2010, A&A, 518, L3

Guilloteau, S., Reboussin, L., Dutrey, A., et al. 2016, A&A, 592, A124

Guzmán, V. V., Bergner, J. B., Law, C. J., et al. 2021, ApJS, 257, 6

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

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

Harsono, D., Alexander, R. D., & Levin, Y. 2011, MNRAS, 413, 423

Hartmann, L., & Kenyon, S. J. 1996, ARA&A, 34, 207

Herbig, G. H. 1960, ApJS, 4, 337

Hornbeck, J. B., Swearingen, J. R., Grady, C. A., et al. 2016, ApJ, 829, 65

Huang, J., Andrews, S. M., Dullemond, C. P., et al. 2018, ApJL, 869, L42

Huang, J., Andrews, S. M., Dullemond, C. P., et al. 2020a, ApJ, 891, 48

Huang, J., Andrews, S. M., Öberg, K. I., et al. 2020b, ApJ, 898, 140

Huang, J., Öberg, K. I., Qi, C., et al. 2017, ApJ, 835, 231

Hughes, A. M., Andrews, S. M., Espaillat, C., et al. 2009, ApJ, 698, 131

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

Ilee, J. D., Boley, A. C., Caselli, P., et al. 2011, MNRAS, 417, 2950

Jorsater, S., & van Moorsel, G. A. 1995, AJ, 110, 2037

Joy, A. H. 1945, ApJ, 102, 168

Kalas, P. 1997, PhD thesis, Univ. of Hawaii

Kama, M., Bruderer, S., van Dishoeck, E. F., et al. 2016, A&A, 592, A83

Kepley, A. A., Tsutsumi, T., Brogan, C. L., et al. 2020, PASP, 132, 024505

Koerner, D. W., Sargent, A. I., & Beckwith, S. V. W. 1993, Icar, 106, 2

Kratter, K. M., Matzner, C. D., & Krumholz, M. R. 2008, ApJ, 681, 375

Kuffmeier, M., Goicovic, F. G., & Dullemond, C. P. 2020, A&A, 633, A3

Kuffmeier, M., Haugbølle, T., & Nordlund, Å. 2017, ApJ, 846, 7

Kurtovic, N. T., Pérez, L. M., Benisty, M., et al. 2018, ApJL, 869, L44

Lam, S. K., Pitrou, A., & Seibert, S. 2015, in Proc. of the Second Workshop on the LLVM Compiler Infrastructure in HPC, LLVM '15, ed. H. Finkel (New York: Association for Computing Machinery), 7

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

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

Le Gal, R., Öberg, K. I., Teague, R., et al. 2021, ApJS, 257, 12

Lesur, G., Hennebelle, P., & Fromang, S. 2015, A&A, 582, L9

Lin, D. N. C., & Papaloizou, J. 1986, ApJ, 307, 395

Liu, H. B., Takami, M., Kudo, T., et al. 2016, SciA, 2, e1500875

Long, F., Pinilla, P., Herczeg, G. J., et al. 2018, ApJ, 869, 17

Luhman, K. L. 2018, AJ, 156, 271

Luhman, K. L., Mamajek, E. E., Shukla, S. J., & Loutrel, N. P. 2017, AJ, 153, 46

Macías, E., Espaillat, C. C., Ribas, Á., et al. 2018, ApJ, 865, 37

Manara, C. F., Morbidelli, A., & Guillot, T. 2019, A&A, 618, L3

McClure, M. K., Bergin, E. A., Cleeves, L. I., et al. 2016, ApJ, 831, 167

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

Miotello, A., Facchini, S., van Dishoeck, E. F., et al. 2019, A&A, 631, A69

Monet, D. G., Levine, S. E., Canzian, B., et al. 2003, AJ, 125, 984

Müller, H. S. P., Schlöder, F., Stutzki, J., & Winnewisser, G. 2005, JMoSt, 742, 215

Müller, H. S. P., Thorwirth, S., Roth, D. A., & Winnewisser, G. 2001, A&A, 370, L49

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

Öberg, K. I., Qi, C., Fogel, J. K. J., et al. 2010, ApJS, 720, 400

Pegues, J., Öberg, K. I., Bergner, J. B., et al. 2020, ApJ, 890, 142

Pérez, S., Casassus, S., & Benítez-Llambay, P. 2018, MNRAS, 480, L12

Pfalzner, S. 2003, ApJ, 592, 986

Piétu, V., Dutrey, A., & Guilloteau, S. 2007, A&A, 467, 163

Pilbratt, G. L., Riedinger, J. R., Passvogel, T., et al. 2010, A&A, 518, L1

Pineda, J. E., Segura-Cox, D., Caselli, P., et al. 2020, NatAs, 4, 1158

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

Qi, C., Kessler, J. E., Koerner, D. W., Sargent, A. I., & Blake, G. A. 2003, ApJ, 597, 986

Qi, C., Wilner, D. J., Aikawa, Y., Blake, G. A., & Hogerheijde, M. R. 2008, ApJ, 681, 1396

Robinson, C. E., & Espaillat, C. C. 2019, ApJ, 874, 129

Rodriguez, J. E., Loomis, R., Cabrit, S., et al. 2018, ApJ, 859, 150

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

Sargent, A. I., & Beckwith, S. 1987, ApJ, 323, 294

Schneider, G., Wood, K., Silverstone, M. D., et al. 2003, AJ, 125, 1467

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

Schwarz, K. R., Calahan, J. K., Zhang, K., et al. 2021, ApJS, 257, 20

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

Simon, M., Dutrey, A., & Guilloteau, S. 2000, ApJ, 545, 1034

Stamatellos, D., Whitworth, A. P., & Hubber, D. A. 2012, MNRAS, 427, 1182

Stapelfeldt, K., Burrows, C. J., Krist, J. E. & WFPC2 Science Team 1997, in IAU Symp. 182, Herbig-Haro Flows and the Birth of Stars, ed. B. Reipurth & C. Bertout (Dordrecht: Springer), 355
Strom, K. M., Strom, S. E., Edwards, S., Cabrit, S., & Skrutskie, M. F. 1989, AJ, 97, 1451

Swam, M. S., & Swade, D. A. 1999, in ASP Conf. Ser. 172, Astronomical Data Analysis Software and Systems VIII, ed. D. M. Mehringer, R. L. Plante, & D. A. Roberts (San Francisco, CA: ASP), 203

Tanaka, H., Takeuchi, T., & Ward, W. R. 2002, ApJ, 565, 1257

Tang, Y. W., Guilloteau, S., Piétu, V., et al. 2012, A&A, 547, A84

Teague, R. 2019, JOSS, 4, 1632

Teague, R. 2020, richteague/keplerian_mask: Initial Release v1.0, Zenodo, doi: 10.5281/zenodo.4321137

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

Thi, W.-F., van Zadelhoff, G.-J., & van Dishoeck, E. F. 2004, A&A, 425, 955

Tripathi, A., Andrews, S. M., Birnstiel, T., & Wilner, D. J. 2017, ApJ, 845, 44

Uyama, T., Currie, T., Christiaens, V., et al. 2020, ApJ, 900, 135

Uyama, T., Hashimoto, J., Kuzuhara, M., et al. 2017, AJ, 153, 106

van der Walt, S., Schönberger, J. L., Nunez-Iglesias, J., et al. 2014, PeerJ, 2, e453

van't Hoff, M. L. R., Harsono, D., Tobin, J. J., et al. 2020, ApJ, 901, 166

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

Vorobyov, E. I., & Basu, S. 2005, ApJL, 633, L137

Vorobyov, E. I., Skliarevskii, A. M., Elbakyan, V. G., et al. 2020, A&A, 635, A196

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

Wilson, T. L. 1999, RPPh, 62, 143

Yen, H.-W., Gu, P.-G., Hirano, N., et al. 2019, ApJ, 880, 69

Yen, H.-W., Takakuwa, S., Chu, Y.-H., et al. 2017, A&A, 608, A134

Yen, H.-W., Takakuwa, S., Ohashi, N., et al. 2014, ApJ, 793, 1

Zapata, L. A., Rodríguez, L. F., Fernández-López, M., et al. 2020, ApJ, 896, 132

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

Zhang, S., Zhu, Z., Huang, J., et al. 2018, ApJL, 869, L47

Zhu, Z., Hartmann, L., & Gammie, C. 2010, ApJ, 713, 1143

Zhu, Z., Hartmann, L., Gammie, C., & McKinney, J. C. 2009, ApJ, 701, 620

Ziegler, C., Law, N. M., Baranec, C., et al. 2018, AJ, 156, 259