A Multifrequency ALMA Characterization of Substructures in the GM Aur Protoplanetary Disk

The 1.1 and 2.1 mm ALMA continuum images, as well as their corresponding azimuthally averaged radial intensity profiles, are shown in Figure 1. The radial profiles are computed by deprojecting the continuum images using P.A. = 5717 and i = 5321, then averaging the pixel intensities within annular bins one au wide. The continuum observations in both bands reveal faint, compact emission at the center of the disk, bright narrow rings at R ∼ 40 and ∼84 au, and faint, diffuse emission beyond R ∼ 100 au. At 1.1 mm, a faint ring is visible at R ∼ 168 au on top of the outer diffuse emission. While the 2.1 mm image does not show an unambiguous counterpart to this ring, perhaps due to the lower S/N of the data, there appears to be a slow rise in emission toward R ∼ 170 au and a steeper falloff outside this radius. In accordance with the nomenclature from Huang et al. (2018b), each ring is labeled with the prefix "B" (for "bright") followed by the radial location of the emission maximum rounded to the nearest whole number of astronomical units. The convention is similar for the gaps, except the prefix "D" (for "dark") is used. B40, D67, and B84, as well as the diffuse outer emission, were previously inferred from 930 μm observations at a resolution of ∼03 in Macías et al. (2018). D15 corresponds to the GM Aur disk's well-known central cavity (e.g., Hughes et al. 2009), although the cavity might be more precisely described as an annular gap given the detection of interior emission. For the sake of continuity with previous works on GM Aur, we refer to this feature as the "cavity" in the rest of the paper.

Zoom In Zoom Out Reset image size

The new observations, which improve upon the resolution of previous GM Aur observations by an order of magnitude, show that B40 and B84 are not radially symmetric. This characteristic is more apparent in radial profiles generated from averaging only the pixels within 20° of the projected disk major axis (Figure 2), especially since the synthesized beam is narrower along the disk major axis. At 1.1 mm, the emission profile of B40 is steeper on the side facing the star. At 2.1 mm, the emission profile of B40 appears to be narrower and more symmetric around the peak compared to the 1.1 mm image. Despite the slightly lower resolution of the 2.1 mm data, the appearance of an outer shoulder makes the two-component nature of B40 clearer. The differing emission profiles at the two wavelengths may either be due to the lower optical depth at 2.1 mm or to the 2.1 mm emission being more sensitive to larger grains, which may be confined to narrower regions by pressure bumps (Whipple 1972). The radial asymmetry of B84 is more subtle and is most easily highlighted by superimposing Gaussian profiles (derived by fitting the emission profiles between 80 and 90 au using the Levenberg–Marquardt minimization implementation in scipy.optimize.curve_fit). In both bands, the B84 emission profile is shallower on the side facing away from the star.

Zoom In Zoom Out Reset image size

The B40 emission has a modest azimuthal asymmetry at 1.1 mm. The southwest side is brighter than the northeast side, with a 0.073 ± 0.014 mJy beam⁻¹ (∼8%) difference in peak intensities. Macías et al. (2018) report a similar brightness asymmetry at the 5σ level in lower-resolution 930 μm data and at the 2σ level in 7 mm data. In our 2.1 mm data, the peak intensity of the northeast side is actually ∼5% brighter than the southwest side, but the difference is not statistically significant. The S/N of the 2.1 mm data is only about half that of the 1.1 mm data, making it more difficult to determine whether asymmetries are present.

Measuring intensities at different frequencies can constrain the optical depth and dust grain properties of disks (e.g., Beckwith & Sargent 1991). To compare the GM Aur disk emission in different bands, the CASA imsmooth task is used to smooth the 1.1 mm continuum image to the same resolution as the 2.1 mm image (57 mas × 34 mas (−132)). The normalized intensity profiles are shown in the top left of Figure 3. While the emission profiles are similar, there are several notable differences. At 2.1 mm, the emission rings are slightly narrower, the contrasts between the ring peaks and gap troughs are slightly larger, the central emission component is brighter relative to the peak intensity, and the relative brightness of the outer diffuse emission is lower.

Zoom In Zoom Out Reset image size

Changes in intensity as a function of frequency are quantified with the spectral index, $d\mathrm{log}{I}_{\nu }/d\mathrm{log}\nu$. It is typically assumed that the intensity scales as ${I}_{\nu }\propto {\nu }^{\alpha }$, in which case

$$\begin{eqnarray}&&\alpha =\displaystyle \frac{\mathrm{log}{I}_{{\nu }_{1}}/{I}_{{\nu }_{0}}}{\mathrm{log}{\nu }_{1}/{\nu }_{0}}.\end{eqnarray} \tag{ 1 }$$

Accounting for the ∼10% systematic flux calibration uncertainty in each band and assuming that the probability distribution is Gaussian, the disk-averaged spectral index between 1.1 and 2.1 mm is α = 2.7 ± 0.2. This is similar to the disk-averaged spectral index (2.94 ± 0.44) that Pinilla et al. (2014) measured for GM Aur between 880 μm and 3 mm.

There is significant radial variation in α, as shown in the bottom left of Figure 3. While the systematic flux calibration uncertainty leads to a large uncertainty (∼0.2, like that for the disk-averaged α) in the absolute offset of the radial α profile, the uncertainty of the profile shape is much smaller because it is determined by the S/N of the observations. Interior to R ∼ 100 au, local minima in α occur near the emission peaks, while local maxima in α occur near continuum gaps. The α variations suggest that either larger dust grains are segregated in the rings or that the rings are optically thick while the gaps are optically thin (e.g., Ricci et al. 2012; Tsukagoshi et al. 2016). We examine the extent to which we can distinguish between these possibilities in Section 4.

In the outer diffuse emission region, the α profile measured from the high-resolution images is noisy. To improve sensitivity, we re-image both bands of data by applying a Gaussian taper of 006 and then using imsmooth to smooth to a common resolution of 110 mas × 90 mas. The corresponding radial intensity and α profiles are shown on the right side of Figure 3. In the inner 100 au, the α variations are muted due to the degraded resolution, but in the outer disk, the improved sensitivity reveals modest radial variations in α. In contrast with the pattern established in the inner disk, the local extrema do not coincide with the locations of D145 and B168.

Interpreting the spectral index profile requires an estimate of the disk optical depth. To do this, we fit for GM Aur's surface brightness profile at each wavelength in the uv plane. To make fitting the data more computationally tractable, each SPW is first frequency-averaged down to a single channel. Since CASA measurement sets store uv coordinates in units of meters, we convert the coordinates to wavelength units (λ) using the frequencies of the individual SPWs. To reduce the data volume further, we follow the example of Hezaveh et al. (2016) for modeling long-baseline ALMA data and bin the visibilities for each ALMA band into 12 kλ × 12 kλ cells in the uv plane (i.e., comparable at 1.1 mm to the ALMA antenna diameter).

An axisymmetric model is adopted, given that the azimuthal variations described in Section 3 are modest. The surface brightness profile is parameterized such that the central emission component is represented by a delta function, B40 is modeled as the sum of two overlapping Gaussian rings to account for the "bump" observed in Band 4, B84 is modeled as an "asymmetric Gaussian" ring, B168 is modeled as a Gaussian ring, and the diffuse outer emission is modeled as a broad Gaussian that is truncated at r = 0. This profile can be written as

