Molecules with ALMA at Planet-forming Scales. XX. The Massive Disk Around GM Aurigae

KAMBER R. SCHWARZ<sup>1,2,*</sup> JENNY K. CALAHAN<sup>3</sup> KE ZHANG<sup>4,3,†</sup> FELIPE ALARCÓN<sup>1</sup>,<sup>3</sup> YURI AIKAWA<sup>1</sup>,<sup>5</sup>

SEAN M. ANDREWS<sup>1,6</sup> JAEHAN BAE<sup>1,7,8,*</sup> EDWIN A. BERGIN<sup>1,3</sup> ALICE S. BOOTH<sup>1,9,10</sup> ARTHUR D. BOSMAN<sup>1,3</sup>

GIANNI CATALDI<sup>1</sup>,<sup>11,5</sup> L. ILSEDORE CLEEVES<sup>1,12</sup> IAN CZEKALA<sup>1,13,14,15,16,17</sup> JANE HUANG<sup>1,3,*</sup> JOHN D. ILEE<sup>1,18</sup>

CHARLES J. LAW<sup>1,6</sup> ROMANE LE GAL<sup>1,6,19,20,21</sup> YAO LIU<sup>1,22</sup> FENG LONG<sup>1,6</sup> RYAN A. LOOMIS<sup>1,23</sup>

ENRIQUE MACÍAS<sup>1,24,25</sup> MELISSA MCCLURE<sup>1,9</sup> FRANÇOIS MÉNARD<sup>1,26</sup> KARIN I. ÖBERG<sup>1,6</sup> RICHARD TEAGUE<sup>1,6</sup>

EWINE VAN DISHOECK<sup>1,9</sup> CATHERINE WALSH<sup>1,10</sup> AND DAVID J. WILNER<sup>1,6</sup>

<sup>1</sup>Lunar and Planetary Laboratory, University of Arizona, 1629 E. University Blvd, Tucson, AZ 85721, USA

<sup>2</sup>Max-Planck-Institut für Astronomie, Königstuhl 17, 69117 Heidel- berg, Germany

<sup>3</sup>Department of Astronomy, University of Michigan, 323 West Hall, 1085 South University Avenue, Ann Arbor, MI 48109, USA

<sup>4</sup>Department of Astronomy, University of Wisconsin-Madison, 475 N Charter St, Madison, WI 53706

<sup>5</sup>Department of Astronomy, Graduate School of Science, The University of Tokyo, Tokyo 113-0033, Japan

<sup>6</sup>Center for Astrophysics | Harvard & Smithsonian, 60 Garden St., Cambridge, MA 02138, USA

<sup>7</sup>Earth and Planets Laboratory, Carnegie Institution for Science, 5241 Broad Branch Road NW, Washington, DC 20015, USA

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

<sup>9</sup>Leiden Observatory, Leiden University, 2300 RA Leiden, the Netherlands

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

<sup>11</sup>National Astronomical Observatory of Japan, Osawa 2-21-1, Mitaka, Tokyo 181-8588, Japan

<sup>12</sup>Department of Astronomy, University of Virginia, 530 McCormick Rd, Charlottesville, VA 22904

<sup>13</sup>Department of Astronomy and Astrophysics, 525 Davey Laboratory, The Pennsylvania State University, University Park, PA 16802, USA

<sup>14</sup>Center for Exoplanets and Habitable Worlds, 525 Davey Laboratory, The Pennsylvania State University, University Park, PA 16802, USA

<sup>15</sup>Center for Astrostatistics, 525 Davey Laboratory, The Pennsylvania State University, University Park, PA 16802, USA

<sup>16</sup>Institute for Computational & Data Sciences, The Pennsylvania State University, University Park, PA 16802, USA

<sup>17</sup>Department of Astronomy, 501 Campbell Hall, University of California, Berkeley, CA 94720-3411, USA*

<sup>18</sup>School of Physics & Astronomy, University of Leeds, Leeds LS2 9JT, UK

<sup>19</sup>IRAP, Université de Toulouse, CNRS, CNES, UT3, Toulouse, France

<sup>20</sup>IPAG, Université Grenoble Alpes, CNRS, IPAG, F-38000 Grenoble, France

<sup>21</sup>IRAM, 300 rue de la piscine, F-38406 Saint-Martin d'Hères, France

<sup>22</sup>Purple Mountain Observatory & Key Laboratory for Radio Astronomy, Chinese Academy of Sciences, Nanjing 210023, China

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

<sup>24</sup>Joint ALMA Observatory, Avenida Alonso de Córdova 3107, Vitacura, Santiago, Chile

<sup>25</sup>European Southern Observatory, Avenida Alonso de Córdova 3107, Vitacura, Santiago, Chile

<sup>26</sup>Univ. Grenoble Alpes, CNRS, IPAG, F-38000 Grenoble, France

(Received; Revised September 15, 2021; Accepted)

Submitted to AJ

ABSTRACT

Gas mass remains one of the most difficult protoplanetary disk properties to constrain. With much of the protoplanetary disk too cold for the main gas constituent, H₂, to emit, alternative tracers such as dust, CO, or the H₂ isotopolog HD are used. However, relying on disk mass measurements from any single tracer requires assumptions about the tracer's abundance relative to H₂ and the disk temperature structure. Using new Atacama Large Millimeter/submillimeter Array (ALMA) observations from the Molecules with ALMA at Planet-forming Scales (MAPS) ALMA Large Program as well as archival
ALMA observations, we construct a disk physical/chemical model of the protoplanetary disk GM Aur. Our model is in good agreement with the spatially resolved CO isotopolog emission from eleven rotational transitions with spatial resolution ranging from 0.15" to 0.46"(24-73 au at 159 pc) and the spatially unresolved HD J = 1 − 0 detection from Herschel. Our best-fit model favors a cold protoplanetary disk with a total gas mass of approximately 0.2 M$_{\odot}$, a factor of 10 reduction in CO gas inside roughly 100 au and a factor of 100 reduction outside of 100 au. Despite its large mass, the disk appears to be on the whole gravitationally stable based on the derived Toomre Q parameter. However, the region between 70 and 100 au, corresponding to one of the millimeter dust rings, is close to being unstable based on the calculated Toomre Q of < 1.7. This paper is part of the MAPS special issue of the Astrophysical Journal Supplement.

Keywords: Astrochemistry

# 1. INTRODUCTION

Protoplanetary disk mass is a fundamental property influencing virtually all aspects of the disk's evolution and the resulting planetary system. It sets a limit on the mass available to forming planets and determines the mechanisms that shape the final system architecture. Gravitational instability (GI) in a protoplanetary disk can result in the formation of massive companions at separations of hundreds of astronomical units from the central star (Boss 1997). Gravitational collapse at early stages results in the formation of close multiple-star systems (Tobin et al. 2018). In other cases, GI can result in the formation of massive planets at wide separation, such as the giant planets proposed to exist in the HD 163296 disk (Isella et al. 2016; Liu et al. 2018; Pinte et al. 2018; Teague et al. 2018).

In a recent study of the HD 163296 system, Booth
et al. (2019) found a total disk mass of 0.31 M<sub>☉</sub> based on
observations of the optically thin CO isotopolog <sup>13</sup>C<sup>17</sup>O.
They further found that the disk currently is stable
against gravitational collapse, though they noted the
disk may have been more massive, and thus unstable,
in the past, having had ample time to accrete mass onto
the central star. Analysis of the HL Tau disk using the
same <sup>13</sup>C<sup>17</sup>O transition indicates this much younger sys-
tem surrounded by a residual envelope has a lower total
disk mass of 0.2 M<sub>☉</sub> but is likely unstable at radii from
50-110 au (Booth & Ilee 2020), spanning several of the
rings and dark bands observed in the millimeter contin-
uum (ALMA Partnership et al. 2015).

Disk gas mass remains one of the most difficult pa-
rameters to constrain. This is because the dominant
gas component, H₂, does not readily emit throughout
most of the disk due to its lack of a permanent dipole
moment. Instead, trace species such as CO and dust are

used to extrapolate to a total gas mass. However, each
tracer relies on assumptions that may not be applicable
to protoplanetary disks (see Bergin & Williams (2017)
for a review).

Converting dust continuum emission to a total dust mass requires assumptions about the grain size distribution, dust composition, scattering, and dust temperature. The dust mass is often then converted to a gas mass assuming a gas-to-dust mass ratio of 100, as measured for the interstellar medium (ISM). However, several processes can change the derived ratio, including differential radial drift for grains of different sizes, dust growth beyond the observable range, accretion onto the central star, and photoevaporative winds. Additionally, observations show that the gas disk often extends far beyond the millimeter dust grains (Guilloteau & Dutrey 1998; Najita & Bergin 2018; Facchini et al. 2019; Trapman et al. 2019).

