Radial Surface Density Profiles of Gas and Dust in the Debris Disk around 49 Ceti

49 Ceti was observed with ALMA on 2013 November 4. The observation incorporated 28 12 m antennas with baseline lengths between 17 and 1284 m using the Band 7 receivers under excellent weather conditions (∼0.4 mm of precipitable water vapor).

Observations cycled every 9 minutes between 49 Ceti (40 minutes total on source) and the nearby quasar J0132–1654, which was used as the primary gain calibrator. Bandpass calibrations were conducted at the beginning of the scheduling block using the quasar J0006–0623, while the absolute flux scale was calibrated using observations of the quasar J0334–401.

The correlator was configured with two broad, low-resolution continuum bands centered at 357 and 344 GHz, each with the maximum bandwidth of 2 GHz and channel width of 15.6 MHz to maximize continuum sensitivity. Another broad-band with 1.9 GHz bandwidth and 488 kHz channel spacing was centered at 355 GHz in order to observe the HCN(4–3) transition at a frequency of 354.50547 GHz; the rest of the band was incorporated into the continuum measurement to improve sensitivity. The fourth spectral window, with a bandwidth of 234 MHz and channel spacing of 61 kHz, was centered on the CO(3–2) transition at a frequency of 345.79599 GHz.

The visibilities were calibrated and imaged using the CASA and MIRIAD software packages, including the standard reduction scripts provided by ALMA. Initial phase correction was accomplished using the water vapor radiometry system, followed by a system temperature calibration to account for changing atmospheric and instrumental conditions. After standard bandpass, flux, and gain calibrations, the resulting calibrated visibilities from the three broadband spectral windows were spectrally averaged to create a single continuum channel with an average frequency of 352 GHz, corresponding to a wavelength of 850 μm. The continuum, HCN(4–3), and CO(3–2) visibilities were imaged using standard Fourier inversion (assuming natural weighting) and the CLEAN algorithm. Synthesized beam sizes (major and minor axis lengths of the clean beam), position angles, and rms noise levels in the resulting maps (imaged at native spectral resolution for the two lines) are listed in Table 1. The flux scale is subject to a standard 10% systematic calibration uncertainty.

We observed 49 Ceti in the 1.3 mm continuum with the Combined Array for Research in Millimeter-wave Astronomy (CARMA) over seven nights in late 2012 and early 2013 for a total of 15 hr. Our observations used the 15-element subarray that consists of six 10.4 m diameter antennas and nine 6.1 m diameter antennas in the compact (D) configuration, providing baseline lengths between 11 and 147 m. Table 2 contains the date of observation, length of observation per night on target, number of antennas used, approximate opacity at 230 GHz (${\tau }_{230\mathrm{GHz}}$), local oscillator (LO) frequency, and derived flux of the gain calibrator, the quasar J0132–169.

Table 2.  Observational Parameters for CARMA Data

Date	Observation Time	Number of	${\tau }_{230\mathrm{GHz}}$	LO Freq (GHz)	Derived Flux
	on 49 Ceti (hr)	Antennas			of J0132–169 (Jy)
2012 Nov 4	2.48	13	0.27–0.37	227.5215	0.858
2012 Nov 14	1.79	13	0.17–0.23	227.5185	0.619
2012 Nov 15	2.23	9	0.20–0.23	227.5184	0.659
2012 Nov 23	2.51	15	0.20–0.23	227.5158	0.637
2012 Nov 25	2.52	14	0.17–0.21	227.5153	0.583
2012 Nov 27	1.56	13	0.17–0.18	227.5147	0.600
2013 Feb 16	1.86	13	0.11–0.26	227.5150	0.375

Download table as:  ASCIITypeset image

The weather was adequate, with 230 GHz opacity varying between 0.11 and 0.37. Data from malfunctioning antennas and from periods of time when the phase rms varied widely were not used in imaging. Due to the southern declination of the source relative to the array's location, the synthesized beam was elongated in the north–south direction. The correlator was configured to take advantage of the maximum continuum bandwidth available, with the LO set to prevent the CO(2–1) line from appearing in any of the 16 spectral windows. Each window had a bandwidth of 487.5 MHz, for an aggregate bandwidth of 7.8 GHz per polarization.

Each night of observation began with short exposures of Uranus (5 minutes) to set the flux scale, followed by the bright radio source 3C84 (10 minutes) to act as a backup flux calibrator and as the bandpass calibrator. Observations cycled through the quasar J0132–169, 49 Ceti, and J0204–170. The quasar J0132–169, located 13 from 49 Ceti, was used as the gain calibrator. J0204–170, another nearby quasar, allowed us to assess the quality of the phase transfer from J0132–169. A point-source fit to the J0204–170 visibilities using the MIRIAD task uvfit returns position offsets of 0056 in α and 015 in δ, indicating that the quality of the phase transfer is good and any absolute astrometric uncertainties are far smaller than the synthesized beam (30 × 22 at a position angle of 28°). A conservative estimate for the systematic uncertainty in the flux scale, set by Uranus, is 10%.

Strong continuum emission from 49 Ceti is detected at a wavelength of 850 μm. The left panel of Figure 1 shows a naturally weighted image of the continuum emission from the disk. In order to derive an initial estimate of the integrated flux density, we use the MIRIAD task uvfit to fit an elliptical Gaussian to the visibilities on baselines $\leqslant 50$ kλ because these correspond to the largest angular scales in the data and are less affected by the small-scale structure introduced by the gap in the center of the disk. This method results in a measured total flux density of 17 ± 3 mJy. Taking into account the ∼10% systematic flux uncertainty that applies to both data sets, there is no statistically significant difference between the ALMA 850 μm flux density and the JCMT/SCUBA2 850 μm flux density of 12.1 ± 2.0 mJy (Greaves et al. 2016; see also Song et al. 2004).

Zoom In Zoom Out Reset image size

In order to estimate the geometry of the ring, we also use uvfit to fit a thin ring to the visibilities (including all baselines, unlike the flux estimate above). This fit results in a major axis diameter of 3822 ± 0003 (233.0 ± 0.2 au) and a minor axis diameter of 0626 ± 0003 (38.2 ± 0.2 au) at a position angle (defined east of north) of $-70\buildrel{\circ}\over{.} 9\pm 0\buildrel{\circ}\over{.} 4$. If we assume a circular geometry for the disk, this axis ratio implies an inclination to the line of sight of $80\buildrel{\circ}\over{.} 6\pm 0\buildrel{\circ}\over{.} 4$, and a radius of approximately 117 au for the brightest part of the disk. These values are marginally consistent with those derived from subarcsecond Keck images of thermal emission from the inner disk by Wahhaj et al. (2007), who measured a position angle of $-55^\circ \pm 10^\circ$ and inclination of $60^\circ \pm 15^\circ$. They are also consistent with the values derived from CO emission in the outer gas disk observed with the Submillimeter Array by Hughes et al. (2008a), who derive a position angle of $-70^\circ \pm 10^\circ$ and an inclination of $90^\circ \pm 5^\circ$, and with the geometry derived from Herschel PACS images by Moór et al. (2015a), yielding a position angle of −706 ± 12 and an inclination of 814 ± 39.

Figure 2 displays the CARMA image of 49 Ceti at a wavelength of 1300 μm. The combination of seven nights (∼15 hr) of data resulted in a 4σ detection (in the image domain) of flux from the disk at the expected position of the star. The phase center of the observations on each night (the origin of the axes in Figure 2) corresponds to the expected position of the star after accounting for the proper motion of 94.8 ± 0.3 mas yr⁻¹ in RA and −3.14 ± 0.19 mas yr⁻¹ in declination (van Leeuwen 2007). The data provide a spatial resolution of ∼22 × 30 when imaged with a Briggs robust parameter of 0.5. The disk is apparently unresolved along both the major and minor axis, according to the results of an elliptical Gaussian fit performed with the MIRIAD task uvmodel. We derive a total flux of 2.1 ± 0.7 mJy by fitting a point source to the visibilities with the MIRIAD command uvfit. In comparison with the size of the continuum ring observed in the ALMA data, it is somewhat surprising that the CARMA data appear unresolved. Given the very low signal-to-noise ratio of the data, it is impossible at this stage to determine whether the morphology of the observed CARMA data reflects a true difference in morphology between grains of different sizes or whether it is simply the result of sensitivity limitations and random noise in low surface brightness data. The surface brightness of the best-fit model that reproduces the ALMA visibilities described in Section 4.1 below predicts a surface brightness at CARMA wavelengths that is essentially undetectable outside the central 1.5 beams at the sensitivity of the data. It is therefore possible that some extended flux is missed by the CARMA continuum data, and the total flux measurement may be an underestimate.

Zoom In Zoom Out Reset image size

We detect strong CO(3–2) emission from the 49 Ceti disk across approximately 250 channels, each with a width of 0.05 km s⁻¹. We illustrate the basic properties of the CO(3–2) emission in Figure 3. The emission is symmetric and regular and appears consistent with Keplerian rotation about the central star.

Zoom In Zoom Out Reset image size

We measured the total flux over the region enclosed by the 3σ contours in the zeroth moment map using the MIRIAD task cgcurs. We derive a total flux density of 6.0 ± 0.1 Jy km s⁻¹. This flux is consistent with the integrated flux from spatially unresolved observations of the CO(3–2) spectrum with the James Clerk Maxwell Telescope (JCMT) by Dent et al. (2005). They report an integrated CO intensity of 0.34 ± 0.07 K km s⁻¹, which corresponds to an integrated flux of 6.6 ± 1.4 Jy km s⁻¹ (Meeus et al. 2012). Both measurements are subject to the typical ∼10% systematic flux uncertainty at this frequency. The ALMA observations therefore recover most, if not all, of the CO flux from the disk. The maximum width of the 3σ contour along the major axis of the disk is 73 ± 01, corresponding to an outer diameter of the CO disk of 430 au at a distance of 59 pc. This indicates that the expected outer radius of CO emission from the disk is approximately 220 au, though the convolution of the beam with the true sky brightness distribution may result in a slight overestimate of this value. Since most of the flux is recovered, it is unlikely that a significant amount of CO flux extends farther from the disk at a level below the sensitivity of the ALMA map. Modeling of the visibilities is required to constrain the outer radius more precisely (see Section 4.2 below).

