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

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

$$\begin{eqnarray}&&{\rm{\Sigma }}(R)={{\rm{\Sigma }}}_{{\rm{c}}}{\left(\frac{R}{{R}_{{\rm{c}}}}\right)}^{-\gamma }\exp \left[-{\left(\frac{R}{{R}_{{\rm{c}}}}\right)}^{2-\gamma }\right],\end{eqnarray} \tag{ 1 }$$

where Σ_c/e is the surface density at a characteristic radius R_c. The 2D density distribution is then

$$\begin{eqnarray}&&\rho (R,Z)=\displaystyle \frac{{\rm{\Sigma }}}{\sqrt{2\pi }{Rh}}\exp \left[-\displaystyle \frac{1}{2}{\left(\displaystyle \frac{Z}{h}\right)}^{2}\right].\end{eqnarray} \tag{ 2 }$$

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

$$\begin{eqnarray}&&h(R)={h}_{\mathrm{ref}}{\left(\displaystyle \frac{R}{{R}_{\mathrm{ref}}}\right)}^{\psi }\end{eqnarray} \tag{ 3 }$$

where h_ref is the characteristic scale height at a radius R_ref and ψ 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., they have the same scale height, and follow the same surface density profile. The scale height of 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 et al. (2021).

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 et al. (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) ∝ a^−3.5 (Mathis et al. 1977) and use the standard dust opacities from Birnstiel et al. (2018), who assumed dust composed of 20% H₂O ice, 32.9% astronomical silicates, 7.4% troilite, and 39.7% refractory organics by mass. This differs from the 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 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.

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 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). The 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 photodesorption, 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 chemical 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 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 ¹³C, ¹⁸O, and deuterium. Instead, isotopologue emission profiles are generated assuming ¹³C/¹²C = 69, ¹⁶O/¹⁸O = 557, ¹⁸O/¹⁷O = 3.6, and D/H = 1.5 × 10⁻⁵ as measured in the local ISM (Linsky 1998; Wilson 1999).

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 density 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 is 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 et al. (2021), who modified 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 et al. (2021). Initially, we use the CO depletion profile of Zhang et al. (2021). We then run additional thermochemical models 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 χ² for the CO isotopologue 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 χ² at each radius. This model, using an updated CO abundance profile but otherwise using the same parameters as the best-fit model from Zhang et al. (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¹⁸O J = 2–1, ¹³CO J = 2–1, ¹²CO J = 2–1, and ¹³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 the models are extracted with the Python package disksurf, ²⁷ using the method presented by Pinte et al. (2018). In regions where the line flux is weak or originates 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 et al. (2021). Finally, we compare the results of each model to the total observed HD flux.

Zoom In Zoom Out Reset image size

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 midplane 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.

Zoom In Zoom Out Reset image size

Figure 8 compares the CO depletion profile after evolution of the chemistry in our best-fit model for 1 Myr to that found by Zhang et al. (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 midplane 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 et al. (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 et al. (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 at roughly 90–150 au, which is not seen to the same extent in the Zhang et al. (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 match a greater number of observations, some variation is to be expected.

Zoom In Zoom Out Reset image size

Figure 9 shows the 2D gas temperature distribution in a subset of our models with the derived temperatures for the ¹²COJ = 2–1, ¹³CO J = 2–1, and ¹³CO J = 3–2 surfaces overplotted. C¹⁸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 et al. (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 ¹³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).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

However, at larger heights our model is over two times warmer than the temperature measured from the ¹²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. 2021a). 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 residual cloud or envelope—the favored explanation for the nonaxisymmetric features seen in ¹²CO J = 2–1 (Huang et al. 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 ¹²CO emits. Additionally, the infalling material is unlikely to have undergone much chemical processing, and thus will have a CO/H₂ abundance 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₂ could result in some ¹²CO emission originating closer to the disk surface.

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.

Zoom In Zoom Out Reset image size

Previous analysis of HD detection by McClure et al. (2016) in GM Aur constrains the disk gas mass to 0.025–0.204 M_⊙. Based on analysis of the millimeter and centimeter continuum, Macías et al. (2018) found a total dust mass of 2 M_J. Assuming a gas-to-dust mass ratio of 100, this corresponds to a total gas mass of 0.19 M_⊙, 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 between 0.01 M_⊙ 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_⊙ and also reasonably reproduces the observed total flux of the 63 μm [OI] line as well as the ¹²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 ¹²CO J = 2–1 line, has a disk mass of 3.3 × 10⁻² M_⊙.

Our best-fit model has a total gas mass of 0.2 M_⊙, consistent with the upper limits of previous works. This 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 ¹²CO 2–1 emission, though less than 3% of the total gas mass is outside of 300 au.

The stability of a rotating disk against gravitational collapse is often characterized using the Toomre Q parameter (Toomre 1964):

$$\begin{eqnarray}&&Q=\frac{{c}_{{\rm{s}}}{\rm{\Omega }}}{\pi G{\rm{\Sigma }}}\end{eqnarray} \tag{ 4 }$$

where c_s is the gas sound speed assuming the midplane gas temperature, Ω is the Keplerian angular velocity of the disk, and Σ is the total gas+dust surface density. For a geometrically thin disk Q ∼ 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.

Zoom In Zoom Out Reset image size

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, 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 a Q of 1.7 while holding the temperature constant. Zhang et al. (2021) note that for four 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 et al. (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 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 by several kelvin as 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 et al. 2021a).

The dip in Toomre Q in our model is due entirely to an overdensity 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 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 development and stabilizing the disk (Rowther et al. 2020).

GI is thought to primarily manifest as nonaxisymmetric features. Previous analysis of the continuum disk emission at far-UV 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 intensity 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 artifact; instead the variation can be attributed to either a nonaxisymmetric or a 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). A 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, observations of HCO⁺ J = 3–2 toward GM Aur do not appear to deviate substantially from Keplerian rotation (Huang 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 a 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.