CO column densities derived from optically thin emis-
sion can be converted into a total gas mass assuming a
CO/H₂ ratio of ~ 10⁻⁴ based on ISM values and cor-
recting for the effects of CO freeze-out in the cold disk
midplane and CO photodissociation in the surface layers
(Miotello et al. 2014; Williams & Best 2014). However,
surveys of protoplanetary disks consistently find a dis-
crepancy between dust-derived and CO-derived disk gas
masses, with the CO-based measurements generally in-
dicating a lower mass (Ansdell et al. 2016; Long et al.
2017). One potential explanation is that gas-phase CO
abundance has been further reduced by processes be-
yond freeze-out and photodissociation, resulting in an
underestimation of the total gas. These processes could
be physical, such as vertical mixing, which preferentially
traps CO in the cold disk midplane (Xu et al. 2017; Krijt
et al. 2018), or chemical, with CO processed into other,
less emissive species in either the gas or ice (Yu et al.
2016; Bosman et al. 2018; Eistrup et al. 2018; Schwarz
et al. 2018).

<footnote>* NASA Hubble Fellowship Program Sagan Fellow</footnote>

<footnote>† NASA Hubble Fellow</footnote>
Alternative tracers are needed to determine the true gas mass in protoplanetary disks. One approach is to use the outer radius of dust emission at different wavelengths to constrain the rate of radial drift (Powell et al. 2017). Masses derived using this technique are significantly larger than those derived from dust or CO line emission (Powell et al. 2019). However, as this analysis does not yet consider disk substructure, which would impact drift timescales, these values are likely upper limits to the true gas mass.

Another approach is to derive the H₂ mass from observations of the H₂ isotopolog HD. The HD/H₂ ratio of 3 × 10⁻⁵ is not subject to the same processes that can change the CO/H₂ and gas-to-dust mass ratios (Linsky 1998). The HD *J* = 1 – 0 transition was detected in three protoplanetary disks using the *Herschel Space Observatory*, including MAPS target GM Aur (Bergin et al. 2013; McClure et al. 2016). Upper limits exist for an additional nineteen systems (McClure et al. 2016; Kama et al. 2020).

Initial analysis of the HD detections has yielded large disk gas masses, more in line with those derived from dust than with those derived from CO (McClure et al. 2016; Trapman et al. 2017). However, the range of disk masses consistent with the observed HD emission strength can be quite large, in some cases spanning more than an order of magnitude. This is due to the strong degeneracy between HD abundance and gas temperature in contributing to the observed HD emission strength. Further, due to the high *J* = 1 upper state energy, the ground state transition of HD does not emit at temperatures lower than roughly 20 K (Bergin et al. 2013). Knowledge of the gas temperature structure in the disk from, e.g., spatially resolved CO observations, can reduce the uncertainty on the HD derived gas mass from over an order of magnitude to approximately a factor of two (Trapman et al. 2017).

Subsequent analysis of the HD toward one source, TW Hya, has used observations of CO isotopologs to constrain the gas temperature and, when combined with HD, the gas mass (Favre et al. 2013; Schwarz et al. 2016; Zhang et al. 2017). Key results include the fact that HD emits primarily within the inner 20 au of the TW Hya disk, with a gas-to-dust mass ratio of 140 in this region, and a CO/H₂ abundance ranging from < 10⁻⁶ in the outer disk to greater than 10⁻⁵ inside the CO snowline, indicating an overall depletion of gas-phase CO in TW Hya as compared to ISM values. Additionally, Calahan et al. (2021) demonstrated a wide range of CO abundance and total gas mass are able to reproduce the observed CO emission profiles, while the additional

constraint of the HD line provides a way to break this degeneracy.

In this paper we focus on the disk around GM Aur. GM Aur is a 1.1 M☉ star hosting a well-known transition disk at a distance of 159 pc (Calvet et al. 2005; Hughes et al. 2009; Gaia Collaboration et al. 2018; Macías et al. 2018). We use the CO isotopologue observations of GM Aur from the Atacama Large Millimeter/submillimeter Array (ALMA) Large Program MAPS, along with the HD detection from Herschel, archival ALMA CO observations, and data on the spectral energy distribution (SED) to construct a 2D thermochemical model of the gas density, temperature, and CO abundance in the GM Aur disk. The observations and data reduction process are summarized in Section 2. Section 3 describes the modeling framework used to fit the data. In Section 4 we present the results of our modeling study. In Section 5 we discuss what our results reveal about the gas mass, temperature, and CO abundance in the GM Aur disk. Finally, our conclusions are summarized in Section 6.

## 2. OBSERVATIONS AND DATA REDUCTION

This study uses the CO isotopolog emission toward GM Aur as part of the MAPS Large Program, covering the ¹³CO, C¹⁸O, and C¹⁷O *J* = 1-0 transitions in Band 3 and the ¹²CO, ¹³CO, and C¹⁸O *J* = 2-1 transitions in Band 6. The full details of the calibration and imaging processes are described in Öberg & MAPS Team (2021) and Czekala & MAPS Team (2021) respectively. Additionally, we augment these data with CO isotopologs in Band 7 and Band 9 from the ALMA Cycle 4 program 2016.1.00565.S (PI K. Schwarz), targeting the ¹³CO and C¹⁸O *J* = 3-2 and ¹³CO, C¹⁸O and C¹⁷O *J* = 6-5 transitions. Observations were obtained in Band 7 on 11 November 2016 with 42 antennas. Observations were obtained in Band 9 on 18 August 2018 with 48 antennas. The continuum data associated with the Band 7 observations was previously analyzed by Macías et al. (2018).

Initial calibration of the archival Band 7 and Band 9 data was carried out by ALMA/NAASC staff using standard procedures. Additionally, phase and amplitude self-calibration were performed in each band using the continuum visibilities. The fixvis task was used to correct the phase center of each data set. Imaging was performed in CASA 5.4 using the imaging scripts developed by Czekala & MAPS Team (2021). The channel widths used in the imaging are 0.2 km s⁻¹ in Band 7 and 0.3 km s⁻¹ in Band 9. The properties of the final CLEAN images, with a robust weighting of 0.5, are given in Table 1, and the moment 0 maps are shown in Figure 1.
All lines are detected, including the C¹⁷O J = 6-5 transition, the first time this transition has been detected toward a protoplanetary disk. To help in constraining the disk gas mass we also consider the spatially unresolved HD J = 1-0 detection from Herschel. McClure et al. (2016) reported a 5σ detection toward GM Aur with a total integrated flux of 2.5 ± 0.5 × 10⁻¹⁸ W m⁻².

### 3. METHODS

#### 3.1. Physical Model

Our physical disk model is based on an axisymmetric viscously evolving disk (Lynden-Bell & Pringle 1974; Andrews et al. 2011). The surface densities of both gas and dust are described by

$$ \Sigma(R) = \Sigma_c \left(\frac{R}{R_c}\right)^{-\gamma} \exp\left[-\left(\frac{R}{R_c}\right)^{2-\gamma}\right], \quad (1) $$

where $\Sigma_c/e$ is the surface density at a characteristic radius $R_c$. The 2D density distribution is then

$$ \rho(R, Z) = \frac{\Sigma}{\sqrt{2\pi R h}} \exp\left[-\frac{1}{2}\left(\frac{Z}{h}\right)^2\right]. \quad (2) $$

The scale height $h$ varies as a function of radius such that

$$ h(R) = h_{\text{ref}} \left(\frac{R}{R_{\text{ref}}}\right)^{\psi} \quad (3) $$

where $h_{\text{ref}}$ is the characteristic scale height at a radius $R_{\text{ref}}$ and $\psi$ is the power-law index characterizing disk flaring. Our disk includes three populations of matter: gas, small grains, and large grains. The small grains are well mixed with the gas, i.e, have the same scale height, and follow the same surface density profile. The scale height for the large grains is smaller than that of the gas and small grains to mimic vertical settling. Further, while the gas and small grains vary smoothly, the large grain surface density includes several depleted regions corresponding to the gaps seen in continuum emission (Huang et al. 2020). For our large grain model we start with the surface density profile from the model of Macías et al. (2018). This model includes rings and gaps and was able to reproduce both the millimeter and centimeter continuum emission. The initial large grain distribution is then adjusted as described by Zhang & MAPS Team (2021).

#### 3.2. SED Fitting