The peak flux density of the two limbs in the zeroth moment map is identical to within the uncertainties: 0.47 ± 0.02 Jy beam⁻¹ for the northwestern peak and 0.46 ± 0.02 Jy beam⁻¹ for the southeastern peak. The separation between these peaks is 143 ±003, corresponding to a spatial separation of 84 au (radius 42 au) at a distance of 59 pc. If the CO is optically thin, the lack of a central emission peak in the CO(3–2) moment map is an indication of a central cavity in the CO emission; however, for optically thick edge-on disks it is possible to artificially create a double-peaked structure via optical depth effects, as illustrated in Wolf et al. (2008).

We perform a simple two-dimensional elliptical Gaussian fit to the visibilities using the MIRIAD task uvfit to constrain the position angle of the CO disk. This method returns a position angle of −720 ± 02, which is consistent to within the uncertainties with the position angle of the continuum data and previous observations of the CO disk. It is also consistent with the position angle of −700 ± 05 derived from the location of the flux peaks in the CO map along the northwest and southeast limbs of the disk. The PV diagram in Figure 3 represents a single-pixel cut across the disk midplane, at a position angle of −72° and passing through the proper motion-corrected star position at the phase center of the observations.

We do not detect any HCN(4–3) emission from the 49 Ceti disk to within the uncertainties in the data. At the full 0.41 km s⁻¹ spectral resolution of the observations, the rms noise in each channel is 3.0 mJy beam⁻¹, with a beam size of 056 × 043. Since the FWHM of the CO line is 8 km s⁻¹, we binned the data by a factor of five in the spectral domain to achieve a higher signal-to-noise per channel, and still detected no HCN emission at the $\gt 3\sigma$ level with an rms noise of 1.7 mJy beam⁻¹ in each 2.5 km s⁻¹-wide channel.

We simultaneously model the spatially resolved 850 μm ALMA visibilities and the unresolved SED in order to constrain geometric properties of the disk and determine the basic characteristics of the constituent dust grains. The modeling process and simplifying assumptions are described in detail in Ricarte et al. (2013) and Steele et al. (2016); we assume a geometrically flat and optically thin dust disk, astrosilicate opacities to determine the temperature of the dust grains as a function of their size, and a simplified emission efficiency of ${Q}_{\lambda }=1-{e}^{-{(\tfrac{\lambda }{2\pi a})}^{-\beta }}$ (following Williams et al. 2004). The hybrid approach of assuming tabulated astrosilicate opacities for the temperature calculation, while also using a parameterized emission efficiency, allows us to approximate the effects of a grain size distribution on the SED while maintaining sufficient computational efficiency to take advantage of the powerful Markov Chain Monte Carlo (MCMC) approach to fitting the data and understanding the uncertainties on the model parameters. A similar approach is used by Pawellek et al. (2014), who also compare their version of the modified blackbody with a true grain size distribution fit and find that the results are generally consistent.

The dust grain properties are characterized by two parameters: the characteristic grain size a and long-wavelength power-law index of grain emission efficiency β. We calculate temperatures on the basis of the characteristic grain size, assuming a compact astrosilicate composition (Draine 2003). We determine the opacity, ${\kappa }_{\mathrm{tot}}(a,\lambda )$, and albedo, $\omega (a,\lambda )$, of the grains using the Mie theory implementation in MCFOST (Pinte et al. 2006, 2009), and use these properties to calculate the temperature of the dust assuming radiative equilibrium.

We compare each model to the SED of the 49 Ceti system at observed photometric points gathered from the literature (displayed in Table 3). We model the SED with two components: (1) a Kurucz–Lejeune model photosphere with a surface gravity of $\mathrm{log}g=4.5$, an effective temperature of T_eff = 10,000 K, and a solar metallicity of Z = 0.01 (Chen et al. 2006), and (2) an extended, spatially resolved debris disk. We calculate a ${\chi }^{2}$ metric as a test of goodness of fit for the SED using only data points at wavelengths greater than 5 μm, since shorter wavelengths are strongly dominated by the stellar photosphere and are not affected by changes to the disk model parameters. We note that more recent sources derive a lower effective temperature of 8900–9500 K and slight reddening toward 49 Ceti (Montesinos et al. 2009; Roberge et al. 2013; Moór et al. 2015a). Since we do not include wavelengths <5 μm in our fit, the details of the model in the optical/UV range do not affect the excess that we derive at longer wavelengths. The ALMA flux is not included in the calculation of the ${\chi }_{\mathrm{SED}}^{2}$ because this flux is implicitly included in the creation of the model image at 850 μm that is compared with the visibilities when calculating the goodness of fit parameter for the visibilities, ${\chi }_{\mathrm{vis}}^{2}$. The Spitzer IRS spectrum (taken from the CASSIS database; Lebouteiller et al. 2011) originally consisted of 360 points between 5 and 35 μm, but was averaged down to 9 for the sake of computational efficiency. Each of the nine points was included as a single data point (with appropriately weighted errors) when calculating the SED ${\chi }^{2}$ value. The errors reported in the table for these points are the absolute calibration uncertainty (assumed to be 10% of the flux measurement) and the statistical uncertainty added in quadrature.

In addition to comparing the model to the observed SED, we also create a high-resolution image of the disk at a wavelength of 850 μm for comparison with the spatially resolved ALMA data. The model image is generated with a spatial resolution of approximately 10% of the spatial scale sampled by the longest baseline in the ALMA data, which corresponds to 3 au per pixel. The model images are sampled at the same spatial frequencies as the ALMA data using the MIRIAD command uvmodel, allowing us to compare our model disk to the data in the visibility domain where the uncertainties are well characterized. We then calculate a ${\chi }_{\mathrm{vis}}^{2}$ statistic as a goodness of fit test of the visibilities.

4.1.1. Surface Density Profiles

Our primary goal in fitting the ALMA data is to characterize the surface density structure of the 49 Ceti debris disk, based on the spatially resolved ALMA visibility data. The most common parameterization of surface density in debris disks is a simple power-law function of radius, with ${\rm{\Sigma }}(r)\propto {r}^{-p}$ and abrupt cutoffs at the inner and outer radii. For this model, we need to specify four variables to generate the surface density profile: the total disk mass, M_disk; the surface density power-law index p; the inner radius, R_in; and the outer radius, R_out. If we normalize to a radius of 100 au, the surface density profile is then

$$\begin{eqnarray}&&{\rm{\Sigma }}(r)={{\rm{\Sigma }}}_{100\mathrm{au}}{\left(\displaystyle \frac{r}{100\mathrm{au}}\right)}^{-p}\end{eqnarray} \tag{ 1 }$$

where

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{100\mathrm{au}}=\displaystyle \frac{{M}_{\mathrm{disk}}(2-p)}{2\pi {(100\mathrm{au})}^{p}({R}_{\mathrm{out}}^{2-p}-{R}_{\mathrm{in}}^{2-p})}.\end{eqnarray} \tag{ 2 }$$

However, the appearance of the naturally weighted image of the 49 Ceti disk in Figure 1, which displays a relatively gradual increase and decrease of surface brightness along the semimajor axis, indicates that perhaps a more complex surface density structure is justified. Motivated by the three models used in the fit to ALMA observations of the debris disk around HD 107146 by Ricci et al. (2015), we also attempt to fit the data with a radially unresolved surface density enhancement (i.e., a thin ring of increased surface density located at a radius R_belt), and a double power-law model in which the surface density profile changes its power-law index p at a transition radius R_T that is constrained to lie between R_in and R_out. The unresolved surface density enhancement requires the addition of two parameters to the default single power-law model: the radius of the extra belt, R_belt, and the mass of the belt, M_belt. We model the width of the belt to be 3 au across the radial dimension, which was chosen to be a small fraction of the spatial scale probed by the longest baseline in the ALMA data set. For the double power-law model, we specify three new parameters: the transition radius, R_T, the power law for the inner surface density profile, p₁, and the power law for the outer surface density profile, p₂. The surface density profile is then