$$\begin{eqnarray}\begin{array}{rcl}{I}_{\nu }(r) & = & {A}_{0}{\displaystyle \int }_{0}^{\infty }\delta (r^{\prime} ){dr}^{\prime} \\ & & +{A}_{1a}\exp \left(-\displaystyle \frac{{\left(r-{r}_{1a}\right)}^{2}}{2{\sigma }_{1a}^{2}}\right)\\ & & +{A}_{1b}\exp \left(-\displaystyle \frac{{\left(r-{r}_{1a}-{\rm{\Delta }}r\right)}^{2}}{2{\sigma }_{1b}^{2}}\right)\\ & & +{A}_{2}\exp \left(-\displaystyle \frac{{\left(r-{r}_{2}\right)}^{2}}{2{\sigma }_{2}^{2}}\right)\\ & & +\sum _{i=3}^{i=4}{A}_{i}\exp \left(-\displaystyle \frac{{\left(r-{r}_{i}\right)}^{2}}{2{\sigma }_{i}^{2}}\right)\end{array}\end{eqnarray} \tag{ 2 }$$

where

$$\begin{eqnarray}\begin{array}{l}{\sigma }_{2}=\left\{\begin{array}{ll}{\sigma }_{2,\mathrm{in}} & r\leqslant {r}_{2}\\ {\sigma }_{2,\mathrm{out}} & r\gt {r}_{2}.\end{array}\right.\end{array}\end{eqnarray} \tag{ 3 }$$

This surface brightness prescription is partly motivated by Macías et al. (2018), which used three concentric Gaussian rings to model lower-resolution observations of GM Aur at 930 μm and 7 mm. The additional free parameters are the disk's P.A., cos i, R.A. offset from the phase center (δ_α), and decl. offset from the phase center (δ_δ), for a total of 21 free parameters. Positive offsets are defined to be north and east of the phase center, respectively.

A model disk image is first generated without the central point source. Using the Python package vis_sample (Loomis et al. 2018), synthetic visibilities ${{ \mathcal V }}_{m}$ are produced by sampling the model image at the same uv coordinates as the observations and performing a phase-shift. The central point source is directly added in the uv plane in the form of a constant A₀ phase-shifted by the model offset. Each model is compared to the observed visibilities ${{ \mathcal V }}_{d}$ using the log-likelihood function

$$\begin{eqnarray}&&\mathrm{log}p({{ \mathcal V }}_{d}| {\rm{\Theta }})=-\displaystyle \frac{1}{2}\sum _{i}\left({W}_{i}| {{ \mathcal V }}_{d,i}-{{ \mathcal V }}_{m,i}{| }^{2}+\mathrm{ln}\displaystyle \frac{2\pi }{{W}_{i}}\right),\end{eqnarray} \tag{ 4 }$$

where W_i is the weight corresponding to visibility ${{ \mathcal V }}_{i}$ and Θ are the model parameters. The weights used in the likelihood calculations are scaled down from the nominal data weights provided in the delivered measurement sets by a factor of 2.667 because CASA's weight averaging procedure during data binning does not account for the effective channel width introduced by Hanning smoothing. This scaling factor was checked by computing the scatter of visibilities close to one another in uv space.

Uniform priors are adopted for the parameters defining I_ν. The bounds of the priors were determined through a combination of considerations, including previous modeling results for GM Aur from Macías et al. (2018), the size of the synthesized beam (i.e., lower or upper bounds for the widths of sources can be set depending on whether the features are well-resolved), and from manual testing to check that the priors are not overly restrictive. Gaussian priors are selected for the parameters governing the disk orientation and phase offset. The priors for the P.A. and $\cos (i)$ are set to the best-fit values derived in Macías et al. (2018), while the standard deviations are set to be a few times wider than the posteriors from Macías et al. (2018) because additional substructures are being modeled and can affect the disk orientation measurements. The priors for δ_α and δ_δ are centered at 0, and each has a standard deviation of 002 (comparable to the scale of the synthesized beam). All priors are listed in Table 3.

Table 3.  Continuum Surface Brightness Model Parameters

Parameter	Prior Type	Band 6 Priora	Band 4 Priora	Band 6 Result	Band 4 Result	Units
A0	Uniform	[0, 2 × 10−4]	[0, 2 × 10−4]	1.1 × 10−4 ± 9 × 10−6	$8.9\times {10}_{-7\times {10}^{-6}}^{-5+6\times {10}^{-6}}$	Jy
A1a	Uniform	[0.3, 0.7]	[0.1, 0.25]	0.50 ± 0.02	${0.167}_{-0.004}^{+0.005}$	Jy arcsec−2
r1a	Uniform	[0.2, 0.25]	[0.2, 0.25]	0.2337 ± 0.0006	0.2361 ± 0.0003	arcsec
σ1a	Uniform	[0.02, 0.06]	[0.01, 0.03]	${0.0257}_{-0.0006}^{+0.0005}$	0.0160 ± 0.0009	arcsec
A1b	Uniform	[0.3, 0.7]	[0.05, 0.25]	0.456 ± 0.007	${0.131}_{-0.003}^{+0.002}$	Jy arcsec−2
Δr	Uniform	[0, 0.1]	[0, 0.1]	${0.0641}_{-0.0013}^{+0.0011}$	0.049 ± 0.002	arcsec
${\sigma }_{1b}$	Uniform	$[0.02,0.06]$	[0.02, 0.05]	${0.0444}_{-0.0008}^{+0.0009}$	${0.0441}_{-0.0012}^{+0.0009}$	arcsec
A2	Uniform	[0.25, 0.3]	[0.05, 0.1]	0.2778 ± 0.0008	0.0812 ± 0.0007	Jy arcsec−2
r2	Uniform	$[0.5,0.55]$	[0.5, 0.55]	0.5202 ± 0.0005	0.5222 ± 0.0012	arcsec
${\sigma }_{2,\mathrm{in}}$	Uniform	$[0.02,0.06]$	[0.02, 0.06]	0.0278 ± 0.0005	0.025 ± 0.001	arcsec
${\sigma }_{2,\mathrm{out}}$	Uniform	$[0.02,0.06]$	[0.02, 0.06]	0.0523 ± 0.0004	0.0449 ± 0.0009	arcsec
A3	Uniform	[0.025, 0.045]	[0, 0.01]	0.0384 ± 0.0002	${0.00597}_{-0.00009}^{+0.0001}$	Jy arcsec−2
r3	Uniform	$[0.4,0.9]$	[0.4, 0.9]	0.553 ± 0.008	0.68 ± 0.03	arcsec
${\sigma }_{3}$	Uniform	$[0.35,0.55]$	[0.35, 0.55]	0.495 ± 0.003	0.431 ± 0.012	arcsec
A4	Uniform	[0.005, 0.02]	[0, 0.005]	0.00954 ± 0.00014	0.00182 ± 0.00012	Jy arcsec−2
r4	Uniform	$[1,1.25]$	[1, 1.25]	1.114 ± 0.002	1.119 ± 0.006	arcsec
${\sigma }_{4}$	Uniform	$[0.05,0.2]$	[0.05, 0.25]	0.135 ± 0.004	0.097 ± 0.012	arcsec
i	Gaussian	[52.8, 1.5]	[52.8, 1.5]	53.20 ± 0.01	53.29 ± 0.03	degree
P.A.	Gaussian	[56.5, 2]	[56.5, 2]	57.16 ± 0.02	57.24 ± 0.05	degree
${\delta }_{x}$	Gaussian	[0, 0.02]	[0, 0.02]	$-{0.00172}_{-0.00004}^{+0.00005}$	0.0072 ± 0.0001	arcsec
${\delta }_{y}$	Gaussian	[0, 0.02]	[0, 0.02]	−0.00428 ± 0.00005	−0.0011 ± 0.0001	arcsec

Note.

^aIf the prior is uniform, the numbers in brackets denote the bounds of the prior. If the prior is Gaussian, the first number corresponds to the center of the Gaussian and the second number corresponds to the standard deviation.

Download table as:  ASCIITypeset image

The posterior probability distributions for the models at each wavelength are sampled with the emcee implementation of the affine invariant MCMC sampler (Goodman & Weare 2010; Foreman-Mackey et al. 2013). Each ensemble employs 96 walkers for 20,000 steps, with the first 10,000 steps discarded as burn-in. Convergence is checked by verifying that the chains are much longer than the estimated autocorrelation times (typically on the order of a few hundred). The median values of the posterior distribution are listed in Table 3, with error bars computed from the 16th and 84th percentiles.

The deprojected, azimuthally averaged visibilities corresponding to the best-fit models are compared to the observed visibilities in Figure 4. The models reproduce the real part of the visibilities well. Because the models are axisymmetric, the imaginary part of the deprojected and azimuthally averaged visibilities by construction should be zero (apart from numerical noise), but non-axisymmetric structure manifests in the data as non-zero imaginary components of the visibilities. The model and residual visibilities are then imaged with CASA in the same manner as the observations. A comparison of the model images and radial profiles to the data, as well as the residual images, are shown in Figure 5. Some significant residuals remain, as large as 8.4σ at 1.1 mm and 6.8σ at 2.1 mm. Part of the residuals are due to the non-axisymmetric structure around B40 discussed in Section 3. The residuals around B84 have a systematic appearance, with the 1.1 mm model overpredicting emission along the major axis and underpredicting emission along the minor axis. Explanations for the discrepancy at B84 are discussed further in Section 4.3. In addition, the B168 ring is not as pronounced in the 1.1 mm model as in the observations. However, the model radial profiles reproduce the observed intensities sufficiently well for the purpose of examining optical depths in Section 4.2.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

4.2.1. Constraints on Optical Depths

We now use our model intensity profiles to determine whether optically thick emission can account for the low spectral index values of GM Aur's dust rings. To assess the optical depths, we first plot the quantity

$$\begin{eqnarray}&&{\tau }_{\mathrm{nominal}}(r)=-\mathrm{ln}\left(1-\displaystyle \frac{{I}_{\nu }(r)}{{B}_{\nu }(r)}\right)\end{eqnarray} \tag{ 5 }$$

in Figure 6. I_ν is the surface brightness model at a given frequency ν, and B_ν is the Planck function evaluated at the midplane dust temperatures derived in Macías et al. (2018) through radiative transfer modeling of GM Aur's SED and resolved millimeter continuum observations. The expression for τ_nominal is typically used to estimate the optical depth in the limit where the dust is optically thin or the scattering opacities are small. The spectral index profile α computed from the best-fit surface brightness models is plotted in the same figure.

Zoom In Zoom Out Reset image size

The dominant source of uncertainty in τ_nominal is the midplane dust temperature. Unfortunately, the uncertainties associated with dust temperatures derived from radiative transfer modeling are usually ill-quantified due to the computational expense of exploring parameter space. However, the 1.1 mm continuum brightness temperatures set a lower bound on the possible midplane dust temperatures—the true midplane temperatures cannot be lower than the model temperatures by more than ∼35%. As shown later in Section 5, the brightness temperatures of optically thick HCO⁺, which emits from a warmer elevated layer, indicate that the true midplane temperatures cannot be more than a factor of two higher than the model temperatures.

With these uncertainties in mind, we can use Figure 6 to examine which parts of the disk are likely to be optically thick or thin. τ_nominal < 1 throughout the disk at both wavelengths. However, one cannot immediately conclude that the disk is completely optically thin—the peak values at B40 and B84 at 1.1 mm are high enough (τ_nominal = 0.94 and 0.70, respectively) such that assuming midplane temperatures a few degrees lower would push the values of τ_nominal above 1. Furthermore, Zhu et al. (2019) point out that when the dust albedo is sufficiently high, τ_nominal can be as low as ∼0.6 at millimeter wavelengths even in an optically thick disk.

On the other hand, τ_nominal ≪ 1 at D67 and beyond R ∼ 100 au. Even with large temperature uncertainties, these regions must be optically thin, and therefore, τ_nominal is a good approximation of the optical depth. The low optical depths are expected given the high spectral indices ($\alpha \gtrapprox 3$) in these regions. The gap at D67 is not completely evacuated—the best-fit model indicates τ_{1.1 mm} ∼ 0.06. Although the inner boundary of the temperature model stops short of the innermost gap, the intensities interior to 30 au are lower than the intensity at D67, which is presumably colder than the inner disk. Thus, the innermost gap should also be optically thin (with perhaps the exception of the central unresolved emission, which is discussed in Section 6.2).

4.2.2. Constraints on Dust Grain Sizes

For optically thin disk regions, dust grain sizes can be constrained by measuring how the optical depth changes between two frequencies. Because of the uncertainties associated with disk temperature and absolute flux calibration, the aim of this section is to comment on which regions of parameter space are consistent with the observations rather than to provide precise measurements.

The frequency dependence of the dust absorption opacity κ_abs is usually quantified with the β parameter, where κ_abs ∝ ν^β. Since ${\tau }_{\mathrm{abs}}={\kappa }_{\mathrm{abs}}{{\rm{\Sigma }}}_{d}/\mu$, where ${{\rm{\Sigma }}}_{d}$ is the dust surface density and $\mu =\cos i$, τ_abs has the same frequency dependence as κ_abs. Thus, we can use the τ_nominal values measured from the 1.1 and 2.1 mm surface brightness models to estimate β in optically thin parts of the disk. The model β radial profile (labeled β_nominal to signify that it is computed by neglecting scattering) is shown in Figure 6.

To connect β values to grain sizes, we adopt the default "DSHARP dust opacities" described in Birnstiel et al. (2018) and use the companion Python package dsharp_opac (Birnstiel 2018) to compute quantities derived from the opacities. The opacities are based on optical constants from Henning & Stognienko (1996), Draine (2003), and Warren & Brandt (2008). Throughout this paper, we assume the dust population follows a power-law size distribution $n(a)\propto {a}^{-p}$ and fix the minimum grain size at 0.1 μm. The specific choice of a_min does not have a large effect on the millimeter wavelength opacities as long as it is much smaller than a_max (e.g., Draine 2006). Figure 7 shows how β varies as a function of the maximum dust grain size a_max for three different power laws: p = 3.5 (i.e., the standard ISM value from Mathis et al. 1977), p = 2.5, and p = 1.5. The shallower power-law distributions may arise via grain growth (e.g., Weidenschilling et al. 1997). The β values measured for GM Aur at D67, D145, and B168 are shaded in gray. They are best matched by a_max values from 1 to 3 mm, with D67 also matching well with a_max values smaller than ∼100 μm due to the non-monotonic shape of β as a function of a_max. However, the 10% absolute flux calibration uncertainty leads to an absolute offset uncertainty of ∼0.3 in GM Aur's β profile, so the β values at these locations are consistent with a large range of grain sizes. Nevertheless, the upward trend in β beyond 200 au indicates that it is unlikely that a_max is significantly less than 1 mm in the optically thin outer disk (R ≳ 100 au). The very high β values near a_max ∼ 0.5 mm are inconsistent with the observations. For values of a_max ≲ 0.5 mm, β is either flat or an increasing function of a_max. In this regime, GM Aur's β profile in the outer disk would imply that a_max is increasing with distance from the star, which would be difficult to reconcile with standard models indicating that larger dust grains preferentially drift inward (e.g., Weidenschilling 1977).

Zoom In Zoom Out Reset image size

The estimated peak optical depths of B40 and B84 are high enough that using β_nominal to constrain grain sizes can lead to significant over-estimates (e.g., Carrasco-González et al. 2019). Instead, the effects of scattering should explicitly be considered. To compute the emergent intensity ${I}_{\nu }^{\mathrm{out}}$, we use the analytic approximations from Sierra et al. (2019) and Carrasco-González et al. (2019), which are summarized in the Appendix. Similar results are obtained using formulae based on the Eddington–Barbier approximation from Birnstiel et al. (2018) and Zhu et al. (2019). As noted in the Appendix, ${I}_{\nu }^{\mathrm{out}}$ depends on five parameters: p, a_max, the dust temperature T_d, Σ_d, and $\mu =\cos i$. The disk inclination is known, and we can use the midplane dust temperatures derived in Macías et al. (2018). This still leaves three free parameters with only two constraints (i.e., the intensities at each wavelength), so we cannot solve for the values of these parameters. Although GM Aur has been observed at other wavelengths (e.g., Macías et al. 2018), the much coarser spatial resolution of earlier observations prevents us from accurately estimating the intensities at the ring peaks.

We can, however, examine how accounting for scattering affects inferences about a_max and Σ_d for several possible values of p. For each of p = 3.5, 2.5, and 1.5, we compute the emergent intensities at 2.1 and 1.1 mm for a grid of Σ_d and a_max values at the temperatures corresponding to the peaks of B40 and B84. One set of calculations ("Absorption-only") is performed by setting the scattering opacities to zero everywhere, while "Absorption and scattering" uses the DSHARP scattering opacities. To account for the flux calibration uncertainty, Figure 8 shades in the combinations of a_max and Σ_d that produce intensities within 10% of the best-fit model intensities at each wavelength at the peak of B40. Figure 9 does the same for B84. Overlapping shaded regions indicate which combinations of a_max and Σ_d are consistent with both bands of the GM Aur observations (leaving aside temperature uncertainties, since our focus is on the qualitative effect of neglecting scattering opacities).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Comparing the two rows of panels in Figure 8 shows that accounting for scattering allows for steeper p values than would be inferred from absorption-only calculations. Both Figures 8 and 9 also demonstrate that accounting for scattering yields solutions for a_max that can be smaller and for Σ_d that can be larger than the solutions derived from absorption-only calculations. However, it is also important to note that the solution space can be discontinuous, and factoring in scattering does not necessarily rule out optically thin dust or grains larger than a millimeter at the ring peaks. For some values of p, the solution space becomes discontinuous somewhere between a_max ∼ 0.01 to 0.1 cm (i.e., at values comparable to the wavelength of the observations) because the albedos are very high for these size distributions, so optically thick emission saturates below the measured intensities. This effect is explained in detail in Zhu et al. (2019). Thus, given the available data, it is ambiguous the extent to which trapping of large dust grains contributes to the low spectral index values measured at B40 and B84.

Despite the aforementioned ambiguities, we argue that GM Aur's millimeter continuum emission is likely tracing some degree of radial variation in dust properties inside R ∼ 100 au. If at least one of the rings is optically thin, then the large spectral index variations across the ring(s) must be a consequence of dust opacity variations due to changing grain sizes or compositions (e.g., Beckwith & Sargent 1991; Testi et al. 2003). If both rings are optically thick and the dust properties are radially uniform, then B84 would have a slightly lower spectral index than B40 due to the lower temperature at a larger radius (note that spectral indices are temperature dependent outside the Rayleigh–Jeans regime). Instead, the spectral index of B84 is higher. Thus, even if both rings are optically thick, the spectral index indicates that grains in B40 have different properties from the grains in B84. On the other hand, since the outer disk (R ≳ 100 au) is optically thin, the spectral index variations point to radial variations in dust properties. However, the local minima and maxima in the spectral index profile do not coincide neatly with D145 and B168, so it is not clear whether B168 is a dust trap.

Despite the large range of a_max values consistent with GM Aur's ring intensities, Figures 8 and 9 also show that for the largest and smallest a_max solutions, the corresponding Σ_d solutions are well above the surface density limit for which the Toomre Q parameter is equal to 1, assuming the standard gas-to-dust ratio of 100 and a stellar mass of 1.32 M_⊙ (Andrews et al. 2018b). The continuum emission does not exhibit any clear signatures of gravitational instability, such as spiral arms, suggesting that Q is greater than 1 at the emission rings. Thus, the gas-to-dust ratios at the rings would have to be below 100 if a_max is smaller than ∼1 mm or larger than a few millimeters.

Of course, as illustrated in Birnstiel et al. (2018), different assumptions about the dust grain composition and structure will affect how multiwavelength intensities are translated to grain sizes. For example, unlike the compact DSHARP grains, the β curves of highly porous grains are nearly flat as a function of a_max (e.g., Kataoka et al. 2014). In this case, the strong spectral index variations in the GM Aur disk indicate that the dust cannot simultaneously be optically thin and highly porous.

Our surface brightness model setup is appropriate for axisymmetric emission originating from a geometrically thin layer. The residuals remaining in Figure 5 suggest that at least one of these assumptions is not wholly appropriate. Examining the breakdown of these assumptions can provide useful clues into disk structure and evolution. The vertical distribution of dust traces how readily large grains settle to the midplane (e.g., Pinte et al. 2016), while departures from axisymmetry may point to the presence of a perturbing body (e.g., Lubow & Ogilvie 2001).

The 1.1 mm residuals near B84 are particularly interesting due to their systematic appearance. As noted in Section 4.1, the surface brightness model overpredicts emission along the projected disk major axis and underpredicts along the minor axis, suggesting that the inner edge of B84 is slightly more elongated along the major axis compared to the outer edge. This is more readily seen by deprojecting the observations, replotting the continuum as a function of radius and azimuthal angle, and drawing the 15σ contours to highlight the edges of B40 and B84 (Figure 10). The major axis of the original image runs between 90° and −90°, while the minor axis runs between 0° and 180/−180°. The contours of B40 and the outer contour of B84 appear to be approximately vertical, indicating that they are at constant radius (i.e., they are axisymmetric). However, the inner contour of B84 undulates, tracing larger radii at azimuthal angles corresponding to the major axis than at angles corresponding to the minor axis. In other words, the inner edge of B84 appears to be at a higher inclination than the outer edge of the ring. (Choosing slightly higher or lower contour levels yields the same effect.) For comparison, Figure 10 also shows a polar plot for the best-fit surface brightness model. All of the contours appear to be at a constant radius, showing that the distorted inner edge of B84 is not a consequence of the uv sampling or imaging artifacts.

Zoom In Zoom Out Reset image size

To explore possible explanations for the geometry of B84, we use RADMC-3D (Dullemond 2012) to perform 3D Monte Carlo radiative transfer calculations and generate synthetic observations of three parametric ring models. As a reference, we first generate an optically and geometrically thin disk model ("Flat RT Model"). Since annular substructures often appear to have different inclinations in scattered light observations due to the projection of the optically thick, flared disk surface (e.g., de Boer et al. 2016), we next generate an optically and geometrically thick model ("Thick RT Model") to determine whether a similar effect could be at play for GM Aur's millimeter continuum emission. Finally, we generate a mildly warped model ("Warped RT model") based on previous hypotheses that the GM Aur disk is warped (e.g., Dutrey et al. 2008; Hughes et al. 2009).

The dust surface density radial profile is modeled as an asymmetric Gaussian ring, analogous to the surface brightness model profile for the B84 ring:

$$\begin{eqnarray}\begin{array}{l}{{\rm{\Sigma }}}_{d}(r)=\left\{\begin{array}{ll}C\exp \left(-\displaystyle \frac{{\left(r-84\mathrm{au}\right)}^{2}}{2{w}_{\mathrm{in}}^{2}}\right) & r\leqslant 84\,\mathrm{au}\\ C\exp \left(-\displaystyle \frac{{\left(r-84\mathrm{au}\right)}^{2}}{2{w}_{\mathrm{out}}^{2}}\right) & r\gt 84\,\mathrm{au}.\end{array}\right.\end{array}\end{eqnarray} \tag{ 6 }$$

The temperature is parameterized as

$$\begin{eqnarray}&&T(r)={T}_{10}{\left(\displaystyle \frac{r}{10\mathrm{au}}\right)}^{-q}.\end{eqnarray} \tag{ 7 }$$

We adopt a vertically isothermal temperature structure because the millimeter continuum emission presumably originates from large dust grains settled in the midplane, and therefore the temperature variation should not be large within this layer. The vertical distribution of the dust is Gaussian, so the dust density is given by

$$\begin{eqnarray}&&{\rho }_{d}(r,z)=\displaystyle \frac{{{\rm{\Sigma }}}_{d}(r)}{\sqrt{2\pi }h(r)}\exp \left(-\displaystyle \frac{{z}^{2}}{2h{\left(r\right)}^{2}}\right),\end{eqnarray} \tag{ 8 }$$

where h(r) is the dust scale height. Due to vertical settling, the dust scale height is assumed to be some constant fraction of the gas pressure scale height $H(r)=\sqrt{\tfrac{{k}_{B}T(r){r}^{3}}{{\mu }_{\mathrm{gas}}{m}_{{\rm{H}}}{{GM}}_{* }}}$, where the mean molecular weight is μm_H = 2.37× the mass of atomic hydrogen. We use a stellar mass of 1.32 M_⊙ computed from stellar evolutionary tracks (Andrews et al. 2018b).

We use a dust population with p = 3.5 and a_max = 1 mm, which corresponds to an absorption opacity of κ_abs = 2.4 cm² g⁻¹ and scattering opacity of κ_sca = 20.6 cm² g⁻¹ at 1.1 mm. A more realistic "two-population" dust model would also include a population of dust grains with a sub-micron a_max in the upper layers of the disk (e.g., D'Alessio et al. 2006), but the grains in the upper layer are not expected to contribute significantly to the millimeter continuum emission because they only constitute a small fraction of the solid mass. The small grains in the upper layers are important for self-consistent thermal structure calculations, but we parameterize our temperature structures because ALMA observations only constrain the properties of dust in the midplane. We note that Section 4.2 suggests that the dust properties in the GM Aur disk could be quite different from what we use in this section. Dust opacities, temperatures, and surface densities are highly degenerate when modeling continuum emission. Thus, the models we present should be taken as illustrative, not as a "best fit" to GM Aur.

For each model, we compute the emission at 1.1 mm using 512 radial cells spaced logarithmically from 0.5 to 175 au, 1024 poloidal cells spaced evenly between $\tfrac{\pi }{6}$ and $\tfrac{5\pi }{6}$ (the midplane is at $\tfrac{\pi }{2}$), 64 azimuthal cells spaced evenly between 0 and 2π, and 10⁸ photon packages. Anisotropic scattering is treated using Müller matrices calculated with the Draine version¹⁰ of the Mie code by Bohren & Huffman (1983). The phase offset is fixed to the best-fit values derived from the 1.1 mm surface brightness modeling, while the overall disk inclination and P.A. are fixed to the weighted averages of the best-fit 1.1 and 2.1 mm models (5321 and 5717, respectively). Synthetic visibilities are generated from the radiative transfer output using vis_sample, and the resulting visibilities are imaged with CASA. The model parameters are listed in Table 4.

Table 4.  Continuum Radiative Transfer Model Parameters

Parameter	Flat	Thick	Warped
T10 (K)	58.1	45	58.1
q	0.525	0.525	0.525
h(r) (au)	0.2H(r)	0.4H(r)	0.2H(r)
C (g cm−2)	0.175	1.5	0.175
win (au)	5	3	3.5
wout (au)	7.5	4.5	7.5
${\rm{\Delta }}{i}_{\max }$ (°)	⋯	⋯	5

Download table as:  ASCIITypeset image

The "Flat RT Model" is relatively settled: h(r) = 0.2H(r). The temperature structure is set by a power-law fit to the GM Aur midplane temperature calculated by Macías et al. (2018). At the peak of the ring, the dust aspect ratio is h/r  ∼ 0.014. The parameters to generate the disk surface density are adjusted manually until the width and height of the ring in the CASA image of the radiative transfer model are comparable to those of B84. As shown in Figure 11(a), the contours trace approximately constant radii in the polar plot, as expected for a geometrically and optically thin disk. Doubling the dust scale height of the "Flat RT" (optically thin) model does not appreciably change the emission geometry of B84.

Zoom In Zoom Out Reset image size

The emitting surface of the "Thick RT Model" is elevated by changing the dust scale height to h(r) = 0.4H(r) and increasing the surface density relative to the "Flat RT Model" such that B84 becomes optically thick. Because increasing the surface density also increases the intensity, we compensate by decreasing T₁₀ and the width of the ring until the width and height of the ring in the CASA image of the radiative transfer model are comparable to those of B84 in GM Aur. At the peak of the ring, the dust disk aspect ratio is h/r ∼ 0.024. Figure 11(b) shows that similarly to the observations of B84, the inner contour of the radiative transfer model traces larger radii near azimuthal angles corresponding to the major axis of the disk image than near angles corresponding to the minor axis. Thus, the apparent inclination of the inner edge is higher than the true inclination. Meanwhile, the outer contour also undulates, but in the opposite direction, indicating that the apparent inclination of the outer edge is lower than the true inclination. Figure 11(c) demonstrates that by deprojecting the same model image with an inclination of 52° rather than the true inclination of 5321, one can obtain a polar plot where the outer contour appears straight on the polar plot but the inner contour undulates, similar to the GM Aur observations.

To produce the "Warped RT model," we start with the parameters from the "Flat RT model" and rotate dust annuli out of the disk plane by some angle Δi around the projected disk major axis. The rotation angle is parameterized as

$$\begin{eqnarray}\begin{array}{l}{\rm{\Delta }}i(r)=\left\{\begin{array}{ll}{\rm{\Delta }}{i}_{\max }\left(1-\displaystyle \frac{(r-75\ \mathrm{au})}{10\ \mathrm{au}}\right) & r\in (75\ \mathrm{au},85\ \mathrm{au})\\ 0 & \mathrm{everywhere}\ \mathrm{else}.\end{array}\right.\end{array}\end{eqnarray} \tag{ 9 }$$

In this parameterization, the inner edge of B84 is mildly misaligned with the plane of the disk, while the outer edge is coplanar. The radiative transfer models are generated such that the inner edge appears more highly inclined than the overall disk from the viewer's vantage point. The parameters are again adjusted until the width of the ring in the CASA image of the radiative transfer model is comparable to that of B84 in GM Aur. As shown in Figure 11, a modest value of i_max = 5° can yield an undulating inner contour and a straight outer contour, mimicking the behavior of the B84 ring.

The radiative transfer models indicate that non-axisymmetric structure and vertical structure are both plausible explanations for the emission geometry of B84. Better sensitivity and resolution could help to distinguish between these and other scenarios. Along the projected minor axis of the disk, the optically and geometrically thick radiative transfer model produces emission that is ∼10% brighter on the far side (southeast) of the B84 ring compared to the near (northwest) side because the warmer inner rim of the far side is more exposed to the viewer. The S/N of the B84 emission in the current GM Aur observations is not high enough to confirm whether such a difference is present (note that the ∼8% brightness difference for B40 along the projected major axis, as described in Section 3, is statistically significant because the inner ring is much brighter). Meanwhile, although our warped model is set up such that the inner edge of B84 is tilted around the projected disk major axis for the sake of simplicity, it is unlikely that GM Aur would have such a coincidental orientation. The case for a warp could be strengthened if better-quality observations also demonstrated that the P.A. of the inner edge of B84 differs from that of the outer edge. High spectral and spatial resolution line observations can also be used to test for the presence of a warp (e.g., Rosenfeld et al. 2014), although warps do not always leave a clear imprint on observed disk kinematics (e.g., Juhász & Facchini 2017). Thus, hydrodynamical simulations with more realistic warp geometries will also be needed to determine whether the observations are compatible with a disk warp. It should also be noted that if the radial thermal profile is more complex than assumed due to additional heating or cooling in the gaps (e.g., Facchini et al. 2018; van der Marel et al. 2018), fully self-consistent radiative transfer modeling may be needed to explain B84's emission geometry.

We use the Python packages bettermoments and eddy (Teague 2019; Teague & Foreman-Mackey 2018) to estimate the height of the HCO⁺ emission layer. The line-of-sight velocity corresponding to the emission line peak at each pixel is computed using a quadratic fit to the intensities as a function of LSRK velocity (Figure 13). This is functionally similar to an intensity-weighted velocity map (also known as a moment 1 map), but without emission from the lower surface biasing the estimate of the velocities of the brighter upper surface. The resulting map is then downsampled by 10 pixels (01) on each side such that there is approximately one pixel remaining for each beam-sized patch. The height of the emitting layer as a function of disk radius is parameterized as

$$\begin{eqnarray}&&z(r)={z}_{0}{\left(\displaystyle \frac{r}{1^{\prime\prime} }\right)}^{\psi }.\end{eqnarray} \tag{ 10 }$$

Assuming Keplerian kinematics, the quantities needed to compute the line-of-sight velocities for a given disk geometry are z₀, ψ, the stellar mass M_*, the systemic velocity v_LSR, the P.A., inclination, and offset from the phase center. The P.A., inclination, and phase offsets are fixed to the values derived from the continuum visibility modeling described in Section 4.1, since the S/N is higher for the continuum data. Broad uniform priors are adopted for the four free parameters and are listed in Table 5.

Zoom In Zoom Out Reset image size

Table 5.  HCO⁺ Emission Height Model

Parameter	Prior	Results	Units
z0	[0, 5]	0.127 ± 0.002	Arcseconds
ψ	[0, 5]	0.81 ± 0.03	Dimensionless
M*	[0.1, 2]	1.206 ± 0.004	${M}_{\odot }$
${v}_{\mathrm{LSR}}$	[0, 12]	5.612 ± 0.003	km s−1

Download table as:  ASCIITypeset image

The eddy backend uses emcee to fit the observed line-of-sight velocity map with the model Keplerian velocity map. The region of the fit is restricted to radii extending from 02 to 2'', where the inner radius is set by angular resolution limitations and the outer radius is set by S/N limitations and to avoid confusion from the dimmer back side of the disk. The posterior probability distributions are sampled with 48 walkers for 3500 steps, with the first 500 steps discarded as burn-in. Convergence is checked by estimating the autocorrelation time for each parameter, which is ∼40 steps. The posterior medians are listed in Table 5, with error bars calculated from the 16th and 84th percentiles. As noted in Keppler et al. (2019), the nominal error bars should be regarded with caution because there may also be systematic uncertainties associated with image-plane fitting and differences between the assumed model structure and true emission behavior.

The best-fit model adequately reproduces the geometry of the HCO⁺ emission, as shown in the comparison of the predicted isovelocity contours to channels tracing the front and back sides of the disk in Figure 14. The emission geometry demonstrates that the northwest side of the disk is tilted toward the observer, consistent with conclusions drawn from scattered light observations (Schneider et al. 2003). While Dutrey et al. (2008) and Hughes et al. (2009) identify a prominent discrepancy between the position angles of the continuum and ¹²CO emission, no such discrepancy is apparent for HCO⁺. This may be a consequence of the CO and HCO⁺ emission coming from different heights in the disk. Interpretation of GM Aur's ¹²CO emission has also been complicated by absorption from either intervening cloud or residual envelope material (e.g., Hughes et al. 2009). No cloud contamination or envelope material is seen in HCO⁺ emission. With a high critical density of ∼10⁶ cm⁻³ (Shirley 2015), the HCO⁺ J = 3 − 2 line should preferentially trace the denser disk material. Based on visual inspection of the channel maps and line-of-sight velocities map, we do not identify clear signatures of a velocity "twist" associated with warps and radial inflows (e.g., Rosenfeld et al. 2014) or a spectrally and spatially localized perturbation from Keplerian velocities (Pinte et al. 2018b). However, since the detection criteria for such features have not been established in detail, our statement on the presence or absence of these features is not intended to be definitive.

Zoom In Zoom Out Reset image size

The isovelocity curves corresponding to the best-fit model are used to generate a Keplerian mask (e.g., Rosenfeld et al. 2013; Salinas et al. 2017) to apply to the image cube before integrating along the spectral axis to produce an integrated intensity map. The integration range, −1.25 to 12.75 km s⁻¹, is selected to encompass channels with disk emission above the 3σ level. The flux is measured by summing all of the unmasked pixels. The flux uncertainty is estimated by generating 50 image cubes consisting only of noise by randomly drawing line-free channels and applying random position shifts, measuring the flux in each cube with the Keplerian mask applied, then taking the standard deviation of these flux measurements. The integrated intensity map and corresponding deprojected, azimuthally averaged radial profile are shown in Figure 15. The HCO⁺ integrated emission peaks at ∼20 au, which lies well interior to the peak of the continuum emission ring B40. The small central cavity in HCO⁺ emission is most likely due to a gas surface density reduction, given that much of the dust is cleared in the inner regions of the disk. Dutrey et al. (2008) inferred a similarly sized hole from modeling low-resolution CO isotopologue emission. The radial profile of the integrated intensity map shows emission decrements near the locations of the continuum rings B40 and B84. These features appear to be largely if not entirely an artifact of continuum subtraction from optically thick line emission (e.g., Boehler et al. 2017). The unresolved inner disk continuum emission, however, is neither extended enough nor bright enough to be responsible for the central HCO⁺ cavity.

Zoom In Zoom Out Reset image size

Spatially resolved, optically thick lines are useful for constraining disk gas temperatures because they require a minimal number of assumptions about disk properties (e.g., Pinte et al. 2018a). Among other things, accurate temperatures are fundamental for measuring disk masses, molecular abundances, and the properties of embedded protoplanets (e.g., Trapman et al. 2017; van der Marel et al. 2018). The brightness temperature of the optically thick HCO⁺ emission in GM Aur can be used to estimate the gas temperature at the best-fit emission height given in Table 5. Using the imaging procedure outlined in Section 2 and the median-stacking procedure described in Huang et al. (2018a), we produce a peak brightness temperature map of HCO⁺ without continuum subtraction in order to avoid an artificial reduction in the peak line intensities (e.g., Boehler et al. 2017; Weaver et al. 2018). The full Planck equation is used to convert peak intensities to brightness temperatures. A deprojected, azimuthally averaged radial profile of the brightness temperature map is produced using the disk coordinates calculated by eddy when fitting the flared HCO⁺ surface, since the peak intensity at each pixel emerges from well above the midplane. (Unless specified, other deprojections in this paper use the geometrically thin approximation.) The brightness temperature map and corresponding radial profile are shown in Figure 15. The profile between R ∼ 30–350 au is approximated well with a power law. Using Levenberg–Marquardt minimization to fit a power law $T(r)={T}_{100}{\left(\tfrac{r}{100\mathrm{au}}\right)}^{-q}$ to the brightness temperature profile and sampling at radii spaced approximately one synthesized beam apart from R = 30–350 au, we obtain best-fit results and 1σ errors of T₁₀₀ = 26.9 ± 3 and q = 0.43 ± 0.01. The absolute flux calibration uncertainty contributes an additional 10% uncertainty on T₁₀₀, although it does not change the shape of the curve.

We can compare the gas temperature estimated from the optically thick HCO⁺ emitting layer to GM Aur disk temperature models. Using a linear interpolation of the two-dimensional temperature structures from McClure et al. (2016) and Macías et al. (2018), we plot the temperatures corresponding to the estimated height of HCO⁺ alongside the brightness temperature profile in Figure 15. The McClure et al. (2016) temperatures are a good match to the HCO⁺ brightness temperature within ∼100 au but increasingly deviate outside 100 au. The Macías et al. (2018) temperatures are ∼30%–40% lower than the HCO⁺ brightness temperatures, a discrepancy that is too large to be explained by the typically quoted ALMA flux calibration uncertainties (∼10%). The key difference between the models is that the McClure et al. (2016) model invokes a much higher depletion of small dust grains in the disk upper layers compared to Macías et al. (2018). The Macías et al. (2018) model is constrained by the SED and by spatially resolved millimeter and centimeter observations, which are sensitive to the inner disk and the disk midplane rather than the intermediate disk layers traced by HCO⁺. Meanwhile, the McClure et al. (2016) model is constrained by the SED and HD line emission, which originates from intermediate disk layers. Thus, the McClure et al. (2016) model has a steep vertical temperature gradient in order to make the intermediate disk layers sufficiently warm to match the HD flux, while the Macías et al. (2018) temperature structure is nearly vertically isothermal from the midplane to z/r ∼ 0.2. (The models are also calculated at slightly different distances, 140 pc for McClure et al. 2016 versus 160 pc for Macías et al. 2018, but this does not account for the large difference in temperature structures.) The difference between these two models highlights the importance of developing an independent method of inferring temperature. Although a single HCO⁺ transition cannot be used to determine the shape of the vertical temperature gradient, this exercise illustrates the potential for constraining GM Aur's thermal structure through high-resolution observations of multiple optically thick lines probing different heights.

It should be noted that the temperature differences between the two models are not as large at the midplane. At a given radius, the midplane temperature from McClure et al. (2016) is typically a few degrees cooler than that of Macías et al. (2018; i.e., the opposite of the behavior at the estimated height of the HCO⁺ layer). Because continuum emission is, in principle, more sensitive to the midplane temperature than HD emission, we used the Macías et al. (2018) thermal structure to analyze the dust properties in Section 4.2. However, had we used the McClure et al. (2016) thermal structure instead, we would still conclude that the extent to which the spectral index profiles signify that the rings are dust traps is ambiguous. The dust optical depth would increase everywhere, but the gaps and the outer disk (R ≳ 100 au) would still be optically thin. Because β is only weakly dependent on temperature, the dust grain size constraints are similar for these optically thin regions. With the increased optical depths at the peaks of the emission rings, it remains possible to produce low spectral indices with relatively small a_max values (e.g., much smaller than 1 mm), but the presence of grains larger than a millimeter is not ruled out. Surface densities remain uncertain due to the possibility that the emission is saturated.

As one of the earliest transition disks to be characterized in detail, the GM Aur disk has long been hypothesized to host at least one giant protoplanet (e.g., Marsh & Mahoney 1992; Rice et al. 2003). To date, no directly imaged protoplanet candidate has been reported for GM Aur, nor for most other disks. However, the interpretation of gaps in protoplanetary disks as being due to planet–disk interactions has been encouraged by the close match between hydrodynamical simulations and observations (e.g., Bae et al. 2017; Zhang et al. 2018), protoplanet detections in the PDS 70 disk cavity (Keppler et al. 2018; Haffert et al. 2019), and identification of kinematic perturbations in the vicinity of disk gaps (e.g., Teague et al. 2018; Pinte et al. 2019).

Using the relationships between gap widths and planet masses derived in Zhang et al. (2018) from a hydrodynamical parameter space study of non-migrating protoplanets in disks, we can perform a back-of-the-envelope estimate of the mass of a planet that might create GM Aur's D67 gap. The relationships are calculated from models at 1.27 mm, so we use the 1.1 mm map of GM Aur to calculate the gap width because the intensity profiles should not significantly differ over this small wavelength range. We use the observations directly rather than our surface brightness model to compute the width, because the Zhang et al. (2018) relationships are given for intensity maps convolved with Gaussian kernels comparable in size to the GM Aur synthesized beam. As defined in Zhang et al. (2018), the normalized width of D67 is Δ = 0.23. Based on the SED and millimeter continuum observations, Macías et al. (2018) estimate a gas surface density of ∼23 g cm⁻² and a midplane temperature of ∼20 K at a radius of 67 au. Assuming hydrostatic equilibrium and adopting a stellar mass of 1.32 M_⊙ (Andrews et al. 2018b), this would correspond to h/r ∼ 0.06 at 67 au. If the maximum grain size is ∼1 mm, a planet with a mass between ∼0.1 and 0.4 M_J would be required to open the gap, assuming that the α viscosity parameter ranges between 10⁻⁴ and 10⁻².

Photoevaporation is another mechanism that can clear the inner regions of disks (e.g., Clarke et al. 2001; Ercolano et al. 2011; Owen et al. 2012). Traditionally, it has not been regarded as a mechanism likely to create GM Aur's inner cavity, given its large radius and the high accretion rate onto the star (e.g., Espaillat et al. 2010). More recently, Ercolano et al. (2018) and Wölfer et al. (2019) have shown that for disks depleted in gas-phase carbon and oxygen, X-ray photoevaporation can open gaps at tens of astronomical units, compared to ∼1 au at standard metallicity. The inner disks, and, therefore, high accretion rates, can be sustained longer when the disk metallicity is decreased.

Annular substructures have also been hypothesized to trace the locations of molecular snowlines because the freezeout of different volatiles is expected to modify the fragmentation and coagulation properties of dust grains (e.g., Zhang et al. 2015; Okuzumi et al. 2016; Pinilla et al. 2017). Qi et al. (2019) found that the inner and outer edges of N₂H⁺ emission in the GM Aur disk coincided with the continuum rings B40 and B84. The inner and outer boundaries of N₂H⁺ were hypothesized to be set by the CO and N₂ snowlines, respectively. An alternative interpretation of the N₂H⁺ observations toward GM Aur could be that ionization is enhanced inside large disk gaps, which was previously invoked by Favre et al. (2019) to explain the rise of DCO⁺ emission inside one of AS 209's disk gaps. This scenario needs to be tested by source-specific thermochemical models, since the emission geometry of N₂H⁺ is highly sensitive to both ionization and the thermal structure (e.g., Qi et al. 2019). While the N₂H⁺ distribution may suggest an association between the dust substructures and snowlines, the appearance of the GM Aur disk in scattered light does not. Models from Pinilla et al. (2017) indicate that snowline-induced dust substructures should be deeper and wider in near-infrared scattered light images compared to millimeter/submillimeter images. However, no clear substructures are observed in the Subaru HiCIAO H-band polarized intensity image of GM Aur from Oh et al. (2016). Another useful test may be to investigate whether the B84 ring is indeed warped, since snowlines are not expected to induce warps while planets may.

The new ALMA observations reveal an unresolved emission source inside the central cavity of the GM Aur disk. The flux of this feature is not straightforward to define due to the more extended low-lying emission in the cavity at 1.1 mm, but we can compare the two ALMA bands by defining a common area over which the flux is measured. Point-spread function (PSF) artifacts are a concern for faint emission features, particularly inside a bright emission cavity observed with ALMA's long baselines, but the observations at the two bands are presumably affected similarly because they were observed with the same array configurations. Based on the continuum images used to generate the high-resolution spectral index profiles in Figure 3, the flux within a diameter of 100 mas is 0.28 ± 0.02 mJy at 1.1 mm and 0.09 ± 0.02 mJy at 2.1 mm. The measurement region is chosen to be slightly larger than the synthesized beam. These fluxes correspond to a spectral index of ${1.9}_{-0.4}^{+0.5}$.

While the absolute uncertainties on the spectral index are large, the normalized continuum profiles in the two bands (Figure 3) indicate that the intensity changes less steeply as a function of frequency at the disk center compared to the emission rings. The presence of larger particles in the inner disk may explain this difference, since large dust grains are expected to drift rapidly toward the star (e.g., Weidenschilling 1977). Given that SED models of the GM Aur disk consistently find that the inner disk is highly dust-depleted (e.g., Chiang & Goldreich 1997; Calvet et al. 2005; Espaillat et al. 2010), it seems less likely that high optical depth can explain the lower spectral index of the inner emission. Nevertheless, if all of the continuum emission is due to thermal dust emission, then SED model results appear to have underestimated the dust content in the cavity. The best-fit SED model of GM Aur from Macías et al. (2018), which is based on Espaillat et al. (2010) but adjusted for the new Gaia distance, infers the presence of an optically thin inner disk extending from 0.17 to 0.85 au. Their best-fit inner disk dust mass of ∼4 × 10⁻¹² M_⊙ (assuming a_min = 0.005 μm, a_max = 0.25 μm, and p = 3.5) yields a 2.1 mm flux of ∼10⁻³ μJy, which is many orders of magnitude lower than the observed flux. McClure et al. (2013) show that it is difficult to determine from SEDs whether millimeter-sized grains are present in the inner disk due to their minor contribution to the near-infrared excess. Thus, GM Aur's SED does not necessarily rule out the presence of a millimeter (or larger) grain population that is responsible for the emission detected by ALMA.

It has been an enduring puzzle how disks such as GM Aur maintain high accretion rates inside strongly depleted dust cavities (e.g., Espaillat et al. 2014; Rosenfeld et al. 2014). Thus, it is interesting to consider how much material would have to be contained in the inner disk to sustain GM Aur's accretion rates, which have been measured to range from 3.9 × 10⁻⁹ M_⊙ yr⁻¹ to 1.96 × 10⁻⁸ M_⊙ yr⁻¹ (Ingleby et al. 2015; Robinson & Espaillat 2019). Assuming a gas-to-dust ratio of 100, maintaining a representative accretion rate of ∼10⁻⁸ M_⊙ yr⁻¹ for 10⁵ yr would require a dust mass of 10⁻⁵ M_⊙ in the inner disk if it is not replenished by the outer disk. This would represent 0.5%–4% of GM Aur's estimated total dust mass (McClure et al. 2016). Using our power-law fit to the temperature structure from Macías et al. (2018), we estimate that an inner disk with a radius of 2 au should have τ_abs = Σ_dκ_abs/cos i ∼ 0.1 to reproduce the measured inner disk flux at 1.1 mm. A dust mass of 10⁻⁵ M_⊙ spread uniformly over this disk would result in τ ≫ 1 at millimeter wavelengths if we use the DSHARP opacities and assume a_max = 1 mm and p = 3.5. However, above a_max ∼ 1 mm, opacities drop rapidly. Values of p at or shallower than 2.5 and a_max of about 1 m can yield sufficiently low emission. Invoking solids this large, though, is challenging insofar as they are expected to drift into the star within hundreds of years (e.g., Weidenschilling 1977). Alternatively, the inner disk might be able to sustain GM Aur's accretion without requiring extremely large solids if the gas-to-dust ratio is significantly enhanced.

If the dust emission at the disk center is optically thin, contributions from non-dust emission can help explain the low spectral index. Using Very Large Array (VLA) centimeter observations of GM Aur, Macías et al. (2016) estimated that free–free emission from ionized gas contributed 76 μJy to the flux measured at 3.0 cm. This emission was attributed to a combination of an ionized radio jet and a photoevaporative wind in the inner disk. Extrapolating the amount of free–free emission expected at millimeter wavelengths from existing data is not straightforward because the spectral index changes as the free–free emission becomes optically thin, which typically occurs between 1 and 10 cm (e.g., Eisner et al. 2018). Furthermore, it is unknown whether the free–free emission measured with the VLA is co-spatial with the inner disk emission measured with ALMA because the synthesized beam for the 3 cm observations is 10× larger. To make a conservative estimate of whether free–free emission could be detectable at ALMA wavelengths, we assume that it becomes optically thin at 3.0 cm. Since optically thin free–free emission scales as ${F}_{\nu ,\mathrm{ff}}\propto {\nu }^{-0.1}$ (e.g., Eisner et al. 2008), the flux due to free–free emission is expected to be ∼58 μJy at 2.1 mm and ∼55 μJy at 1.1 mm. This represents about two-thirds and one-fifth of the inner disk flux measured at 2.1 and 1.1 mm, respectively. The spectral index of the inner disk is within the range computed for optically thick winds (e.g., Reynolds 1986), but it is unlikely that the inner disk millimeter emission of GM Aur originates from optically thick free–free emission. Using the ${\alpha }_{\mathrm{free}-\mathrm{free}}=0.75$ from centimeter observations of GM Aur in Macías et al. (2016), a direct extrapolation results in a millimeter wavelength flux that is several times higher than the measured inner disk flux. Magnetic reconnection in the stellar corona is also hypothesized to generate significant emission at millimeter wavelengths (e.g., Massi et al. 2006; Salter et al. 2008), but this mechanism is usually associated with variability on timescales of a few hours, whereas the GM Aur fluxes are consistent between execution blocks. Ultimately, to disentangle the contributions of dust emission from other sources of emission in the inner disk, high angular resolution observations at other frequencies are necessary.

6.3.1. Ring Properties

6.3.2. Central Cavity Properties

6.3.3. Spectral Index Behavior