Our initial values for the disk model were determined by fitting the SED. The general procedure for modeling the SEDs of all MAPS sources is described in detail by Zhang & MAPS Team (2021) and we use the same dust

surface density models for GM Aur as in that work. Briefly, data is fit using RADMC-3D (Dullemond et al. 2012). The large grains range in size from 0.005 µm to 1 mm with an MRN distribution $n(a) \propto a^{-3.5}$ (Mathis et al. 1977) and use the standard dust opacities from Birnstiel et al. (2018), who assumed a dust composed of 20% H₂O ice, 32.9% astronomical silicates, 7.4% troilite, and 39.7% refractory organics by mass. This differs from 38% graphite and 62% silicate composition used in previous SED modeling of GM Aur (Espaillat et al. 2011), necessitating some modification of the dust surface density profile from Macías et al. (2018). The small grains have a size range of 0.005 µm to 1 µm with an MRN distribution and are assumed to be composed of equal parts silicates and refractory organics by mass. At large scale heights, where most of the dust mass is in small grains, models suggest water is removed from grains via photodesorption (Hogerheijde et al. 2011). For a given large and small grain distribution, we used RADMC-3D to compare our model to the observed SED and ALMA continuum image.

#### 3.3. Thermochemical models

After using the SED to constrain the dust distribution we pass our disk density model to the 2D thermochemical modeling code RAC2D to model the molecular line emission. RAC2D self-consistently computes the chemistry as well as the evolving balance between heating and cooling in the disk. The modeling framework is described in detail by Du & Bergin (2014). For each model run RAC2D first calculates the dust thermal structure, cosmic-ray attenuation, and the radiation field, taking into account photon scatter and absorption. We assume a cosmic-ray ionization rate at the disk surface of 1.38 × 10⁻¹⁸ s⁻¹ per H₂, consistent with cosmic-ray modulation by stellar winds (Cleeves et al. 2013). Chemical evolution and gas temperature structure are then solved simultaneously.

The chemical network is based on the gas-phase network from the UMIST 2006 database (Woodall et al. 2007) and the grain surface network of Hasegawa et al. (1992). We consider a total of 5830 reactions among 484 species. The chemical network includes two body gas-phase reactions, photodissociation, including Ly-α dissociation of H₂O and OH, adsorption of species onto grain surfaces, thermal desorption, UV photo-desorption, and cosmic-ray induced desorption, as well as a limited network of two body grain surface reactions. The default initial chemical composition is given in Table 2.

The chemistry is run for 1 Myr. The main consequence of changing the run time is the amount of chem-
<table><caption>**Table 1.** GM Aur image parameters</caption><thead><tr><th rowspan="2">Molecular Transition</th><th>Beam</th><th>rms</th><th>Channel spacing</th><th rowspan="2">Program ID</th></tr><tr><th>(" × ", °)</th><th>(mJy beam<sup>-1</sup>)</th><th>(km s<sup>-1</sup>)</th></tr></thead><tbody><tr><td>C<sup>18</sup>O J=1-0</td><td>(0.30 × 0.30, 75.0)</td><td>0.514</td><td>0.5</td><td>2018.1.01055.L</td></tr><tr><td><sup>13</sup>CO J=1-0</td><td>(0.30 × 0.30, 81.1)</td><td>0.498</td><td>0.5</td><td>2018.1.01055.L</td></tr><tr><td>C<sup>17</sup>O J=1-0 F=<sup>3</sup>&frasl;<sub>2</sub> - <sup>5</sup>&frasl;<sub>2</sub>,</td><td>(0.30 × 0.30, -29.8)</td><td>0.629</td><td>0.5</td><td>2018.1.01055.L</td></tr><tr><td>C<sup>17</sup>O J=1-0 F=<sup>7</sup>&frasl;<sub>2</sub> - <sup>5</sup>&frasl;<sub>2</sub>,</td><td>(0.30 × 0.30, -29.8)</td><td>0.629</td><td>0.5</td><td>2018.1.01055.L</td></tr><tr><td>C<sup>17</sup>O J=1-0 F=<sup>5</sup>&frasl;<sub>2</sub> - <sup>5</sup>&frasl;<sub>2</sub>,</td><td>(0.30 × 0.30, -29.8)</td><td>0.629</td><td>0.5</td><td>2018.1.01055.L</td></tr><tr><td>C<sup>18</sup>O J=2-1</td><td>(0.15 × 0.15, 54.8)</td><td>0.484</td><td>0.2</td><td>2018.1.01055.L</td></tr><tr><td><sup>13</sup>CO J=2-1</td><td>(0.15 × 0.15, 72.4)</td><td>0.660</td><td>0.2</td><td>2018.1.01055.L</td></tr><tr><td>CO J=2-1</td><td>(0.15 × 0.15, 66.3)</td><td>0.730</td><td>0.2</td><td>2018.1.01055.L</td></tr><tr><td>C<sup>18</sup>O J=3-2</td><td>(0.37 × 0.27, 6.7)</td><td>9.2</td><td>0.2</td><td>2016.1.00565.S</td></tr><tr><td><sup>13</sup>CO J=3-2</td><td>(0.38 × 0.27, 7.5)</td><td>8.3</td><td>0.2</td><td>2016.1.00565.S</td></tr><tr><td>C<sup>18</sup>O J=6-5</td><td>(0.46 × 0.26, 3.8)</td><td>48.7</td><td>0.3</td><td>2016.1.00565.S</td></tr><tr><td><sup>13</sup>CO J=6-5</td><td>(0.46 × 0.26, 3.7)</td><td>56.1</td><td>0.3</td><td>2016.1.00565.S</td></tr><tr><td>C<sup>17</sup>O J=6-5</td><td>(0.45 × 0.26, 4.7)</td><td>57.1</td><td>0.3</td><td>2016.1.00565.S</td></tr></tbody></table>

<table><caption>**Table 2.** Standard Initial Chemical abundances</caption><thead><tr><th></th><th>Abundance Relative<br>to Total H</th></tr></thead><tbody><tr><td>H<sub>2</sub></td><td>5 × 10<sup>-1</sup></td></tr><tr><td>He</td><td>0.09</td></tr><tr><td>CO</td><td>1.4 × 10<sup>-4</sup></td></tr><tr><td>N</td><td>7.5 × 10<sup>-5</sup></td></tr><tr><td>H<sub>2</sub>O ice</td><td>1.8 × 10<sup>-4</sup></td></tr><tr><td>S</td><td>8 × 10<sup>-8</sup></td></tr><tr><td>Si<sup>+</sup></td><td>8 × 10<sup>-9</sup></td></tr><tr><td>Na<sup>+</sup></td><td>2 × 10<sup>-8</sup></td></tr><tr><td>Mg<sup>+</sup></td><td>7 × 10<sup>-9</sup></td></tr><tr><td>Fe<sup>+</sup></td><td>3 × 10<sup>-9</sup></td></tr><tr><td>P</td><td>3 × 10<sup>-9</sup></td></tr><tr><td>F</td><td>2 × 10<sup>-8</sup></td></tr><tr><td>Cl</td><td>4 × 10<sup>-9</sup></td></tr></tbody></table>

ical processing of CO that takes place, as the gas and dust temperatures converge on much shorter timescales. However, we adjust our initial CO abundance in order to match the observed CO line emission profiles, as described in Section 3.4. Thus, the decision on how long to run the chemistry does not impact our final results.

Line radiative transfer calculations assuming local thermal equilibrium are also carried out using RAC2D. Line and continuum emission are modeled together using

<table><caption>**Table 3.** Source Properties</caption><thead><tr><th></th><th>Value</th><th>Reference</th></tr></thead><tbody><tr><td>distance (pc)</td><td>159</td><td>Gaia Collaboration et al. (2018)</td></tr><tr><td><i>i</i> (°)</td><td>53.2</td><td>Huang et al. (2020)</td></tr><tr><td>PA (°)</td><td>57.2</td><td>Huang et al. (2020)</td></tr><tr><td>T<sub>eff</sub> (K)</td><td>4350</td><td>Espaillat et al. (2011)</td></tr><tr><td>M<sub>*</sub> (M<sub>☉</sub>)</td><td>1.1</td><td>Macías et al. (2018)</td></tr><tr><td>R<sub>*</sub> (R<sub>☉</sub>)</td><td>1.9</td><td>Macías et al. (2018)</td></tr></tbody></table>