$$\begin{eqnarray}{\rm{\Sigma }}(r)=\left\{\begin{array}{ll}{{\rm{\Sigma }}}_{{\rm{T}}}{\left(\tfrac{r}{{R}_{{\rm{T}}}}\right)}^{-{p}_{1}} & \,\mathrm{for}\,{R}_{\mathrm{in}}\leqslant r\leqslant {R}_{{\rm{T}}}\\ {{\rm{\Sigma }}}_{{\rm{T}}}{\left(\tfrac{r}{{R}_{{\rm{T}}}}\right)}^{-{p}_{2}} & \,\mathrm{for}\,{R}_{{\rm{T}}}\leqslant r\leqslant {R}_{\mathrm{out}}\end{array}\right.\end{eqnarray} \tag{ 3 }$$

where

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{{\rm{T}}}=\displaystyle \frac{{M}_{\mathrm{disk}}{j}_{1}{j}_{2}}{2\pi ({m}_{1}{j}_{2}{R}_{{\rm{T}}}^{{p}_{1}}+{m}_{2}{j}_{1}{R}_{{\rm{T}}}^{{p}_{2}})}\end{eqnarray} \tag{ 4 }$$

and ${j}_{1}=2-{p}_{1}$, ${j}_{2}=2-{p}_{2}$, ${m}_{1}={R}_{{\rm{T}}}^{{j}_{1}}-{R}_{\mathrm{in}}^{{j}_{1}}$, and ${m}_{2}\,={R}_{\mathrm{out}}^{{j}_{2}}-{R}_{{\rm{T}}}^{{j}_{2}}$.

The three surface density profiles are illustrated schematically in the right column of Figure 5, where dark green represents high surface density. To ensure that ${R}_{\mathrm{in}}\,\leqslant {R}_{{\rm{T}}}\leqslant {R}_{\mathrm{out}}$, within the code we specify the transition radius R_T as ${R}_{\mathrm{in}}+{\rm{\Delta }}{R}_{{\rm{T}}}$ and the outer radius R_out as ${R}_{\mathrm{in}}\,+{\rm{\Delta }}{R}_{{\rm{T}}}+{\rm{\Delta }}{R}_{\mathrm{out}}$.

4.1.2. MCMC Fitting

In order to find the best fit to the data and explore the uncertainties (including degeneracies) associated with each parameter, we utilize the affine-invariant MCMC fitting technique as described by Goodman & Weare (2010) and implemented in Python as emcee by Foreman-Mackey et al. (2013). We fit the sum of the SED and visibility ${\chi }^{2}$ values, assuming uniform priors on all parameters (though the assumption of uniform priors on those parameters sampled logarithmically is equivalent to the assumption of a log-normal prior). Using an MCMC method allows us to probabilistically sample the full parameter space described by our models, obtain a best-fit result, and place uncertainties on the model parameters by characterizing the full posterior distribution function for each parameter. Several runs were attempted with different starting positions to ensure that the solution was a global rather than a local minimum. Most of our MCMC chains included approximately 100 walkers that were allowed to move for 1000 steps, and the "burn-in" phase during which the ${\chi }^{2}$ values rapidly decreased as the walkers approached the best-fit region of parameter space typically comprised of order ∼100 steps. The burn-in phase was discarded from each MCMC run.

We experimented with fitting the spatially resolved visibilities and SED both simultaneously and independently. These runs made it clear that it was not possible to fit both the SED and visibilities with only a single characteristic dust grain size, a, and long-wavelength emission efficiency power law, β. For example, using the Akaike Information Criterion (AIC) test, the best-fit simultaneous SED and visibility model using a single power-law surface density distribution was a poorer fit to the visibilities than the best visibility-only fit at the $\gt 10\sigma$ level. This is unsurprising, since comparisons of the mid-IR imaging (Wahhaj et al. 2007) with the ALMA image in Figure 1 reveal two spatially disparate grain distributions: the small grains that dominate the thermal IR emission are concentrated within 30–60 au from the star, while the large grains that dominate the millimeter-wavelength emission exhibit peak surface brightness at a radius of more like 100 au from the star and extend out to several hundred astronomical units. These observations are consistent with modeling results described in Wahhaj et al. (2007), which indicate the need for two grain populations to reproduce both the SED and the resolved imaging in the mid-IR. We therefore introduce an "inner disk" component into our model, with a characteristic grain size of $a=0.1\,\mu {\rm{m}}$ to maintain consistency with Wahhaj et al. (2007), variable inner radius ${R}_{\mathrm{Inner}\mathrm{disk}}$, and an outer radius that coincides with the inner radius of the outer disk, R_In. Since the inner disk surface density profile is not spatially resolved enough to constrain the power-law index as a function of radius in the Wahhaj et al. (2007) data, we assume ${{\rm{\Sigma }}}_{\mathrm{Inner}\mathrm{disk}}\propto {r}^{-1}$ and vary the total mass of the inner disk ${M}_{\mathrm{Inner}\mathrm{disk}}$ and the total mass of the outer disk ${M}_{\mathrm{Outer}\mathrm{disk}}$ separately. For the model with the unresolved surface density enhancement, we also separately vary the mass of the ring, M_Ring. We find that R_in is essentially unconstrained by the spatially resolved ALMA data in the case of the double power law (since the spatial resolution is not high enough to distinguish between a sharp and gradual inner disk edge), so we therefore fix the location of R_in to be 20 au greater than ${R}_{\mathrm{Inner}\mathrm{disk}}$ for the double power-law model.

Table 4 reports both the median and $\pm 1\sigma$ intervals from the posterior distribution, as well as the global best-fit set of model parameters from each MCMC run, for the three classes of surface density structure that were used for the fits. The best-fit values are generally within $\pm 1\sigma$ of the median. For all models, disk mass parameters (${M}_{\mathrm{Inner}\mathrm{disk}}$, ${M}_{\mathrm{Outer}\mathrm{disk}}$, and M_Ring) and the grain size are sampled logarithmically, while the rest of the parameters are sampled linearly. It should be noted that since the parameter β is constrained primarily by the slope of the SED between the ALMA and CARMA flux points, the low sensitivity of the CARMA data to extended low-surface-brightness flux might bias our measurement of β to be steeper than the true value.

Table 4.  Best-fit Parameters for Three Models of the Dust Disk

	Single Power Law		Power Law w/Ring		Double Power Law
Parameter	Median $\pm 1\sigma$	Best Fit	Median $\pm 1\sigma$	Best Fit	Median $\pm 1\sigma$	Best Fit
${R}_{\mathrm{Inner}\mathrm{Disk}}$ [au]	${5}_{-1}^{+2}$	6	${4.6}_{-1.2}^{+1.5}$	3.5	${7.0}_{-1.5}^{+2.0}$	6.7
RIn [au]	${73}_{-3}^{+3}$	73	${62}_{-5}^{+5}$	60	⋯	26.7
${\rm{\Delta }}{R}_{{\rm{T}}}$ [au]	⋯	⋯	⋯	⋯	${68}_{-11}^{+7}$	68.6
${\rm{\Delta }}{R}_{\mathrm{Out}}$ [au]	${310}_{-9}^{+9}$	310	${302}_{-8}^{+9}$	303	${215}_{-9}^{+12}$	219
RRing [au]	⋯	⋯	${114}_{-3}^{+3}$	114	⋯	⋯
$\mathrm{log}({M}_{\mathrm{Inner}\mathrm{disk}}$ [M⊕])	$-{3.1}_{-0.2}^{+0.2}$	−3.0	$-{3.3}_{-0.2}^{+0.2}$	−3.4	$-{3.33}_{-0.19}^{+0.14}$	−3.40
$\mathrm{log}({M}_{\mathrm{Outer}\mathrm{disk}}$ [M⊕])	$-{1.17}_{-0.10}^{+0.08}$	−1.16	$-{1.21}_{-0.10}^{+0.09}$	−1.21	$-{1.07}_{-0.05}^{+0.04}$	−1.09
$\mathrm{log}({M}_{\mathrm{Ring}}$ [M⊕])	⋯	⋯	$-{2.26}_{-0.18}^{+0.15}$	−2.26	⋯	⋯
$\mathrm{log}(a[\mu {\rm{m}}$])	${0.08}_{-0.14}^{+0.11}$	0.09	${0.09}_{-0.13}^{+0.10}$	0.06	${0.23}_{-0.05}^{+0.05}$	0.22
β	${1.17}_{-0.05}^{+0.05}$	1.18	${1.17}_{-0.05}^{+0.05}$	1.17	${1.23}_{-0.03}^{+0.03}$	1.22
p	${1.29}_{-0.11}^{+0.11}$	1.27	${0.80}_{-0.17}^{+0.19}$	0.75	⋯	⋯
p1	⋯	⋯	⋯	⋯	$-{2.5}_{-2.2}^{+0.8}$	−2.7
p2	⋯	⋯	⋯	⋯	${1.44}_{-0.20}^{+0.15}$	1.50
i [°]	${79.4}_{-0.4}^{+0.4}$	79.3	${79.3}_{-0.4}^{+0.4}$	79.3	${79.0}_{-0.2}^{+0.2}$	79.2
PA [°]	$-{71.4}_{-0.5}^{+0.4}$	−71.3	$-{71.4}_{-0.4}^{+0.4}$	−71.2	$-{71.6}_{-0.5}^{+0.4}$	−71.7

Download table as:  ASCIITypeset image

Figure 5 displays a graphical summary of the fitting results, including the best-fit model and residuals for each of the three classes of models. A schematic of each model is presented in the right column. The left column compares the observed SED (black points; see Table 3) with the various components of the model SED and the Spitzer IRS spectrum. All three of the models are consistent with the data to within the uncertainties, as illustrated by the absence of any statistically significant ($\gt 3\sigma$) features in the residual image. This lack of residual flux also indicates that the surface density structure of the dust disk is consistent with axisymmetry, despite the slight apparent mismatch of peak flux levels in the ansae of the image derived from the ALMA data (Figure 1). Figure 4 emphasizes this point by plotting the tapered residuals, which also show no significant features, as would be expected if there were a more subtle but significant asymmetry on the scale of a few beams rather than the 04 scale of the naturally weighted image. Plots of the best-fit surface density profile for the three classes of models are presented in Figure 6. The three models are consistent in showing an overall trend of decreasing surface density with radius, with the inner radius of the outer disk located between ∼50 and 70 au and the outer radius close to 300 au (the inner and outer radius are most uncertain with the double power-law model, due to degeneracy between the slope of the power law and the location of the edges). The double power law and the single power law with a ring both indicate a maximum in the surface density profile that falls at a radius of ∼90–110 au.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

While the two-dimensional map of the model and residuals indicates no significant asymmetry, subtle differences in the distribution of flux are evident in the three model images, and it is possible that the ensemble of visibilities may be able to distinguish statistically between the three modeled surface density distributions (akin to boosting the signal-to-noise ratio by considering an azimuthally averaged version of the surface density distribution rather than only point-to-point variations in flux). To explore this possibility, we employ both the AIC and the related Bayesian Information Criterion (BIC), which allow us to make statistical comparisons between the likelihood functions of models with different numbers of parameters. A key assumption of the AIC is that the posterior distribution functions are Gaussian or nearly Gaussian (Liddle 2007); the BIC is a similar test that typically imposes a larger penalty for additional free parameters. All of the parameters from the MCMC chain are well constrained and fill out posterior distributions that fit the criteria. According to the AIC, there is marginally significant improvement in the fit to the data when a local surface density enhancement is included in the model: the double power-law model is a better fit than the single power-law model at the 2.2σ level, while the power law with a thin ring is a better fit than the single power-law model at the 2.9σ level. By contrast, the AIC detects no significant difference between the double power law and the single power law with a ring (0.5σ). The BIC returns similar results, except that ΔBIC is $\gt 10$ when comparing the single power-law fit with both the double power law and the power law with a ring, corresponding to "very strong" evidence that the models that include an enhancement in density at a radius of ∼100 au. By contrast, ΔBIC is 4 for the comparison between the double power law and the single power law with a ring, which falls into the category of "positive" but not "strong" evidence that the power law with a ring is preferred. There is also a philosophical point to be made that the tendency of both the double power-law model and the single power law with a ring to converge on a solution with a density increase at ∼100 au is suggestive of a density enhancement, since either model could instead have converged to the single power-law solution. On the whole, we take these statistical tests to indicate that models of the surface density that include a surface density enhancement at a distance of ∼100 au from the central star better represent the physical nature of 49 Ceti's dust disk than a single power law. The current data cannot distinguish between a radially broad enhancement like the double power-law model and a radially narrow enhancement like the power law plus ring.

In order to understand the underlying temperature and density structure of the gas disk, we fit the new ALMA CO(3–2) visibilities in tandem with previously published SMA CO(2–1) observations of the gas disk (Hughes et al. 2008a). The combination of two lines provides constraints on the excitation temperature of the gas in the disk, while the superior sensitivity and angular resolution of the ALMA data provides the best constraints on the brightness profile and kinematics of the gas disk.

Because previous observations of gas-bearing debris disks have derived low excitation temperatures (e.g., Kóspál et al. 2013; Flaherty et al. 2016), there is a high probability that the disk is not in local thermodynamic equilibrium (LTE; see also Matrà et al. 2015). The integrated line ratio of 3.0 in the 49 Ceti system, corresponding to a disk-averaged excitation temperature of 32 K, is also consistent with a non-LTE scenario. We therefore model the gas disk using the non-LTE Line Modeling Engine (LIME) radiative transfer code written by Brinch & Hogerheijde (2010). Since each set of model parameters takes on the order of an hour to converge, it is not feasible to fit the data using the MCMC method that we employed for the continuum emission. Instead, we carry out a simple ${\chi }^{2}$ minimization via a grid search technique, which is good at locating the best-fit combination of model parameters in a well-behaved space, but does not efficiently characterize the uncertainties on each parameter. To avoid local minima, we start by searching a broad swath of parameter space using a coarsely sampled grid, and gradually narrow down our search to a smaller, more finely sampled grid. The ${\chi }^{2}$ value that is minimized is the sum of the raw ${\chi }^{2}$ values for the SMA and ALMA data.

Since the underlying H₂ abundance is unknown, we run two sets of models: the first assumes a constant, nominal CO/H₂ abundance of 10⁻⁴, consistent with ISM abundances. The second assumes a constant CO/H₂ abundance of 1, intended to be representative of a secondary gas origin. We recognize that the chemistry of the 49 Ceti disk could in fact deviate substantially from these values, but a comparison between the model outputs should provide some guidance as to which results are independent of assumptions about the origin of the gas in the system.

We assume a simple power-law index for the temperature structure of the CO disk, and a similarity solution for the surface density (see, e.g., Hughes et al. 2008b). In cases where the CO emission is optically thin, the power-law indices describing the radial variations of temperature and density may be degenerate. We therefore assume that the temperature decreases with distance from the star R as $1/\sqrt{R}$, which is a standard assumption based on blackbody equilibrium calculations that may or may not apply to the unusual physical conditions in the disk around 49 Ceti (we relaxed this assumption for the condition of low H₂ abundance, in which the emission is more optically thick, although the best-fit value for the temperature power law remained at 0.5). We also assume that the temperature structure of the gas is vertically isothermal, and depends only on distance from the star: $T(R)={T}_{100}{\left(\tfrac{R}{100\mathrm{au}}\right)}^{-0.5}$. For the case of CO/H₂ = 10⁻⁴, we use the default assumption of a standard ISM composition for the gas that is 80% molecular hydrogen by mass with a mean molecular weight ${m}_{0}=2.37{m}_{{\rm{H}}}$, where m_H is the mass of the hydrogen atom; for the case of CO/H₂ = 1, we assume a mean molecular weight of 14, consistent with the dominant mass component being the C and O atoms released from CO photodissociation (Matrà et al. 2017). We then calculate the sound speed as a function of radius: ${c}_{s}(R)=\sqrt{\tfrac{{k}_{B}T(R)}{{m}_{0}}}$. The scale height is then $h(R)=\tfrac{\sqrt{2}{c}_{s}(R)}{{\rm{\Omega }}(R)}$, where ${\rm{\Omega }}(R)=\sqrt{\tfrac{{{GM}}_{* }}{{R}^{3}}}$ and ${M}_{* }$ is the stellar mass. We parameterize the surface density Σ as a function of radius as

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

where R_c is the characteristic radius and γ describes the power-law drop-off of surface density with radius when $R\ll {R}_{c}$. We also include an inner cutoff in the disk surface density at radius R_in. The volume mass density as a function of radius and height above the midplane z is then

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

We do not include a prescription for CO freeze-out since the freeze-out timescale is inversely proportional to the surface area of grains per unit volume (Flower et al. 2005), which is likely to be extremely small in the low-density environment of the debris disk around 49 Ceti. We further assume a constant turbulent line width of 10 m s⁻¹, which is substantially smaller than the spectral resolution of the data, and cylindrical Keplerian rotation, where $v(R)=\sqrt{\tfrac{{{GM}}_{* }}{R}}$.

Initially, we varied the following parameters: $\{{M}_{* },{M}_{\mathrm{disk}},{R}_{c},{R}_{\mathrm{in},\mathrm{gas}}$, γ, T₁₀₀, i}. To minimize the number of free parameters, we fixed the position angle PA at the best-fit value of −714 from the fits to the continuum visibilities, which is intermediate between the two best-fit CO values of −720 (from the two-dimensional Gaussian fit to the integrated CO visibilities) and −700 (from the position of the two flux peaks in the CO map). After the initial coarse ${\chi }^{2}$ minimization, we performed a one-dimensional ${\chi }^{2}$ minimization on the systemic velocity, considering values between 2 and 3.5 km s⁻¹. We settle on a best-fit LSR velocity of 2.78 ± 0.13 km s⁻¹.

Strong residuals visible along the extreme edges of the disk minor axis suggested that the vertical extent of the CO emission was perhaps incorrect, for both classes of models. We therefore introduced an additional multiplicative factor h₀ that modifies the scale height so that $h(R)={h}_{0}\tfrac{\sqrt{2}{c}_{s}(R)}{{\rm{\Omega }}R}$. The scale height factor h₀ is a simple way of accounting parametrically for a disk structure that may be flatter or puffier than expected based on the standard Gaussian scale height approximation to hydrostatic equilibrium. A physical motivation for this additional factor is that, especially for a disk that is not in LTE, the kinetic gas temperature may differ significantly from the excitation temperature, causing us to calculate an inappropriate scale height since the temperature of the model is constrained by the line ratio rather than the kinematics. Furthermore, the Gaussian approximation to the vertical structure of the disk is not rigorous, and to be truly internally self-consistent, we would have to iteratively calculate the vertical structure by solving the equations of hydrostatic equilibrium for each new disk structure that is proposed, which would increase the computational resources required even further. Introducing this parameter decreased the ${\chi }^{2}$ value of the best fit by 2117 to 2,598,848 (reduced ${\chi }^{2}$ 1.35), which represents a significantly better fit at the $\gt 10\sigma$ level according to the AIC test.

In the final fit, we therefore varied the following parameters: {${M}_{* }$, M_disk, R_c, ${R}_{\mathrm{in},\mathrm{gas}}$, γ, T₁₀₀, i, h₀}. The best-fit values, minimum and maximum grid values considered, and the finest step size in the final grid are presented in Table 5. The reported best-fit disk masses are the sum of the total H2 and CO mass. This best fit results in a reduced ${\chi }^{2}$ value of 1.166 for the SMA CO(2–1) data and 1.094 for the ALMA CO(3–2) data (the differences between the CO/H₂ = 10⁻⁴ and CO/H₂ = 1 cases are insignificant at this level of precision). The derived stellar mass of 2.1 M_⊙ is in good agreement with that derived by Roberge et al. (2013) using isochrone fitting. Figure 7 presents the residual moment zero map for the case of CO/H₂ = 10⁻⁴ (contours) overlaid on the moment zero map of the original data (colors), with the fixed position angle of 714 overplotted as a line (the residuals are very similar for the case of CO/H₂ = 1. The residuals in the integrated intensity map peak at roughly 16% of the peak value of the observed moment zero map, but interestingly they display a departure from axisymmetry. Along the southeast limb of the disk, the peak in the residuals is offset to the north of the disk major axis, whereas along the northwest limb of the disk the residual peak is offset to the south—a geometry suggestive of spiral arms, a warp, a second disk with a different position angle (like β Pictoris Matrà et al. 2017), or other non-axisymmetric structure in the disk. Running the fit with a variable position angle yields a best-fit PA of −779, which differs significantly from the continuum position angle (−714 ± 05) as well as the two best-fit position angle values obtained from fitting the integrated CO map with a 2D Gaussian (−720 +/− 02) and from fitting the position angle of the two flux peaks in the integrated CO map (−700 +/− 05). While the position angle of the integrated CO map is statistically indistinguishable from that of the continuum map, the results of the channel map fit support a small but statistically significant deviation from axisymmetry in the CO map that manifests as a difference in position angle.

Zoom In Zoom Out Reset image size

Table 5.  Best-fit Parameters for Parametric Model Fit to ALMA CO(3–2) and SMA CO(2–1) Data

Parameter	Best Fit	Range of Values	Best Fit	Range of Values	Finest Spacing
	CO/H2 = 10−4		CO/H2 = 1
${M}_{* }$ [M⊙]	2.1	1.5 to 3.0	2.1	1.5 to 3.0	0.1
Log(Mdisk [M⊙])	−6.5	−8 to −5	−9.0	−10 to −7	0.2
Rc [au]	140	50 to 200	140	50 to 200	20
${R}_{\mathrm{in},\mathrm{gas}}$ [au]	20	10 to 100	20	10 to 100	10
γ	−0.5	−2.0 to 1.6	−1.8	−2.0 to 1.6	0.2
T100 [K]	40	10 to 100	14	10 to 100	2
i [°]	79.8	60 to 90	79.1	60 to 90	0.1
h0	0.75	0.2 to 1.2	2.1	0.2 to 3.0	0.1

Download table as:  ASCIITypeset image

An important point to consider is the optical depth of the CO emission from 49 Ceti. The best-fit LIME models are at least marginally optically thin across a significant area of the disk for both classes of models: for the case of CO/H₂ = 10⁻⁴ fewer than 0.4% of the pixels with flux greater than 1% of the peak value in the CO(3–2) model have $\tau \gt 1$, whereas for the case of CO/H₂ = 1 the optical depth is higher with ∼20% of the pixels with flux $\gt 1$% of the peak exhibiting $\tau \gt 1$ (in this latter case the optical depth rises above 10 across only 10% of the disk area). It is also possible to estimate the optical depth of the emission from the 3D data cube. If we assume LTE and calculate the optical depth under that assumption from the location exhibiting the peak flux per beam per spectral resolution element in the CO(3–2) position–position–velocity channel maps (since this location presumably has the highest optical depth), following, e.g., Matrà et al. (2017):

$$\begin{eqnarray}&&{\tau }_{\nu }=\displaystyle \frac{h\nu }{4\pi {\rm{\Delta }}\nu }\left(\displaystyle \frac{{x}_{l}}{{x}_{u}}{B}_{{lu}}-{B}_{{ul}}\right)\displaystyle \frac{4\pi {d}^{2}{F}_{{ul}}}{h\nu {A}_{{ul}}{\rm{\Delta }}A}\end{eqnarray} \tag{ 7 }$$

where h is Planck's constant, ν is the frequency of the spectral line, ${\rm{\Delta }}\nu$ is the line width (1.7 km s⁻¹ at the location of peak flux), x_l and x_u are the fractional populations in the lower and upper level of the transition, respectively, B_lu, B_ul, and A_ul are the Einstein coefficients for the transition, d is the distance to the star, F_ul is the integrated flux density of the line (0.28 Jy km s⁻¹ at the location of peak flux in the 3D data cube, integrating across velocity space), and ${\rm{\Delta }}A$ is the linear area of sky intersected by the emission along the line of sight (in this case, one synthesized beam ellipse with major and minor axis lengths of 047 and 032, respectively). Assuming that the level populations are thermalized at a kinetic temperature of 40 K (appropriate for the ∼100 au distance of the emission peak from the star based on our two-line fit to the excitation temperature), this calculation yields a peak optical depth of 0.5. While the calculation is sensitive to the assumed temperature, even in the case where the excitation temperature corresponds to the 32 K value indicated by the integrated line ratio the peak optical depth is only slightly higher, with a value of 0.66. Departures from LTE can raise the value further, but the estimate is generally consistent with the output of the LIME models for the system.

The radial extent of the CO emission is also noticeably smaller than that of the dust continuum emission. While it is difficult to compare directly, the power-law surface density profile that we assumed for the dust with the similarity solution profile that we assumed for the gas, it is nevertheless clear that the observed outer radius of CO emission (220 au when convolved with the beam, as we reported in Section 3.3) is substantially smaller than the ∼310 au outer radius of the best-fit models for the continuum emission. In the context of our assumption of a CO/H₂ ratio of 10⁻⁴, the best-fit midplane density structure of the CO disk falls below the critical density of the CO(3–2) transition at a radius of ∼105 au. A substantially lower H₂ density would result in the density structure reaching this value at an even smaller radius. This calculation is also affected by the departure from LTE and the uncertain kinetic versus excitation temperatures in the disk.

Comparison between the model with primordial abundances (CO/H₂ = 10⁻⁴) and the model with low H₂ content (CO/H₂ = 1) results in parameters that obtain consistent results for stellar mass, inner radius, characteristic radius, inclination, and a scale height that is spatially unresolved (though the scale height factor is larger as expected for the higher mean molecular weight gas). The disk mass of the low-H₂ model is lower, which is easily understood due to the lack of (invisible) H₂ that populated the primordial model and made up the bulk of the mass in that condition. The surface density power law is steeper for the low-H₂ model; the reason is not entirely clear, but it may be related to the higher optical depth of the low-H₂ model, which requires a more extreme density slope to generate a sufficient flux contrast between the inner and outer disks (there is also likely to be a higher uncertainty on γ in this more optically thick case). The temperatures listed in Table 5 are the kinetic temperatures that LIME takes as its input; the unreasonably low value of 14 K for the kinetic temperature of the low-H₂ model is probably not physical but may be related to the need to fine-tune the collision partner densities, which is beyond the scope of this paper's attempts to present a broad initial characterization of the gas and dust properties of the disk. Nevertheless, the comparison between the two classes of models is instructive in revealing that the increasing surface density as a function of radius, along with the basic geometrical and kinematic properties of the system, is a result that is robust to large changes in the collision partner density in the disk.

The resolution and sensitivity now available with ALMA are making it possible to characterize the radial surface density profile of the large grains in debris disks for the first time. Unlike the smaller, μm size grains probed by optical and near-infrared observations of debris disks, the ALMA images are most likely primarily sensitive to large grains that trace the gravitational potential of the debris disk system and are not strongly affected by stellar radiation pressure. They are therefore thought to trace the location of the parent planetesimal belts in the disk, and to provide an indication of the collision rate (which controls the dust production rate) between the parent planetesimals. Thus far, only a handful of debris disks have had their radial surface density profiles (i.e., whether the surface density increases or decreases with radius) characterized in detail using high-sensitivity ALMA data. Some disks, like that around β Pic (Dent et al. 2014), have been imaged in sufficient detail but are not sufficiently axisymmetric to permit the derivation of a meaningful dependence of surface density on radius. Dent et al. (2014) do derive a nonparametric estimate for the surface density profile of the disk, which is comparable to the derived 49 Ceti profile in that it first increases and then decreases with radius, showing a local enhancement at a distance of ∼100 au from the central star. Others, including Fomalhaut, HD 181327, η Corvi, HD 21997, and a sample of the brightest debris disks in the Sco Cen region, are either intrinsically too radially narrow or have only been observed at too coarse an angular resolution (${\rm{\Delta }}R/d\lesssim \theta$, where θ is the angular resolution of the data and d is the distance to the system) for the dependence of surface density on radius to be reliably measured (Boley et al. 2012; Moór et al. 2013; Lieman-Sifry et al. 2016; Marino et al. 2016, 2017; White et al. 2016). For example, while the power-law index p describing the dependence of surface density on radius was included in models of the continuum emission from HD 21997, the uncertainty on the measured value of −2.4 ± 4.5 is too large to meaningfully distinguish between rising and falling surface density as a function of radius in that system. The other systems that have been observed to date that are intrinsically radially broad enough, axisymmetric enough, and have been observed at sufficiently high angular resolution to permit a detailed characterization of their dependence of surface density on radius include AU Mic and HD 107146.

The M-type star AU Mic, a member of the β Pic moving group with an age of ∼20 Myr, hosts an edge-on debris disk that displays no CO emission (A. M. Hughes et al. 2017, in preparation). Millimeter-wavelength observations with the SMA by Wilner et al. (2012) and with ALMA by MacGregor et al. (2013) show it to be consistent with a dust surface density increasing steeply with radius as ${\rm{\Sigma }}\propto {r}^{2.8}$ between radii of 10–40 au before being abruptly truncated. This surface density profile is consistent with predictions for a self-stirred disk with ongoing planet formation (Kennedy & Wyatt 2010), but the timescale required to initiate the collisional cascade is much longer than its age (see Kennedy & Kenyon 2008). Strubbe & Chiang (2006) suggest that a birth ring exists at a radius of ∼40 au, in which grains of all sizes are created through collisions between parent planetesimal bodies. The smaller grains, highly susceptible to stellar winds and radiation pressure, are blown out on eccentric orbits forming an extended halo visible in scattered light (Kalas et al. 2004; Krist et al. 2005; Graham et al. 2007; Schüppler et al. 2015), and exhibit time-dependent behavior (Boccaletti et al. 2015).

The 100 Myr old G star HD 107146 hosts a debris disk that has been previously spatially resolved at millimeter wavelengths with CARMA (Corder et al. 2009) and the SMA (Hughes et al. 2011). Cycle 0 observations with ALMA at an angular resolution of 1 arcsec have spatially resolved the radial surface density profile (Ricci et al. 2015). The ALMA continuum emission between radii of 30 and 150 au show an overall trend of increasing surface density with radius, similar to AU Mic. However, mirroring the surface density enhancement that we have measured in the 49 Ceti disk, there is evidence for a local decrement of surface density at a radius near 80 au. The simplest explanation for the clearing in the dust disk is that a planet of a few ${M}_{\oplus }$ has coalesced and cleared a nearly circular orbit at a distance of ∼80 au from the central star. A planet on a high-eccentricity orbit with smaller semimajor axis has also been proposed as the cause (Pearce & Wyatt 2015). The increasing surface density of the outer disk is only marginally consistent with a self-stirring scenario, requiring an initial mass surface density roughly 10 times that of the minimum mass solar nebula to initiate stirring on a timescale comparable to the age of the system.

49 Ceti's best-fit surface density power law of p = 1.29 in the single power-law model is comparable to the typical surface density power law of $p\sim 1$ observed for optically thick protoplanetary dust disks in nearby young star-forming regions (see, e.g., Andrews & Williams 2007), perhaps suggesting that 49 Ceti's dust disk is a scaled-down (and hollowed-out) remnant of the initial protoplanetary disk. In this scenario, the local enhancement at a radius of ∼110 au might correspond to the outer boundary of a self-stirred region, with the outer component perhaps composed of a significant fraction of primordial grains.

Whatever its origin, the decreasing surface density profile in the outer regions of the 49 Ceti dust disk stands in sharp contrast to the increasing surface density profiles of AU Mic and HD 107146, which are consistent with theoretical predictions for self-stirred debris disks (e.g., Kenyon & Bromley 2002; Kennedy & Wyatt 2010). While only three debris disks have yet been observed at millimeter wavelengths with sufficient sensitivity and angular resolution to provide reliable constraints on the radial surface brightness profile of continuum emission, the results so far are suggestive that 49 Ceti's gas-bearing dust disk may be undergoing a fundamentally different physical process than that of the gas-poor debris disks around AU Mic and HD 107146, resulting in the qualitatively different surface density profiles. With the new gas-bearing debris disks recently identified by Moór et al. (2015b) and Lieman-Sifry et al. (2016), it should soon be possible to determine whether or not this trend persists for a larger sample of objects.

A striking feature of the best-fit CO disk structure is that the surface density increases with distance from the central star (Section 4.2 and Table 5), independent of the density of the collision partner. The difference is statistically significant; when the surface density power law γ is fixed at zero while the other parameters are varied, the best-fit model has a ${\chi }^{2}$ value that exceeds the best-fit model with negative gamma by a difference of ${\rm{\Delta }}{\chi }^{2}=1631$, which represents $\gt 10\sigma$ significance according to the AIC. This result, along with the central gap in the CO disk, is a clear departure from the typical surface density structure of protoplanetary disks around pre-main-sequence stars, which exhibit surface density profiles that decrease with distance from the central star—very consistently in continuum emission, but also in the few cases for which the gas surface density profile has been measured using low optical depth molecular lines (e.g., Piétu et al. 2005; Panić et al. 2008; Qi et al. 2011; Flaherty et al. 2015; Schwarz et al. 2016). While increasing gas surface density profiles have been observed within the central cavity of transition disks (Carmona et al. 2014; van der Marel et al. 2016), the same disks exhibited evidence for decreasing surface density of gas in the outer disk. It should be noted that the quantity directly probed by the data is the surface brightness profile, $p+q=0$ for our assumed temperature power law of q = 0.5, but nonetheless the surface brightness is approximately constant with distance from the central star, in stark contrast with previous observations of gas disks.

At the time of writing we are not aware of any theoretical papers that consider or predict gas surface density profiles in debris disks that increase with radius, making the interpretation of this result challenging. In a steady-state scenario in which the gas is produced and destroyed at a constant rate, the CO volume density at a given radius will depend on the gas production rate and the destruction rate. In the case of negligible shielding of CO from photodissociation, the photodissociation rate should be proportional to the UV flux. In regions of the disk where the stellar contribution dominates over the interstellar component of the UV radiation field, the UV flux and therefore the photodissociation rate would dilute with radius as ${R}^{-2}$. If the gas production rate follows a power-law profile described by ${R}^{-x}$, then the final radial density distribution of CO would be described by a power law of ${R}^{2-x}$. In this scenario, our best-fit radial volume density distribution of CO in the midplane ${\rho }_{0}\propto \tfrac{{\rm{\Sigma }}(R)}{h(R)}\propto {R}^{-0.75}$ would imply a gas production rate (in units of surface density) that varies with radius as ${R}^{-1.5}$, comparable to the dust surface density in the 49 Ceti system and to a typical protoplanetary disk surface density profile. (This only applies to the case of negligible shielding. In the presence of sufficient quantities of neutral carbon shielding might change the dependence of photodissociation rate on radius; see Kral et al. 2016). If the dynamical timescale (∼130 years at a radius of 100 au) is comparable to or shorter than the CO photodissociation timescale, then the produced CO might be plausibly axisymmetrically distributed throughout the disk by Keplerian shear.

It is also important to consider whether the quantity of gas in the system is sufficient to cause noticeable effects on the dust due to interactions between the two components. There have been several theoretical investigations of the effects of gas dynamics on the solid component of debris disks. The pioneering work of Takeuchi & Artymowicz (2001) predicted that for A stars surrounded by a total gas mass ranging from a fraction to dozens of Earth masses (within which range our observations place the dust mass of the 49 Ceti disk, with the exact value depending on the assumed CO abundance relative to H₂), small dust grains (with sizes of ∼10–200 μm) will start to orbit at sub-Keplerian speeds due to Poynting–Robertson drag and subsequently experience a tailwind from the (slightly less sub-Keplerian) gas that can eject them out to orbital radii beyond the edge of the gas disk. This dust migration can result in either a ring just beyond the outer edge of the gas disk or a bimodal radial dust distribution. In a related work by Klahr & Lin (2005), they demonstrate a general result that gaps and sharp edges can form in debris dust due to instabilities within optically thin gas disks where the gas temperature is an increasing function of the dust surface density. Subsequent studies by Herrmann & Krivov (2007) and Besla & Wu (2007) added photophoresis and photoelectric effects, respectively, to the calculation. Herrmann & Krivov (2007) demonstrated that photophoresis generally increases the stability radius for a grain of a given size (as does a higher gas density, a lower absolute value of the temperature power-law index q, or higher stellar luminosity). They also point out that the outer ring-like structure described by Takeuchi & Artymowicz (2001) might not be so sharp in a disk with a less abruptly truncated outer edge. This might help to explain the more gradual decrease in dust density beyond the outer edge of the gas disk that we observe in the 49 Ceti system. The observation that the radial extent of the CO disk is smaller than that of the dust continuum emission is both in sharp contrast with observations of protoplanetary disks (e.g., Andrews et al. 2012) and consistent with an emerging trend in the radial profiles of other gas-bearing debris disks as well. Flaherty et al. (2016) and White et al. (2016) report observations showing that the CO emission from the HD 141569 disk is confined to within the smaller of the two scattered light rings, while Lieman-Sifry et al. (2016) show that for a sample of 23 debris disks in the Sco-Cen region, the three gas-bearing disks in the sample exhibit preferentially large radii for their dust continuum emission. However, these effects would require that the dust grain size distribution in 49 Ceti be such that the dominant source of opacity in the 850 μm image of the disk be grains smaller than ∼200 μm that are sensitive to Poynting–Robertson drag, which may pose additional problems for describing the surface brightness distribution (e.g., Thébault & Augereau 2005) and is likely inconsistent with the observed millimeter spectral index of the system.

Besla & Wu (2007) point out that in disks where the gas-to-dust mass ratio is $\gtrsim 1$ and the dust mass is at least ∼0.005 M_⊕, photoelectric heating effects will enhance any overdensity of dust grains (e.g., due to a recent catastrophic collision). This instability tends to create narrow rings of enhanced dust density, where the ring location is not sensitive to the radial profile of the gas disk. A similar mechanism is also described with more detailed simulations by Lyra & Kuchner (2013), who predict a clumping instability that organizes dust into narrow, eccentric rings within gas-bearing debris belts (again, in systems with gas-to-dust mass ratios $\gtrsim 1$). Such an instability might help to explain the local enhancement of dust surface density at a radius of ∼110 au in the 49 Ceti disk, perhaps indicating the recent production of dust through a destructive collision between planetesimals. Although the gas-to-dust mass ratio in the 49 Ceti system is not well constrained, it seems plausible that it might lie in the $\gtrsim 1$ regime required by the Besla & Wu (2007) and Lyra & Kuchner (2013) theories. The assumption of an ISM-like CO/H₂ ratio of 10⁻⁴ yields a global gas-to-dust mass ratio of six for the double power-law model, while an assumption that CO provides the dominant contribution to the gas mass results in a gas-to-dust mass ratio of 0.008. It is also important to consider that the local gas-to-dust mass ratio might also exceed these values (due to unresolved local vertical or radial density variations), though it is difficult to determine precisely where this might occur since the spatial resolution of our observations is insufficient to resolve the vertical structure of the gas disk.

One puzzling feature of the dust emission that might be explained by dynamical interactions between gas and dust is the apparent presence of a population of dust grains substantially smaller than the blowout size for the 49 Ceti system, previously remarked upon by Wahhaj et al. (2007). The presence of an inner disk of small (∼0.1 μm) dust grains is echoed by a recent survey for gas emission in 10 Myr old debris disks in the Sco-Cen region, which demonstrated that gas-bearing debris disks appear more likely than gas-poor debris disks to demonstrate evidence for a population of grains smaller than the blowout size, based on a simultaneous analysis of their SEDs and spatially resolved visibilities (Lieman-Sifry et al. 2016). In that survey, which included three gas-bearing debris disks, two of the three disks with evidence for a population of dust below the blowout size were gas-bearing, and the third gas-bearing debris disk exhibited a characteristic grain size only slightly larger than the blowout size. It is possible that drag forces from the substantial quantity of molecular gas in these systems is preventing the removal of small grains by radiation pressure. It is also possible that these systems have undergone more recent, active collisions that have generated both the gas and the small dust grains that have not yet been removed from the system; however, as Lieman-Sifry et al. (2016) observe, such a scenario would tend to predict that the high gas masses would be accompanied by concomitantly high dust masses, yet the survey reveals no obvious bias toward high dust mass among the gas-bearing debris disks in the sample.

The lack of a clear silicate emission feature in the Spitzer IRS spectrum of 49 Ceti is somewhat puzzling in the context of the small grain population. There are at least two plausible explanations for its absence. First, silicate features are most readily apparent when emitting grains are substantially smaller than ∼1 μm (Papoular & Pegourie 1983), though Natta et al. (2007) show that a 10 μm silicate feature is also apparent in the emission of 200–600 K grains smaller than a few microns. Our best-fit models suggest an inner radius of ∼5 au for the inner dust component, with a substantial amount of degeneracy in the absence of spatially resolved data. At this location, or if the true location of the warm component is somewhat further from the central star, the grains may not be quite hot enough to excite a silicate emission feature. Alternatively, it is possible that the grains are predominantly carbonaceous in nature rather than silicate, which would also alter the density structure of the disk. This scenario is supported by the low albedo of 49 Ceti in previous HST NICMOS observations searching for scattered light in the outer disk (Weinberger et al. 1999; Choquet et al. 2017), and by the evidence for carbon overabundance in the atomic gas disk from line-of-sight UV absorption spectroscopy (Roberge et al. 2014).

The ALMA observations also allow us to place the best constraints to date on the vertical structure of the gas disk. We find a best-fit scale height factor h₀ of 0.75 for the primordial abundance scenario, indicating that the observed scale height is substantially smaller than that expected for a typical H₂-rich disk in hydrostatic equilibrium. A scale height factor less than one cannot be easily explained by the vertically isothermal approximation (any appreciable optical depth would cause the disk surface to be hotter than the midplane, puffing up the disk beyond the height expected from the midplane temperature), nor by a scenario in which the kinetic temperature exceeds the excitation temperature (if the kinetic temperature is significantly hotter than the excitation temperature, we would expect a more, not less, vertically extended disk). A departure from LTE could be responsible for the difference only if the kinetic temperature were lower than the excitation temperature, e.g., in the case of strong UV pumping in a very low-density gas that is not thermalized by collisions with atomic or molecular hydrogen. The scale height factor should also be considered an upper limit, since the best-fit disk structure exhibits a scale height approximately equal to the spatial resolution at the largest radii of the gas disk. The simplest explanation for a scale height smaller than that predicted for hydrostatic equilibrium in the primordial abundance scenario is that the primary constituent of the gas disk is not molecular hydrogen. Since the scale height is proportional to the sound speed, which scales as $1/\sqrt{\mu }$, if the gas mass is dominated by CO instead of H₂, we would expect a scale height factor of $\sqrt{2.37/28}=0.3$, which is consistent with our upper limit for the primordial abundance scenario, and with the best-fit value of 2.1 derived for the low-H₂ scenario (which should also be considered an upper limit since it too is spatially unresolved). A small-scale height is interesting because it increases the potential for CO self-shielding along the disk midplane and might play a role in enabling longer lifetimes for volatile-rich molecular gas disks despite the energetic radiation from the central A star. A lower H₂ gas mass would also indicate that dynamical interactions between gas and dust are probably not significant, and would provide some support for a second-generation origin for the CO. Higher-resolution observations of the 49 Ceti disk are necessary to better constrain the vertical structure of the gas disk, and thereby understand the potential for CO self-shielding.

Since it is unclear whether or not there is sufficient gas to affect dust dynamics in the 49 Ceti system, we consider dust stirring mechanisms in the context of predictions for typical gas-poor debris disks, with the caveat that depending on the composition of the gas it might have an effect on the apparent properties of the dust disk. The two main mechanisms proposed for generating dust in debris disk systems are self-stirring, in which ensembles of Pluto-size planetesimals excite eccentricity (e.g., Kenyon & Bromley 2004), and stirring by a planetary-mass object (e.g., Mustill & Wyatt 2009). Other possibilities include stochastic collisions or steady-state collisional evolution of pre-stirred bodies (for an overview, see Wyatt 2008).

The self-stirring mechanism begins to operate with the formation of the first objects of ∼1000 km in scale. Kennedy & Kenyon (2008) derive a timescale of

$$\begin{eqnarray}&&t\simeq 145\,{x}_{m}^{-1.15}{\left(\displaystyle \frac{R}{80\mathrm{au}}\right)}^{3}{\left(\displaystyle \frac{2{M}_{\odot }}{{M}_{* }}\right)}^{3/2}\,[\mathrm{Myr}]\end{eqnarray} \tag{ 8 }$$

for the formation of objects sufficiently large to excite self-sirring, where R is the semimajor axis of the dust belt, ${M}_{* }$ is the mass of the central star, and x_m is a scaling factor used to describe the initial surface density of the protoplanetary disk ${{\rm{\Sigma }}}_{i}$:

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{i}(a)={{\rm{\Sigma }}}_{0}({M}_{* })\,{x}_{m}{\left(\displaystyle \frac{R}{{R}_{0}}\right)}^{-3/2}\,[{\rm{g}}\,{\mathrm{cm}}^{-2}].\end{eqnarray} \tag{ 9 }$$

The typically assumed reference surface density is the minimum mass solar nebula (MMSN; Weidenschilling 1977), for which ${{\rm{\Sigma }}}_{0}\simeq 0.18$ g cm⁻² at ${R}_{0}=30\,\mathrm{au}$. Kennedy & Kenyon (2008) assume that this surface density scales linearly with the mass of the central star, such that ${{\rm{\Sigma }}}_{30\mathrm{au}}({M}_{* })\simeq 0.18\,\left(\tfrac{{M}_{* }}{{M}_{\odot }}\right)$ g cm⁻². Solving for the scaling factor needed to produce a ring at 110 au within the 40 Myr age of 49 Ceti yields a factor ${x}_{m}\simeq 8$, or an initial surface density of solids around 49 Ceti of roughly 16×  the MMSN. Mustill & Wyatt (2009) demonstrate that x_m must be $\lt 10$ to avoid gravitational instability for a canonical 100:1 gas-to-dust mass ratio, implying that 49 Ceti's initial protoplanetary disk would need to be close to the point of gravitational instability in order for self-stirring to be a viable mechanism to produce the observed surface density enhancement.

The ability of a planet to excite dust production through its gravitational influence generally depends on its mass (M_pl), orbital eccentricity (e_pl) and semimajor axis (R_pl), as well as the mass of the central star (${M}_{* }$). The required relative excitation velocity for destructive impacts between equal-size bodies in the outer disk, where weak conglomerations of ices tend to dominate, is derived by Mustill & Wyatt (2009) as

$$\begin{eqnarray}&&{v}_{\mathrm{rel}}(a)={\left[0.8{\left(\displaystyle \frac{a}{80{\rm{m}}}\right)}^{-0.33}+0.2{\left(\displaystyle \frac{a}{80{\rm{m}}}\right)}^{1.2}\right]}^{0.83}\,[{\rm{m}}\,{{\rm{s}}}^{-1}]\end{eqnarray} \tag{ 10 }$$

where a is the size of the body. Destructive collisions between grains of size $a\sim 1\mbox{--}10\,\mathrm{cm}$ can produce the 100 μm–1 mm particles that are thought to dominate the emission at the wavelength of the ALMA image, though larger bodies such as kilometer-size comets may also be the parent bodies of the observed dust (and gas)—either parent body requires collisional velocities of order 5–10 m s⁻¹ to result in destructive collisions. For an internal perturber, the semimajor axis R* within which a planet can excite grains of radius a to v_rel is derived by Mustill & Wyatt (2009) as

$$\begin{eqnarray}{R}^{* }(a) & = & 3.8\,{\left(\displaystyle \frac{{e}_{\mathrm{pl}}}{0.1}\right)}^{2/3}{\left(\displaystyle \frac{{M}_{* }}{1{M}_{\odot }}\right)}^{1/3}\\ & & \times {\left(\displaystyle \frac{{R}_{\mathrm{pl}}}{1\mathrm{au}}\right)}^{2/3}{\left(\displaystyle \frac{{v}_{\mathrm{rel}}(a)}{1\mathrm{km}{{\rm{s}}}^{-1}}\right)}^{-2/3}\,[\mathrm{au}]\end{eqnarray} \tag{ 11 }$$

while the timescale for stirring to occur is given by

$$\begin{eqnarray}t & \sim & 1.53\times {10}^{3}\displaystyle \frac{{(1-{e}_{\mathrm{pl}}^{2})}^{3/2}}{{e}_{\mathrm{pl}}}{\left(\displaystyle \frac{R}{10\mathrm{au}}\right)}^{9/2}\\ & & \times {\left(\displaystyle \frac{{M}_{* }}{1{M}_{\odot }}\right)}^{-1}{\left(\displaystyle \frac{{M}_{\mathrm{pl}}}{{M}_{\odot }}\right)}^{-1}{\left(\displaystyle \frac{{R}_{\mathrm{pl}}}{1\mathrm{au}}\right)}^{-3}\,[\mathrm{years}].\end{eqnarray} \tag{ 12 }$$

Figure 8 displays the eccentricity required for a planet of semimajor axis R_pl and mass M_pl to stir the planetesimal belt at 110 au in less than 49 Ceti's age of 40 Myr. For example, a Jupiter-like planet with ${e}_{\mathrm{pl}}\sim 0.05$ and ${M}_{\mathrm{pl}}\sim {10}^{-3}{M}_{\odot }$ located at a distance of 40 au from the central star could excite planetesimals to collisional velocities within the dust ring in less than the age of the system.

Zoom In Zoom Out Reset image size

49 Ceti joins a handful of gas-bearing debris disks that have been imaged with ALMA, and therefore have high-quality images of both the molecular gas disk and thermal dust continuum emission: the others are HD 21997 (Kóspál et al. 2013; Moór et al. 2013), β Pictoris (Dent et al. 2014), HD 141569 (White et al. 2016), HD 181327 (Marino et al. 2016), HIP 76310, HIP 84881, and HIP 73145 (Lieman-Sifry et al. 2016). Among these systems, 49 Ceti is the first with sufficiently axisymmetric structure and sufficiently high-resolution and high-sensitivity observations that the surface density profile of both gas and dust can be studied in detail. However, previous observations have provided clues to the origin of the gas: the investigations center around whether the gas is remnant primordial material that has somehow managed to survive despite the destructive influence of high-energy radiation from the central star, or whether the gas is second-generation material that was previously incorporated into planetary embryos, planetesimals, or the icy mantles of dust grains.

Structurally, β Pictoris is unique among the disks imaged in exhibiting a very strong asymmetry of both gas and dust emission, with a large clump of emission appearing along one limb of the disk (though since the β Pic disk is viewed nearly edge-on to the line of sight, it is impossible to uniquely deproject the radial structure of the feature). The structural asymmetries, along with the surprising presence of a number of atomic species revealed by absorption spectroscopy, provided strong hints that the gas and dust were likely of secondary origin, resulting from a recent collision between Mars-size bodies or collisional evaporation of icy planetesimals at resonances induced by the presence of a giant planet (Dent et al. 2014). A recent analysis of optically thin CO line ratios from the system determines conclusively that the gas must be of second-generation origin, thanks to a non-LTE analysis demonstrating that the H₂ densities must be too low to shield CO from photodissociation, and disfavors the collision of Mars-size bodies as the origin of the gas (Matrà et al. 2017). However, it is not clear to what extent β Pictoris is typical or atypical compared to other gas-bearing debris disks, since all the other systems have shown at most minor deviations from axisymmetry at least at the ∼10% level enabled by the contrast ratios of current observations and are predominantly characterized by ring- or disk-like structures in both gas and dust.

Among the more regular rings and disks observed around the other stars, arguments about the origin of the gas tend to focus on the following characteristics.

Relative location and structure of gas and dust—An initial clue to the physical mechanisms underlying the origin of gas and dust comes from their relative locations in the disk. If both the gas and dust are second-generation, produced in similar events involving the collisional destruction of planetesimals, then one might expect their physical locations in the disk to be similar. If, on the other hand, the gas is primordial while the dust is second-generation, the gas might follow the radial profile of the primordial protoplanetary disk while the dust exhibits ring-like structures more characteristic of an archetypal debris disk. Muddying the waters a bit, both configurations—and in 49 Ceti, a configuration that does not obviously match either scenario—have already been observed even in the small sample of disks with both thermal continuum and molecular gas emission images. Marino et al. (2016) note that the CO emission, while faint and difficult to image, appears to follow the ring-like structure observed in dust continuum emission in the disk around HD 181327, and cite this co-location as potential evidence for a second-generation origin. By contrast, Kóspál et al. (2013) observe a far more radially extended CO configuration than the dust continuum ring observed in the HD 21997 disk, and cite this difference as evidence that the gas may be primordial while the dust is second-generation in that system. 49 Ceti's gas and dust configurations exhibit different radial dependences as well as different inner and outer extents, perhaps suggesting that their physical origins differ as well—though it is also suggestive to consider the possibility that the underlying gas production rate is diluted by the ${R}^{-2}$ UV flux, implying that the underlying surface density profile might in fact match that of the dust disk—and that of a typical protoplanetary disk rather than a typical self-stirred debris disk. Similarly high-resolution and high-sensitivity observations of the other gas-bearing debris disks will provide a more thorough understanding of the relative location of gas and dust in the known gas-bearing debris disk systems.

Molecular gas densities and lifetimes—Another cornerstone of the discussion has been the true quantity of molecular gas in the disk, which centers around the question of the gas composition, and particularly the CO-to-H₂ abundance and whether or not it differs from the primordial value of 10⁻⁴. Interstellar abundances of H₂ might be capable of shielding CO from destruction by the harsh stellar UV radiation fields, implying that primordial CO might be able to survive to the 10–40 Myr ages of the stars surrounded by substantial molecular gas disks. In the absence of a substantial H₂ population, by contrast, the CO photodissociation timescales are quite rapid in these systems, of the order of hundreds to thousands of years depending on the details of the spatial distribution of CO (Visser et al. 2009). Circumstantial evidence in several systems has implied that the H₂ content of the disks may be quite low, including low excitation temperatures derived from multi-transition and/or multi-transition analysis of CO that imply a lack of collision partners for the CO molecules (Kóspál et al. 2013; Flaherty et al. 2016; Matrà et al. 2017). There is also more conclusive evidence of a non-primordial CO/H₂ ratio in both β Pic and AU Mic (Lecavelier des Etangs et al. 2001; Roberge et al. 2005). While the excitation temperature in the 49 Ceti disk is not obviously lower than the kinetic temperature, the spatially unresolved vertical structure in this nearly edge-on system provides circumstantial evidence of a mean molecular weight substantially higher than the ISM value of ∼2.37, lending support to the non-ISM CO/H₂ abundance and therefore the second-generation origin of material in the disk. While the apparent trend of dust disks larger than their CO counterparts might be neatly explained by the presence of gas in sufficient quantities to create the tailwind effect predicted by Takeuchi & Artymowicz (2001), this would also require that the dust grain size distribution be sufficiently weighted toward small ($\lt 200\,\mu {\rm{m}}$) grains that they dominate the opacity even at millimeter wavelengths.

Required production rates for secondary material—One of the most significant challenges to the second-generation gas disks scenario is the question of whether it is plausible to sustain collision rates high enough to maintain the gas disk for a substantial fraction of the star's lifetime. Standard calculations based on the rapid photodissociation timescales of hundreds to thousands of years suggest that it would be necessary to evaporate a large (Hale-Bopp-size) comet every few minutes to sustain the necessary gas production rates (Zuckerman & Song 2012; Kóspál et al. 2013; Dent et al. 2014; White et al. 2016). Whether or not such collision rates are reasonable given the likely reservoirs of cometary material populating the collisional cascade is still a subject of discussion within the literature. Kóspál et al. (2013) argue that the rate is at best marginally consistent with expected reservoir, while Zuckerman & Song (2012) argue that such a high cometary collision rate is consistent with the high dust luminosity of the HD 21997 system. Dent et al. (2014) calculate the collision rate without commenting on its plausibility. White et al. (2016) calculate that there is potentially sufficient cometary material in the HD 141569 disk to produce the observed CO flux, though significant shielding is required to maintain the observed abundance for a significant length of time. Some disks do exhibit demonstrable consistency with the mass-loss rate produced by a steady-state collisional cascade based on solar system cometary abundances (Matrà et al. 2015, 2017; Marino et al. 2016). Chance collisions of larger, Mars-size bodies in the disk have the potential to release larger quantities of CO, but such statistically rare events are unlikely to explain the population of "typical" axisymmetric gas-bearing debris disks given observation that approximately half of intermediate-mass debris disk host stars with ages of ∼10 Myr harbor gas-bearing debris disks (Lieman-Sifry et al. 2016). A recent model suggests that evaporation from grains in the collisional cascade may be able to explain the quantity of gas observed in the 49 Ceti system (Kral et al. 2017), although this will need to be reconciled with the lack of correlation between dust and gas mass in debris disks observed by Lieman-Sifry et al. (2016). The best unexplored avenue for disentangling a cometary from a primordial gas disk origin is by examining the chemical composition of the molecular gas component. While HCN seemed a promising molecule due to its relatively high abundance both in typical protoplanetary disks (Öberg et al. 2010) and its order-of-magnitude higher abundance relative to CO in solar system comets (Biver et al. 2002), unfortunately the search for HCN in the current study of 49 Ceti was fruitless. It is difficult to predict how the composition of second-generation material produced through collisions at large distances from the star might compare to that observed for photodesorption of ices on comets in the inner solar system. However, other molecules that are not as readily photodissociated may be more likely to yield useful chemical insights that can distinguish the composition of primordial gas from that of comet-like material.

At the moment, the most promising future avenues for disentangling the origin of the gas in debris disks would appear to be a non-LTE multi-line CO analysis of a larger sample of disks, along the lines of that presented in Matrà et al. (2017), perhaps combined with a chemical search for targeted molecular lines that might help to distinguish the composition of protoplanetary gaseous material from that of cometary material. Combining these insights with measurements of the detailed spatial distribution of gas and dust in these systems will yield insight into the physical origins of the constituent material.