the source properties in Table 3. The resulting image cubes are then continuum subtracted and convolved with an elliptical Gaussian beam with the same size and orientation as the corresponding observation before comparison to the observations. By treating the dust continuum and line emission simultaneously, we account for any extinction of the line emission by the dust. Our chemical network is not fractionated to include species such as <sup>13</sup>C, <sup>18</sup>O, and deuterium. Instead, isotopolog emission profiles are generated assuming <sup>13</sup>C/<sup>12</sup>C = 69, <sup>16</sup>O/<sup>18</sup>O = 557, <sup>18</sup>O/<sup>17</sup>O = 3.6, and D/H = 1.5 × 10<sup>-5</sup> as measured in the local ISM (Linsky 1998; Wilson 1999).

### 3.4. *Parameter study*

The mass in dust is constrained by the continuum imaging and the SED. In attempting to fit the line emission we limit our parameter study to the total gas mass and the variables that determine the gas and dust den-
<figure><img src="image_1.png" /><figcaption>**Figure 1.** Integrated intensity maps with a logarithmic color stretch for the 11 CO isotopolog rotational transitions considered in this work. The C$^{17}$O *J* = 1 − 0 map is a sum of the three hyperfine components. The *J* = 1 − 0 and *J* = 2 − 1 transitions were observed as part of MAPS. The *J* = 3 − 2 and 6-5 transitions were observed as part of ALMA program 2016.1.00565.S.</figcaption></figure>

<table><caption>**Table 4.** Range of parameter values considered</caption><thead><tr><th></th><th>Gas</th><th>Small Dust</th><th>Large Dust</th></tr></thead><tbody><tr><td>Mass (M<sub>☉</sub>)</td><td>0.02-0.41</td><td>1.03 &times; 10<sup>-4</sup></td><td>5.94 &times; 10<sup>-4</sup></td></tr><tr><td>&psi;<sup>a</sup></td><td>1-2</td><td>1-2</td><td>1-2</td></tr><tr><td>&gamma;</td><td>0.3-1.5</td><td>0.3-1.5</td><td>&hellip;</td></tr><tr><td>R<sub>c</sub> (au)</td><td>100-176</td><td>100-176</td><td>&hellip;</td></tr><tr><td>h<sub>ref</sub> (au)</td><td>5-12</td><td>5-12</td><td>0.75-12</td></tr><tr><td>R<sub>in</sub> (au)</td><td>0.5-27</td><td>&hellip;</td><td>&hellip;</td></tr><tr><td>CO/H</td><td>7 &times; 10<sup>-7</sup> &ndash; 1.4 &times; 10<sup>-3</sup></td><td>&hellip;</td><td>&hellip;</td></tr></tbody></table>

$^{a}$ψ for the gas and small grains are varied together, while ψ for the large grains is changed independently.

sity distribution: γ, R$_{c}$, h$_{ref}$, and ψ, as well as the initial CO abundance and the gas temperature in the inner disk, as discussed in Section 4. In total we generate 145

unique models. The range of parameters considered are given in Table 4. Due to the long run times required for each model, it is unrealistic to use a systematic parameter study using, e.g., a Markov Chain Monte Carlo method to find the best-fit model. Instead, parameters are changed one at a time in order to achieve a reasonable fit.

We modify the CO abundance before modeling the chemistry, in contrast with the CO abundance study of Zhang & MAPS Team (2021), who modify the CO abundance after modeling the chemistry. Because RAC2D includes line processes when calculating the gas temperature, our choice to remove CO initially impacts the gas temperature structure and, by extension, the strength of the HD *J* = 1 − 0 emission. In constructing our initial model we use the same disk parameters as the best-fit model from Zhang & MAPS Team (2021). Initially, we use the CO depletion profile of Zhang & MAPS Team (2021). We then run additional thermochemical models
<table><caption>**Table 5.** Gas and Dust Population Parameters: Initial Model Values</caption><thead><tr><th></th><th>Gas</th><th>Small Dust</th><th>Large Dust</th></tr></thead><tbody><tr><td>Mass (M<sub>☉</sub>)</td><td>0.2</td><td>1.03 × 10<sup>-4</sup></td><td>5.94 × 10<sup>-4</sup></td></tr><tr><td>ψ</td><td>1.35</td><td>1.35</td><td>1.35</td></tr><tr><td>γ</td><td>1.0</td><td>1.0</td><td>1.0</td></tr><tr><td>R<sub>c</sub> (au)</td><td>176</td><td>111</td><td>...</td></tr><tr><td>R<sub>ref</sub> (au)</td><td>100</td><td>100</td><td>100</td></tr><tr><td>h<sub>ref</sub> (au)</td><td>7.5</td><td>7.5</td><td>3.75</td></tr><tr><td>R<sub>in</sub> (au)</td><td>1.0</td><td>1.0</td><td>34</td></tr><tr><td>R<sub>out</sub> (au)</td><td>650</td><td>650</td><td>310</td></tr></tbody></table>

with the depletion profile multiplied by a constant factor ranging from 0.1 to 2. To quantify how well a given model fits the data we calculate the reduced $\chi^2$ for the CO isotopolog emission radial profiles, comparing the model and observed emission at half beam spacing. We construct a new initial CO abundance profile taking the best fit based on the reduced $\chi^2$ at each radius. This model, using an updated CO abundance profile but otherwise using the same parameters as the best-fit model from Zhang & MAPS Team (2021), serves as our initial model.

To constrain the model in the vertical direction we compare the extracted emission surfaces from the observations and models for the C$^{18}$O $J = 2-1$, $^{13}$CO $J = 2-1$, $^{12}$CO $J = 2-1$, and $^{13}$CO $J = 3-2$ lines. The signal-to-noise ratios of the other transitions are too low to meaningfully constrain the emission height. Emission surfaces for both the observations and models are extracted with the Python package `disksurf`<fnref n="1" />, using the method presented by Pinte et al. (2018). In regions where the line flux is weak or originating from a large vertical range, i.e., is optically thin, there is greater uncertainty in the derived emission surface for both the models and the observations. For a detailed discussion of this technique as it applies to the MAPS data see Law & MAPS Team (2021). Finally, we compare the results of each model to the total observed HD flux.

## 4. RESULTS

The parameters for our initial model, based on the best-fit model of Zhang & MAPS Team (2021), are given in Table 5. This model fits the radial intensity profiles for the majority of the observed lines within $1\sigma$ outside of 160 au ($1''$) (Figure 2). Inside of 160 au, the

model under-predicts the line flux for nearly all transitions. The integrated HD $J = 1-0$ flux in our initial model is $1.9 \times 10^{-18} \text{ W m}^{-2}$ compared to the observed $2.5 \pm 0.5 \times 10^{-18} \text{ W m}^{-2}$, just below the $1\sigma$ uncertainty. In the outer disk the model emission surfaces for the $^{13}$CO lines are below the $1\sigma$ uncertainty of the surfaces derived from observations, while for C$^{18}$O the model over-predicts the emission surface (Figure 3).

To raise outer disk emission surfaces in our model, we modify the surface density profile of the gas and small grains by changing $R_c$ and $\gamma$, thus shifting more mass to larger radii, and the disk flaring by adjusting $\psi$. We also consider models with varying total gas mass (see Table 4) but find that holding the gas mass at 0.2 M$_{\odot}$ provides the best fit to the data. After adjusting the gas surface density, we modify the initial CO depletion profile as described in Section 3.4 in order to match the observed radial emission profiles. The resulting model brings the C$^{18}$O 2-1, $^{13}$CO 2-1, and $^{13}$CO 3-2 model emission surfaces within the $1\sigma$ uncertainty of the surfaces derived from observations throughout much of the disk (Figure 3). Further, the model HD flux increases from $1.9 \times 10^{-18} \text{ W m}^{-2}$ in our initial model to $2.5 \times 10^{-18} \text{ W m}^{-2}$, in excellent agreement with the observed flux of $2.5 \pm 0.5 \times 10^{-18} \text{ W m}^{-2}$.

In the inner disk the model under-predicts the strength of the $^{12}$CO 2-1 emission by a factor of 2.6, while providing a reasonable fit to the emission from the less abundant CO isotopologs (Figure 4). The emission surface profiles derived from observations show that $^{12}$CO emission originates from higher in the disk than the other lines. To fit the $^{12}$CO $J = 2-1$ emission in the inner disk we increase the temperature inside of 32 au and for $Z/R > 0.1$. Increasing the gas temperature by a factor of 10 in this region, from several tens of kelvin to several hundred kelvin, greatly improves the agreement between the model $^{12}$CO $J = 2-1$ emission and the observations without significantly increasing emission from the other lines (Figure 5). Possible sources of this extra heating are discussed below. Increasing the temperature in the inner disk has a negligible effect on the model HD flux. This model, with an increased temperature in the inner disk, provides the best fit to the data. The input values for our best-fit model are shown in Table 6 and the 2D hydrogen gas distribution, dust temperature, and gas temperature are shown in Figure 6. However, because of our sparsely sampled parameter space we cannot rule out the possibility of other model solutions fitting the data equally well.

In order to match the $^{12}$CO 2-1 flux in the inner disk, we need to increase the gas temperature in the surface layers of the inner disk. The primary source of dust

<footnote n="1">¹ https://github.com/richteague/disksurf</footnote>
<figure><img src="image_2.png" /><figcaption>**Figure 2.** Red lines show the deprojected, azimuthally averaged radial emission profiles for our model using the same disk density parameters as Zhang & MAPS Team (2021). Blue lines show the observations. Blue shading indicates the 1σ uncertainty. Light gray lines show the profiles for all models.</figcaption></figure>

<table><caption>**Table 6.** Gas and Dust Population Parameters: best-fit model Values</caption><thead><tr><th></th><th>Gas</th><th>Small Dust</th><th>Large Dust</th></tr></thead><tbody><tr><td>Mass (M<sub>☉</sub>)</td><td>0.2</td><td>1.03 × 10<sup>-4</sup></td><td>5.94 × 10<sup>-4</sup></td></tr><tr><td>ψ</td><td>1.5</td><td>1.5</td><td>2.0</td></tr><tr><td>γ</td><td>0.59</td><td>0.59</td><td>1.0</td></tr><tr><td>R<sub>c</sub> (au)</td><td>111</td><td>111</td><td>...</td></tr><tr><td>R<sub>ref</sub> (au)</td><td>100</td><td>100</td><td>100</td></tr><tr><td>h<sub>ref</sub> (au)</td><td>7.5</td><td>7.5</td><td>1.0</td></tr><tr><td>R<sub>in</sub> (au)</td><td>15</td><td>1.0</td><td>0.45</td></tr><tr><td>R<sub>out</sub> (au)</td><td>650</td><td>650</td><td>600</td></tr></tbody></table>

heating in RAC2D is radiation from the central star. The gas temperature is initially assumed to be the same as the dust temperature and allowed to evolve due to a number of heating and cooling processes, including photoelectric heating, endothermic and exothermic chemical reactions, and viscous dissipation (Du & Bergin 2014).

However, in the surface layers of the inner disk photoelectric heating of polycyclic aromatic hydrocarbons (PAHs) can be an important contributor to the local gas temperature (Kamp & Dullemond 2004; Woitke et al. 2016). While heating by PAHs is included in RAC2D, the PAH abundance in the disk surface layers is uncertain. We assume a PAH abundance relative to H of $1.6 \times 10^{-7}$ (Du & Bergin 2014). A higher PAH abundance could result in increased photoelectric heating and thus in a warmer disk surface.

Alternatively, mechanical heating mechanisms such as stellar winds and accretion onto the central star can raise the temperature of the inner disk. When including mechanical heating, both Glassgold et al. (2004) and Najita & Ádámkovics (2017) found temperatures of several hundred kelvin for vertical column densities equivalent to those in the region of our model where we artificially increase the temperature, though these models looked at full disks without a large central dust cavity as seen in the GM Aur disk. Calahan & MAPS Team (2021) also found that their RAC2D model requires additional
<figure><img src="image_3.png" /><figcaption>**Figure 3.** Extracted emission surface profiles for the ¹²CO 2-1, ¹³CO 3-2, ¹³CO 2-1, and C¹⁸O 2-1 lines in our observations (blue) and three models (red). Shading indicates the uncertainty based on the scatter of points in the extraction before averaging. From left to right the columns show our results for the model using the disk parameters from Zhang & MAPS Team (2021), our best-fit model, and our best-fit model with small grain settling.</figcaption></figure>

heating to match the observed ¹²CO J = 2 – 1 flux in the HD 163296 disk inside of 32 au, further suggesting that under-predicting emission in the surface layers of the inner disk is due to limitations in the code.

## 4.1. CO Abundance

Figure 7 shows the CO abundance relative to H₂ as a function of height and radius in our best-fit model. CO is effectively frozen out from the gas near the mid-plane from 75 to 250 au. Outside the millimeter dust disk, the gas temperature remains below the CO freeze-out temperature. However, nonthermal desorption by UV photons allows some CO to remain in the gas (e.g., Öberg et al. 2015). In the region directly above the CO freeze-out layer, gas-phase CO has been converted into CO₂ ice (e.g., Reboussin et al. 2015; Bosman et al. 2018). Closer to the disk surface, where the temperature exceeds the CO₂ freeze-out temperature, CO remains in the gas.

Figure 8 compares the CO depletion profile after evolving the chemistry in our best-fit model for 1 Myr to that found by Zhang & MAPS Team (2021). Both results follow the general trend of a roughly constant, high level of CO depletion outside of roughly 100 au, with the inner disk less depleted in CO. The location of the mid-plane CO snowline in our model, here defined as where the CO gas and ice abundances are equal, is 31 au, consistent with the derived CO snowline of 30 ± 5 au from Zhang & MAPS Team (2021). Our model has greater CO depletion inside 200 au, as well as a more abrupt return of CO in the inner disk compared to the Zhang & MAPS Team (2021) results. The abrupt change in the CO column is due to the conversion of CO gas into CO₂ ice and CH₄ ice near the midplane from roughly 90-150 au, which is not seen to the same extent in the Zhang & MAPS Team (2021) model. Given that the two approaches remove CO from the disk at different points in the modeling process, and that this work attempts to
<figure><img src="image_4.png" /></figure>

**Figure 4.** Red lines show the deprojected, azimuthally averaged radial emission profiles for our best-fit model prior to adjusting the inner disk temperature are. Blue lines show the observations. Blue shading indicates the $1\sigma$ uncertainty. Light gray lines show the profiles for all models.

match a greater number of observations, some variation is to be expected.

## 5. DISCUSSION

### 5.1. 2D Temperature & CO Distribution

Figure 9 shows the 2D gas temperature distribution in a subset of our models with the derived temperatures for the $^{12}$CO $J = 2-1$, $^{13}$CO $J = 2-1$, and $^{13}$CO $J = 3-2$ surfaces over-plotted. The C$^{18}$O $J = 2-1$ is not included as this emission is optically thin and thus is not a good temperature tracer. The temperature extraction follows the same process described by Law & MAPS Team (2021). Briefly, the gas temperature is determined from the peak surface brightness at a given radius for the non-continuum subtracted line image cubes using the full Planck function and assuming the line emission is optically thick. Using non-continuum subtracted data to measure the temperature ensures the temperature is not underestimated in the case of optically thick dust emission (Weaver et al. 2018). In the layers traced by the $^{13}$CO $J = 2-1$ and $J = 3-2$ lines our best-fit

model temperature is in reasonable agreement with the data, with the model gas temperature varying by less than 10 K from the measured temperature at most radii (Figure 10).

However, at larger heights our model is over two times warmer than the temperature measured from the $^{12}$CO $J = 2-1$ emission surface. Past models of CO emission in protoplanetary disks have cooled the upper layers of the disk by increasing dust settling (McClure et al. 2016; Calahan et al. 2021). We test this solution by decreasing the scale height of the small dust grains in our model. Decreasing the small grain scale height from 7.5 to 5 au does in fact decrease the gas temperature in some regions of the disk, particularly in the upper layers beyond 400 au (Figure 9). However, the regions traced by CO emission are warmer in this model, increasing the discrepancy between the model and the observations (Figure 10).

An alternative explanation is that the upper layers of the disk are warmer and more CO rich than we assume in our models. If material is falling onto the disk from a
<figure><img src="image_5.png" /><figcaption>**Figure 5.** Red lines show the deprojected, azimuthally averaged radial emission profiles for our best-fit model after adjusting the inner disk temperature. Blue lines show the observations. Blue shading indicates the $1\sigma$ uncertainty. Light gray lines show the profiles for all models.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 6.** Maps of the H number density (left), dust temperature (center), and gas temperature (right) in our best-fit model.</figcaption></figure>

residual cloud or envelope, the favored explanation for the nonaxisymmetric features seen in $^{12}$CO $J = 2-1$ (Huang & MAPS Team 2021), the infalling material is expected to produce shock heating (Sakai et al. 2014). This heating will enhance the gas temperature in the surface layers from which $^{12}$CO emits. Additionally, the infalling material is unlikely to have undergone much chemical processing, and thus will have a CO/H$_2$ abun-

dance ratio similar to that of the dense ISM. While the amount of CO supplied to the disk by a residual envelope is likely to be small, the combination of increased temperature and elevated CO/H$_2$ could result in some $^{12}$CO emission originating closer to the disk surface.
<figure><img src="image_7.png" /><figcaption>**Figure 7.** Map of CO abundance relative to H₂ in our best-fit model.</figcaption></figure>

Figure 11 shows the distribution of the HD emission in our best-fit model. Seventy-five percent of the HD emission originates from the inner 100 au. In particular, the hot, low density gas inside the millimeter dust inner radius at 32 au contributes a non-negligible amount to the HD flux. In comparison, only 47% of the total disk gas mass is inside 100 au. Since HD does not readily emit at temperatures less than ~20 K, much of the disk beyond 100 au is not well traced by the HD *J* = 1–0 emission. There may be more mass in the outer disk than accounted for in our best-fit model. Additional analysis, e.g., using CS emission (Teague et al. 2018), is needed to better constrain the gas density in the outer disk.

Previous analysis of the HD detection by McClure et al. (2016) in GM Aur constrains the disk gas mass to 0.025-0.204 $M_{\odot}$. Based on analysis of the millimeter and centimeter continuum, Macías et al. (2018) found a total dust mass of 2 $M_{\odot}$. Assuming a gas-to-dust mass ratio of 100, this corresponds to a total gas mass of 0.19 $M_{\odot}$, at the high end of the range given by McClure et al. (2016). Trapman et al. (2017) reanalyzed the HD detection in GM Aur by comparing to the HD line flux in a grid of generic 2D thermochemical models. They constrained the mass of the GM Aur disk to be between 0.01 $M_{\odot}$ and a few tenths of a solar mass. Woitke et al. (2019) have also built a model of the GM Aur disk as part of the DIANA project. Their model based on the observed SED has a total disk mass of 0.11 $M_{\odot}$ and also reasonably reproduces the observed total flux of the 63 µm [OI] line as well as the $^{12}$CO 2-1 and HCO$^{+}$ 3-2 lines. Their final model, independent of the SED fit and based on Submillimeter Array (SMA) observations of the $^{12}$CO *J* = 2-1 line, has a disk mass of 3.3 × 10$^{-2}$ $M_{\odot}$.

Our best-fit model has a total gas mass of 0.2 $M_{\odot}$, consistent with the upper limits of previous works. This

<figure><img src="image_8.png" /><figcaption>**Figure 8.** Top: CO column density profile for our best-fit model and for the best-fit model of Zhang & MAPS Team (2021). Bottom: CO depletion factor as a function of radius for our best-fit model after evolving the chemistry for 1 Myr, as well as the best-fit model of Zhang & MAPS Team (2021) relative to the initial CO abundance. Vertical gray line in the midplane CO snowline at 31 au in our model. Horizontal dashed lines indicate depletion factors of 10 and 100.</figcaption></figure>

high value comes in part from our overall low disk temperature, needed to match the gas temperature from CO observations. Additionally, previous works have set the outer radius of the disk to 250-300 au. Here we set the disk outer radius to 650 au based on the observed extent of the $^{12}$CO 2-1 emission, though less than 3% of the total gas mass is outside of 300 au.

### 5.3. Disk Stability

The stability of a rotating disk against gravitational collapse is often characterized using the Toomre $Q$ pa-
<figure><img src="image_9.png" /><figcaption>Figure 9. Temperatures of the emission surfaces derived from ¹²CO, and ¹³CO 2-1 the ¹³CO 3-2 observations overlaid on the gas temperature map from our best-fit model.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Figure 10. Temperatures of the emission surfaces derived from ¹²CO 2-1, ¹³CO 2-1, and ¹³CO 3-2 observations (blue) compared to the temperature of the corresponding location in our models (red).</figcaption></figure>

parameter (Toomre 1964):

$$Q = \frac{c_s \Omega}{\pi G \Sigma} \qquad (4)$$

where $c_s$ is the gas sound speed assuming the midplane gas temperature, $\Omega$ is the Keplerian angular velocity of the disk, and $\Sigma$ is the total gas+dust surface density. For a geometrically thin disk $Q \sim 1$ is needed for density perturbations to develop. However, numerical simulations demonstrate that instabilities can develop for $Q < 1.7$ in systems not well approximated by a geometrically thin disk (Helled et al. 2014). Figure 12 shows the Toomre $Q$ radial profile using values from our best fit disk model.

Our calculated $Q$ value for GM Aur is greater than 1.7 throughout much of the disk. In the outer disk, where spiral-like features are seen in the ¹²CO 2-1 emission, our calculated $Q$ is extremely high, ranging from ~7 at 250 au to ~700 at 500 au. The two values in

our model that determine $Q$ are temperature and surface density. Without changing the surface density, the midplane temperature at 250 au would need to be unrealistically low, ~1 K, to reach a $Q$ of 1.7.

As discussed in the previous section, the HD detection does not constrain the gas surface density in the outer disk. The model surface density would need to be increased by a factor of four at 250 au to reach $Q$ of 1.7 while holding the temperature constant. Zhang & MAPS Team (2021) note that for four out of the five MAPS sources, including GM Aur, the CO column density profiles are very shallow, consistent with a viscously evolving disk. This supports the conclusion of Huang & MAPS Team (2021), who argue against the nonaxisymmetric features seen in ¹²CO $J = 2-1$ being driven exclusively by disk instability based on the kinematics.

However, $Q$ dips below 1.7 from 70 to 100 au, corresponding to the location of one of the bright rings seen in
<figure><img src="image_11.png" /><figcaption>**Figure 11.** Top: HD *J* = 1 − 0 emitting region in our best-fit model overlayed on a map of the gas temperature. The gray contour contains the middle 75% of the HD emission. Middle: Deprojected, azimuthally averaged radial profile for the HD *J* = 1 − 0 emission. Since the observed HD emission is spatially unresolved, the model is not convolved with a beam. Bottom: Plot showing the total HD flux (blue) and gas mass (red) interior to a given radius. The HD emission preferentially originates from the warm inner disk.</figcaption></figure>

the continuum. The concentration of large dust grains in
this region increases the disk opacity and thus decreases
the temperature of both the gas and the dust. This lower
temperature, in turn, results in a lower sound speed and
thus in a lower Toomre *Q* value. The presence of dust

rings and gaps can lower the midplane temperature in
a dust ring by several kelvin compared to a disk with a
smoothly varying surface density profile (Facchini et al.
2018; van der Marel et al. 2018; Alarcón et al. 2020;
Calahan & MAPS Team 2021).

The dip in the Toomre *Q* in our model is due entirely to an over-density of dust. Our model is able to fit the CO emission profiles in this region without a corresponding increase in the gas density. Decreasing the disk surface density to 88% of our assumed value between 70 and 100 au would bring the dust ring into a gravitationally stable regime. Alternatively, a warmer temperature than assumed would also lead to a higher *Q*. Between 70 and 100 au our model midplane temperature is ~ 12 K. The temperature derived from the observed ¹³CO 2-1 line is ~ 22 K and can be considered an upper limit on the midplane temperature. We use the temperature from ¹³CO because the C¹⁸O is optically thin at these radii and therefore not a good tracer of temperature. Taking the midplane temperature to be 22 K increases the *Q* value to 2.0. Interestingly, recent smoothed particle hydrodynamics modeling shows that a migrating planet can increase the local disk temperature, suppressing spiral structure and stabilizing the disk (Rowther et al. 2020).

GI is thought to primarily manifest as nonaxisymmet-
ric features. Previous analysis of the continuum disk
emission at FUV wavelengths does not indicate any
such features in the 70-120 au range (Hornbeck et al.
2016). No nonaxisymmetric substructure is seen in the
CO emission profiles in this region. However, the in-
tensity of the bright continuum ring at 40 au at 7 mm
shows a low signal-to-noise (~2σ) asymmetry (Macías
et al. 2018). While the ring at 84 au is not detected at
7 mm with high enough sensitivity to enable a similar
analysis, Huang et al. (2020) note that the 84 au ring
is nonaxisymmetric at 1.1 mm. The ring appears wider
along the major axis of the disk, which as Huang et al.
(2020) demonstrate is unlikely to be an imaging arti-
fact; instead the variation can be attributed to either
nonaxisymmetric or vertical structure within the ring.

An interesting point of comparison is the HL Tau disk,
which, though less evolved than GM Aur, has a similar
total gas mass and a region of instability centered on a
dust gap (Booth & Ilee 2020). Spiral structure is also
seen in the HCO⁺ *J* = 3-2 emission toward HL Tau (Yen
et al. 2019). While this spiral structure was originally
attributed to the infalling envelope, Booth & Ilee (2020)
note that the feature could also be associated with the
region of instability in the disk. Conversely, observa-
tions of HCO⁺ *J* = 3-2 toward GM Aur do not appear
to deviate substantially from Keplerian rotation (Huang
<figure><img src="image_12.png" /></figure>

**Figure 12.** Top: Calculated Toomre Q value. Dotted line indicates the gravitationally unstable threshold of 1.7 for a geometrically thick disk. Bottom: Midplane temperature as a function of radius in our best-fit model. Dashed line is the gas temperature derived from the ¹³CO 2-1 observations. Gray regions indicate the locations of the observed gaps in the millimeter continuum (Macías et al. 2018). The dips in temperature, and, by extension, Q, corresponds to bright rings in the continuum.

et al. 2020). It is possible that the GM Aur disk is in the process of stabilizing after a period of infall or planet formation. Given the limitations of using Toomre Q to determine the stability of the non-geometrically thin disk, a more detailed study of the kinematics is required to determine the stability of the GM Aur disk from 70 to 100 au.

## 6. SUMMARY AND CONCLUSIONS

In this work we use observations of CO isotopologs in the GM Aur disk taken as part of the MAPS ALMA

Large Program along with archival observations of CO from ALMA and HD from Herschel to build a model of the disk gas density and temperature structure. Based on our results we conclude the following:

* Much of the disk (32% by mass) is cooler than 20 K. As such the HD emission only traces the inner 200 au, while the gas disk extends to 650 au based on observations of ¹²CO.

* We constrain the gas mass of the GM Aur disk to be ~ 0.2 M☉. While the total mass in the outer disk remains somewhat uncertain, only 15% of the mass in our best-fit model is beyond 200 au. Any variation in mass in the outer disk will likely have only a small effect on the total disk mass.

* The CO gas abundance relative to H₂ is reduced by approximately one order of magnitude with respect to the ISM values inside 100 au and by two orders of magnitude outside 100 au. This is consistent with the analysis of Zhang & MAPS Team (2021). Our model also shows CO gas returning to the midplane outside of the millimeter dust disk due to nonthermal desorption.

* Based on the calculated Toomre Q parameter, the outer disk appears stable against gravitational collapse. However, Q dips into the unstable regime between 70 and 100 au, corresponding to the second bright ring seen in millimeter dust emission. While there is some evidence for nonaxisymmetric features in the dust continuum at these radii, a more detailed study is needed to determine if the GM Aur disk is gravitationally unstable.

## ACKNOWLEDGMENTS

This paper makes use of the following ALMA data: ADS/JAO.ALMA#2016.1.00565.S and ADS/JAO.ALMA#2018.1.01055.L. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), MOST and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. This work is based on observations made with Herschel, a European Space Agency Cornerstone Mission with significant participation by NASA.
K.R.S., K.Z., J.B., J.H., and I.C. acknowledge the support of NASA through Hubble Fellowship Program grants HST-HF2-51419.001, HST-HF2-51401.001, HST-HF2-51427.001-A, HST-HF2-51460.001-A, and 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.K.C. acknowledges support from the National Aeronautics and Space Administration FINESST grant, under Grant no. 80NSSC19K1534. C.J.L. and J.K.C. acknowledge funding from the National Science Foundation Graduate Research Fellowship under Grant No. DGE1745303 and DGE1256260. 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. Y.A. acknowledges support by NAOJ ALMA Scientific Research Grant Code 2019-13B, and Grant-in-Aid for Scientific Research Nos. 18H05222 and 20H05847. S.A. and J.H. 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., A.D.B., and F.A. acknowledge support from NSF AAG Grant #1907653. A.S.B. acknowledges the studentship funded by the Science and Technology Facilities Council of the United Kingdom (STFC). G.C. is supported by NAOJ ALMA Scientific Research grant Code 2019-13B. L.I.C. gratefully acknowledges support from the David and Lucille Packard Foundation and Johnson & Johnson's WiSTEM2D Program. J.D.I. acknowledges support from the Science and Technology Facilities Council of the United Kingdom (STFC) under ST/T000287/1. R.L.G. acknowledges support from a CNES research fellowship. F.L. and R.T. acknowledge support from the Smithsonian Institution as Submillimeter Array (SMA) Fellows. Y.L. acknowledges the financial support by the Natural Science Foundation of China (Grant No. 11973090). 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"). K.I.O. acknowledges support from the Simons Foundation (SCOL #321183) and an NSF AAG Grant (#1907653). C.W. acknowledges financial support from the University of Leeds, STFC and UKRI (grant numbers ST/R000549/1, ST/T000287/1, MR/T040726/1).

Facilities: ALMA, Herschel

Software: analysisUtils (https://casaguides.nrao.edu/index.php/Analysis_Utility), astropy (Astropy Collaboration et al. 2013), CASA (McMullin et al. 2007), `diskprojection` (https://github.com/richteague/disksurf), `matplotlib` (Hunter 2007), ` numpy` (van der Walt et al. 2011), `RAC2D` (Du & Bergin 2014), `RADMC-3D`, (Dullemond et al. 2012), `SciPy` (Jones et al. 2001-),
REFERENCES

Alarcón, F., Teague, R., Zhang, K., Bergin, E. A., & Barraza-Alfaro, M. 2020, ApJ, 905, 68, doi: 10.3847/1538-4357/abc1d6

ALMA Partnership, Brogan, C. L., Pérez, L. M., et al. 2015, ApJL, 808, L3, doi: 10.1088/2041-8205/808/1/L3

Andrews, S. M., Wilner, D. J., Espaillat, C., et al. 2011, ApJ, 732, 42, doi: 10.1088/0004-637X/732/1/42

Ansdell, M., Williams, J. P., van der Marel, N., et al. 2016, ApJ, 828, 46, doi: 10.3847/0004-637X/828/1/46

Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33, doi: 10.1051/0004-6361/201322068

Bergin, E. A., & Williams, J. P. 2017, The Determination of Protoplanetary Disk Masses, ed. M. Pessah & O. Gressel, Vol. 445, 1, doi: 10.1007/978-3-319-60609-5_1

Bergin, E. A., Cleeves, L. I., Gorti, U., et al. 2013, Nature, 493, 644, doi: 10.1038/nature11805

Birstiel, T., Dullemond, C. P., Zhu, Z., et al. 2018, ApJL, 869, L45, doi: 10.3847/2041-8213/aaf743

Booth, A. S., & Ilee, J. D. 2020, MNRAS, 493, L108, doi: 10.1093/mnrasl/slaa014

Booth, A. S., Walsh, C., Ilee, J. D., et al. 2019, ApJL, 882, L31, doi: 10.3847/2041-8213/ab3645

Bosman, A. D., Walsh, C., & van Dishoeck, E. F. 2018, A&A, 618, A182, doi: 10.1051/0004-6361/201833497

Boss, A. P. 1997, Science, 276, 1836, doi: 10.1126/science.276.5320.1836

Calahan, J., & MAPS Team. 2021, ApJ, 0, 0, doi: 0

Calahan, J. K., Bergin, E., Zhang, K., et al. 2021, ApJ, 908, 8, doi: 10.3847/1538-4357/abd255

Calvet, N., D'Alessio, P., Watson, D. M., et al. 2005, ApJL, 630, L185, doi: 10.1086/491652

Cleeves, L. I., Adams, F. C., & Bergin, E. A. 2013, ApJ, 772, 5, doi: 10.1088/0004-637X/772/1/5

Czekala, I., & MAPS Team, A. 2021, ApJ, 0, 0, doi: 0

Du, F., & Bergin, E. A. 2014, ApJ, 792, 2, doi: 10.1088/0004-637X/792/1/2

Dullemond, C. P., Juhasz, A., Pohl, A., et al. 2012, RADMC-3D: A multi-purpose radiative transfer tool. http://ascl.net/1202.015

Eistrup, C., Walsh, C., & van Dishoeck, E. F. 2018, A&A, 613, A14, doi: 10.1051/0004-6361/201731302

Espaillat, C., Furlan, E., D'Alessio, P., et al. 2011, ApJ, 728, 49, doi: 10.1088/0004-637X/728/1/49

Facchini, S., Pinilla, P., van Dishoeck, E. F., & de Juan Ovelar, M. 2018, A&A, 612, A104, doi: 10.1051/0004-6361/201731390

Facchini, S., van Dishoeck, E. F., Manara, C. F., et al. 2019, A&A, 626, L2, doi: 10.1051/0004-6361/201935496

Favre, C., Cleeves, L. I., Bergin, E. A., Qi, C., & Blake, G. A. 2013, ApJL, 776, L38, doi: 10.1088/2041-8205/776/2/L38

Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2018, A&A, 616, A1, doi: 10.1051/0004-6361/201833051

Glassgold, A. E., Najita, J., & Igea, J. 2004, ApJ, 615, 972, doi: 10.1086/424509

Guilloteau, S., & Dutrey, A. 1998, A&A, 339, 467

Hasegawa, T. I., Herbst, E., & Leung, C. M. 1992, ApJS, 82, 167, doi: 10.1086/191713

Helled, R., Bodenheimer, P., Podolak, M., et al. 2014, in Protostars and Planets VI, ed. H. Beuther, R. S. Klessen, C. P. Dullemond, & T. Henning, 643, doi: 10.2458/azu_uapress_9780816531240-ch028

Hogerheijde, M. R., Bergin, E. A., Brinch, C., et al. 2011, Science, 334, 338, doi: 10.1126/science.1208931

Hornbeck, J. B., Swearingen, J. R., Grady, C. A., et al. 2016, ApJ, 829, 65, doi: 10.3847/0004-637X/829/2/65

Huang, J., & MAPS Team. 2021, ApJ, 0, 0, doi: 0

Huang, J., Andrews, S. M., Dullemond, C. P., et al. 2020, ApJ, 891, 48, doi: 10.3847/1538-4357/ab711e

Hughes, A. M., Andrews, S. M., Espaillat, C., et al. 2009, ApJ, 698, 131, doi: 10.1088/0004-637X/698/1/131

Hunter, J. D. 2007, Computing in Science and Engineering, 9, 90, doi: 10.1109/MCSE.2007.55

Isella, A., Guidi, G., Testi, L., et al. 2016, PhRvL, 117, 251101, doi: 10.1103/PhysRevLett.117.251101

Jones, E., Oliphant, T., Peterson, P., et al. 2001-, SciPy: Open source scientific tools for Python. http://www.scipy.org/

Kama, M., Trapman, L., Fedele, D., et al. 2020, A&A, 634, A88, doi: 10.1051/0004-6361/201937124

Kamp, I., & Dullemond, C. P. 2004, ApJ, 615, 991, doi: 10.1086/424703

Krijt, S., Schwarz, K. R., Bergin, E. A., & Ciesla, F. J. 2018, ApJ, 864, 78, doi: 10.3847/1538-4357/aad69b

Law, C., & MAPS Team. 2021, ApJ, 0, 0, doi: 0

Linsky, J. L. 1998, SSRv, 84, 285

Liu, S.-F., Jin, S., Li, S., Isella, A., & Li, H. 2018, ApJ, 857, 87, doi: 10.3847/1538-4357/aab718

Long, F., Herczeg, G. J., Pascucci, I., et al. 2017, ApJ, 844, 99, doi: 10.3847/1538-4357/aa78fc

Lynden-Bell, D., & Pringle, J. E. 1974, MNRAS, 168, 603, doi: 10.1093/mnras/168.3.603

Macías, E., Espaillat, C. C., Ribas, Á., et al. 2018, ApJ, 865, 37, doi: 10.3847/1538-4357/aad811

Mathis, J. S., Rumpl, W., & Nordsieck, K. H. 1977, ApJ, 217, 425, doi: 10.1086/155591
McClure, M. K., Bergin, E. A., Cleeves, L. I., et al. 2016, ApJ, 831, 167, doi: 10.3847/0004-637X/831/2/167

McMullin, J. P., Waters, B., Schiebel, D., Young, W., & Golap, K. 2007, in Astronomical Society of the Pacific Conference Series, Vol. 376, Astronomical Data Analysis Software and Systems XVI, ed. R. A. Shaw, F. Hill, & D. J. Bell, 127

Miotello, A., Bruderer, S., & van Dishoeck, E. F. 2014, A&A, 572, A96, doi: 10.1051/0004-6361/201424712

Najita, J. R., & Ádámkovics, M. 2017, ApJ, 847, 6, doi: 10.3847/1538-4357/aa8632

Najita, J. R., & Bergin, E. A. 2018, ApJ, 864, 168, doi: 10.3847/1538-4357/aad80c

Öberg, K. I., Furuya, K., Loomis, R., et al. 2015, ApJ, 810, 112, doi: 10.1088/0004-637X/810/2/112

Öberg, K. I., & MAPS Team. 2021, ApJ, 0, 0, doi: 0

Pinte, C., Price, D. J., Ménard, F., et al. 2018, ApJ, 860, L13, doi: 10.3847/2041-8213/aac6dc

Powell, D., Murray-Clay, R., Pérez, L. M., Schlichting, H. E., & Rosenthal, M. 2019, ApJ, 878, 116, doi: 10.3847/1538-4357/ab20ce

Powell, D., Murray-Clay, R., & Schlichting, H. E. 2017, ApJ, 840, 93, doi: 10.3847/1538-4357/aa6d7c

Reboussin, L., Wakelam, V., Guilloteau, S., Hersant, F., & Dutrey, A. 2015, A&A, 579, A82, doi: 10.1051/0004-6361/201525885

Rowther, S., Meru, F., Kennedy, G. M., Nealon, R., & Pinte, C. 2020, ApJL, 904, L18, doi: 10.3847/2041-8213/abc704

Sakai, N., Oya, Y., Sakai, T., et al. 2014, ApJL, 791, L38, doi: 10.1088/2041-8205/791/2/L38

Schwarz, K. R., Bergin, E. A., Cleeves, L. I., et al. 2016, ApJ, 823, 91, doi: 10.3847/0004-637X/823/2/91
—. 2018, ApJ, 856, 85, doi: 10.3847/1538-4357/aaae08

Teague, R., Bae, J., Bergin, E. A., Birnstiel, T., & Foreman-Mackey, D. 2018, ApJ, 860, L12, doi: 10.3847/2041-8213/aac6d7

Tobin, J. J., Looney, L. W., Li, Z.-Y., et al. 2018, ApJ, 867, 43, doi: 10.3847/1538-4357/aae1f7

Toomre, A. 1964, ApJ, 139, 1217, doi: 10.1086/147861

Trapman, L., Facchini, S., Hogerheijde, M. R., van Dishoeck, E. F., & Bruderer, S. 2019, A&A, 629, A79, doi: 10.1051/0004-6361/201834723

Trapman, L., Miotello, A., Kama, M., van Dishoeck, E. F., & Bruderer, S. 2017, A&A, 605, A69, doi: 10.1051/0004-6361/201630308

van der Marel, N., Williams, J. P., Ansdell, M., et al. 2018, ApJ, 854, 177, doi: 10.3847/1538-4357/aaaa6b

van der Walt, S., Colbert, S. C., & Varoquaux, G. 2011, Computing in Science and Engineering, 13, 22, doi: 10.1109/MCSE.2011.37

Weaver, E., Isella, A., & Boehler, Y. 2018, ApJ, 853, 113, doi: 10.3847/1538-4357/aaa481

Williams, J. P., & Best, W. M. J. 2014, ApJ, 788, 59, doi: 10.1088/0004-637X/788/1/59

Wilson, T. L. 1999, Reports on Progress in Physics, 62, 143, doi: 10.1088/0034-4885/62/2/002

Woitke, P., Min, M., Pinte, C., et al. 2016, A&A, 586, A103, doi: 10.1051/0004-6361/201526538

Woitke, P., Kamp, I., Antonellini, S., et al. 2019, PASP, 131, 064301, doi: 10.1088/1538-3873/aaf4e5

Woodall, J., Agúndez, M., Markwick-Kemper, A. J., & Millar, T. J. 2007, A&A, 466, 1197, doi: 10.1051/0004-6361:20064981

Xu, R., Bai, X.-N., & Öberg, K. 2017, ApJ, 835, 162, doi: 10.3847/1538-4357/835/2/162

Yen, H.-W., Gu, P.-G., Hirano, N., et al. 2019, ApJ, 880, 69, doi: 10.3847/1538-4357/ab29f8

Yu, M., Willacy, K., Dodson-Robinson, S. E., Turner, N. J., & Evans, II, N. J. 2016, ApJ, 822, 53, doi: 10.3847/0004-637X/822/1/53

Zhang, & MAPS Team. 2021, ApJ, 0, 0, doi: 0

Zhang, K., Bergin, E. A., Blake, G. A., Cleeves, L. I., & Schwarz, K. R. 2017, Nature Astronomy, 1, 0130, doi: 10.1038/s41550-017-0130