Detection of Exocometary CO within the 440 Myr Old Fomalhaut Belt: A Similar CO+CO₂ Ice Abundance in Exocomets and Solar System Comets

Figure 1 (top right) shows the CO J = 2-1 moment-0 map of the Fomalhaut ring, spectrally integrated within ±3.5 km s⁻¹ around the expected radial velocity of the star (${v}_{\star ,\mathrm{LSR}}=5.8\,\pm 0.5$ km s⁻¹; Gontcharov 2006). No statistically significant emission is observed above a noise distribution that is well approximated by a Gaussian and that has an rms of 4 mJy km s⁻¹ beam⁻¹. The few 3–4σ peaks observed are in line with the expectation from the Gaussian noise distribution and the large number of resolution elements in the image. The contour lines indicate the sky region where continuum emission is detected above the 4σ level (MacGregor et al. 2017), corresponding to an area comprising a number of beams of ${N}_{\mathrm{beams}}\sim 40$. As done for previous Band 7 CO observations of the Fomalhaut system (Matrà et al. 2015), we now proceed by assuming that any CO present in the system should be of secondary origin and colocated with the dust from which it is produced, and later on we test these assumptions (Sections 3.5 and 4.1). This allows us to apply a spatial filter, that is, to spatially integrate the CO data cube over this region where the continuum is detected, yielding an improvement on the signal-to-noise ratio (S/N) equal to $\sqrt{{N}_{\mathrm{beams}}}\sim 6.4$. The obtained 1D spectrum, shown in Figure 1 (left), presents a 3.3σ peak at a velocity consistent with the expected ${v}_{\star ,\mathrm{LSR}}$.

Zoom In Zoom Out Reset image size

In order to test the robustness and improve the significance of the detection, we make the further assumption that the CO gas lies in a vertically flat disk in Keplerian rotation around a star of 1.92 M_⊙, with the same orbital parameters as the dust ring, obtained from dust continuum imaging (MacGregor et al. 2017) and consistent with scattered light imaging (Kalas et al. 2013). Then, for each pixel location $({x}_{\mathrm{sky}},{y}_{\mathrm{sky}})$, where $(0,0)$ is the accurately known location of the star (detected in the continuum data set), we can determine the expected radial velocity field ${v}_{z}({x}_{\mathrm{sky}},{y}_{\mathrm{sky}})$ of the gas, both for the case in which the NW ansa is moving away from us (Figure 1, center right) or toward us (Figure 1, bottom right). Using this predicted velocity field, we apply a spectral filter; that is we shift the 1D spectrum at each $({x}_{\mathrm{sky}},{y}_{\mathrm{sky}})$ pixel location by $-{v}_{z}({x}_{\mathrm{sky}},{y}_{\mathrm{sky}})$, moving any CO emission from the predicted disk velocity to the stellar velocity ${v}_{\star ,\mathrm{LSR}}$. If we then sum together the contributions from all pixels where the dust continuum is detected (or in other words, if we apply the spatial filter), we obtain a spectrum that is free from noise originating from spectral channels as well as spatial locations where no CO emission is expected. As shown in Figure 1 (center left, for the velocity field in center right), this method boosts the S/N of the peak to $5.4\sigma$.

The spectral filter yields kinematic information on the orbiting gas, since the S/N of the CO line is maximized only for the correct velocity field. This is because, for an incorrect velocity field, the spectral filter spreads CO emission over a larger number of velocity channels, reducing the S/N (see Figure 1, bottom left and right). This proves that CO gas in the NW ansa is moving away from us, while gas in the southeast (SE) ansa is coming toward us (Figure 1, center right). This is consistent with the kinematics of the star, whose SE part was also found to be moving toward us (Le Bouquin et al. 2009). We note that this information is still insufficient to derive the sense of rotation of the ring on-sky (clockwise or counterclockwise), and in turn which side of the ring is closer to Earth and forward scattering (see discussion in Min et al. 2010; Kalas et al. 2013). This is because both the star and the CO only inform us on the sky-projected rotation axis, which is perpendicular to the sky-projected major axis of the ring and points in the northeast (NE) direction. However, the crucial missing piece of information is whether the ring's rotation axis is pointing toward us (i.e., inclination of $\sim +65^\circ$ to the plane of the sky) or away from us ($\sim -65^\circ$ inclination). Therefore, the current data are insufficient to determine whether the ring is rotating in a clockwise or counterclockwise direction on-sky, and whether the brighter NE side observed by HST is in front of or behind the sky plane.

In order to measure the integrated CO J = 2-1 line flux and the velocity centroid, we fit a Gaussian to the spectrospatially filtered spectrum derived above and shown in Figure 1 (center left). The best-fit Gaussian width is consistent with the line being unresolved, as expected from the application of our filtering method, and the best-fit velocity is 6.1 ± 0.2 km s⁻¹, consistent with the expected stellar velocity of 5.8 ±0.5 km s⁻¹. The best-fit integrated line flux is 68 ±16 mJy km s⁻¹ (including the absolute flux calibration uncertainty, added in quadrature to the uncertainty from the Gaussian fit), or $(5.2\pm 1.2)\times {10}^{-22}$ Wm⁻².

Assuming that the CO is optically thin to the line of sight at the observed frequency (which we verify in the next paragraph), we use a nonlocal thermodynamic equilibrium molecular excitation analysis (Matrà et al. 2015) to derive constraints on the total CO mass in the system. The derived mass value depends on two unknown parameters, the kinetic temperature of the gas and the density of the main collisional partners. We take the main colliders to be electrons, since these are most likely to be the dominant species for which CO collision rates are known if the gas is of exocometary origin (Matrà et al. 2015). We cover the full electron density parameter space between the radiation-dominated and local thermodynamic equilibrium regimes, and we probe a wide range of kinetic temperatures between 10 and 1000 K. The CO mass derived is in the range $(0.65\mbox{--}42)\times {10}^{-7}\,{M}_{\oplus }$, where the boundaries were obtained from the $\pm 1\sigma$ limits on the integrated line flux. This is the lowest CO mass detected in any circumstellar disk to date, which is readily understood by noting that Fomalhaut, at a distance of 7.7 pc, is the nearest circumstellar disk where CO has been searched for by ALMA to date. We will discuss the implication of this mass measurement in Section 4.1.

We note that the CO excitation model does not yet account for the effect of infrared or UV pumping, that is, transitions between vibrational and electronic levels within the molecule. These are excited in the presence of a strong infrared or UV radiation field, as may be the case around an A star such as Fomalhaut. As molecules relax to the ground electronic or vibrational levels, higher rotational levels may be more populated than predicted through CMB radiation alone, affecting the rotational level populations. Nonetheless, this would only influence the molecule in the radiation-dominated regime, or in other words it would only change our upper limit on the derived range of CO masses. Since less mass would be needed to produce the same flux in the presence of significant UV/IR pumping, our derived CO mass range can be taken as a conservative estimate and is likely narrower than derived above.

We now carry out a check to probe the optical thickness of the CO line. We assume that the CO density is uniform in a disk with inner radius at $R=136.3\,\mathrm{au}$ and ${\rm{\Delta }}R=13.5\,\mathrm{au}$ wide, with a constant aspect ratio h = H/R (height above the midplane divided by the radius) equal to the best-fit mean proper eccentricity of the planetesimals from continuum observations (0.06, MacGregor et al. 2017). This way, we can derive the ring volume and hence measure an average CO number density for our range of CO masses, obtaining a range of number densities in the range $(2\mbox{--}75)\times {10}^{-2}$ cm⁻³. We then estimate the maximum column density using the longest path length along the line of sight passing through our simple model ring described above. We neglect any column density enhancement expected at the two ansae due to projection effects (see Section 4.6). For the range of number densities above, this yields maximum column densities of $(0.04\mbox{--}1.67)\times {10}^{14}$ cm⁻² for a maximum line-of-sight path length of ∼15 au. Finally, we calculate the maximum optical thickness through its definition (Equation (3) in Matrà et al. 2017) for the full range of electron densities and kinetic temperatures probed above, and for an intrinsic line width taken to be equal to the Doppler broadening expected from each of the temperatures considered. The maximum, worst-case-scenario optical thickness we obtain is $\tau \lesssim 0.18$ (with values ranging down to 10⁻⁴ depending on the parameters), which confirms our optically thin assumption and our CO mass measurement.

We now check that our detection of the CO J = 2-1 line is consistent with the CO J = 3-2 upper limit available from previous ALMA Band 7 observations (Matrà et al. 2015). For a well-characterized millimeter radiation field such as that around Fomalhaut (see Figure 3 in Matrà et al. 2015), line ratios of optically thin molecular transitions only depend on two free parameters determining the level populations of the upper level of each transition. These are the gas kinetic temperature ${T}_{{\rm{k}}}$ and the density of collisional partners ${n}_{{{\rm{e}}}^{-}}$. Taking our J = 2-1 integrated line flux of $(5.2\pm 1.2)\times {10}^{-22}$ Wm⁻² and the 3σ upper limit of 18 × 10⁻²² Wm⁻² on the J = 3-2 line flux from Matrà et al. (2015), we obtain an upper limit on the average 3-2/2-1 line ratio in the Fomalhaut ring of 3.5 or 6.7, depending on whether we consider our J = 2-1 measurement or its 3σ lower limit. This line ratio upper limit of 3.5 in Fomalhaut will therefore trace a line in ${T}_{{\rm{k}}}$–${n}_{{{\rm{e}}}^{-}}$ space (Figure 2), the latter two quantities being degenerate. For comparison, we show the average line ratio measured in the β Pictoris disk (Matrà et al. 2017). While strictly speaking we can only exclude 3-2/2-1 line ratios higher than 6.7, we find that the archival nondetection of the J = 3-2 transition is in agreement with a β Pictoris-like gaseous environment (average line ratio of 1.9 ± 0.3) and fully consistent with the new J = 2-1 detection. While this is purely a consistency check, we remind the reader that there is no reason to assume that the electron density or CO excitation conditions in the Fomalhaut belt should be the same as around β Pictoris. Furthermore, we note once again that the introduction of UV/IR pumping in our model may influence Figure 2 by increasing the minimum possible line ratio in the radiation-dominated regime of excitation (left-hand side in the figure).

Zoom In Zoom Out Reset image size

In order to extract information on the spatial distribution of CO, we relax the assumption that it must be colocated with the dust millimeter continuum. In particular, we apply the spectral filter method to the entire data cube, shifting spectra in each spatial pixel to align CO emission with the stellar velocity. Instead of spatially integrating across the area where the continuum is detected, we examine the channel map corresponding to the stellar velocity (as shown in Figure 3). As expected, due to spatial dilution of the emission, no obvious significant emission is seen along the whole dust ring. For a given semimajor axis in the ring's orbital plane, taking once again the orbital parameters from dust continuum fitting, we can obtain a sky-projected orbit, and we average all azimuthal contributions along it to obtain an intensity profile as a function of ring semimajor axis. The black line in Figure 4 shows the semimajor axis profile between 50 and 220 au obtained by azimuthally averaging around all true anomalies, showing a CO detection at the semimajor axis and width consistent with the dust continuum ring. This confirms that CO is indeed colocated with dust in the Fomalhaut belt.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Finally, we test for any indication of azimuthal asymmetry in the CO emission. To begin with, we construct the same semimajor axis profile for each of six azimuthal regions, or "wedges," as displayed through different colors in Figure 3. These profiles are shown by the corresponding colors in Figure 4. For a perfectly axisymmetric on-sky surface brightness distribution of the ring, due to averaging across only about one-sixth of the azimuths, we expect the S/N in each wedge to decrease by a factor of $\sim \sqrt{6}=2.45$ compared to the profile averaged across all azimuths and to lie at the $1\mbox{--}2\sigma$ level. However, the majority of CO emission appears to originate from the south (S) wedge (blue), where CO is detected at the $6.2\sigma$ level. The latter is also apparent in Figure 3, where we see that this region of emission is located near the ring's pericenter location, as determined by ALMA dust continuum fitting ($\omega =22\buildrel{\circ}\over{.} 5\pm 4\buildrel{\circ}\over{.} 3$, black cross in Figure 3; MacGregor et al. 2017). We find that the flux in the S wedge contributes to $(49\pm 27)$% of the flux integrated across all azimuths; given the low S/N levels, we cannot distinguish whether all of the emission originates from this wedge or whether other regions also contain CO emission, as however hinted at by the low levels of positive emission observed across all wedges.

To analyze this azimuthal variation in more detail, Figure 5 shows the azimuthal profile obtained by integrating emission between semimajor axes of 120 and 150 au for all true anomalies f. We recover the enhancement at pericenter (true anomaly $f=0^\circ$) and another $2.6\sigma$ peak of emission at $f\sim 45^\circ$, also evident in Figure 3. We then fit a cosine function to this azimuthal profile; this is a good representation of the expectation from our steady-state release model discussed in Section 4.3, for the case where a pericenter enhancement of the CO mass with respect to apocenter is predicted. We use Levenberg–Marquardt least-squares minimization to find the best-fit phase, mean intensity, and pericenter-to-apocenter intensity enhancement ($1-{I}_{\mathrm{apo}}/{I}_{\mathrm{peri}}$) of our model cosine function. The phase obtained with respect to pericenter ($20^\circ \pm 25^\circ$) is consistent with the model having a maximum near pericenter. The mean (0.68 ± 0.16 mJy beam⁻¹ per 1.3 km s⁻¹ channel) being significantly different from zero is further confirmation of our detection at a level consistent with the radial profile in Figure 4, bottom, and with the spectrum in Figure 1, center left. This cosine fit leads to an estimate of a pericenter enhancement with respect to apocenter of $(88\pm 25)$%. We will compare such an enhancement to model predictions in Section 4.6.

Zoom In Zoom Out Reset image size

The origin of gas remains to be found for most of the known gas-bearing debris disks. A primordial versus secondary origin dichotomy has emerged in the past years, due to the youth ($\lesssim 40\,\mathrm{Myr}$) of the detected systems (see Section 1). Aside from a tentative detection of CO emission in the 1–2 Gyr old η Corvi system (Marino et al. 2017), which is not colocated with the outer dust belt and remains to be confirmed, Fomalhaut is the most evolved debris belt to host CO gas, at an age of (440 ± 40) Myr.

Can the observed CO have survived since the protoplanetary phase of evolution? CO survival requires shielding from the interstellar UV radiation field, which otherwise rapidly photodissociates the molecule on a timescale ${t}_{\mathrm{phd}}$ of 120 years. This shielding can be produced by CO itself and other molecules such as H₂, where the latter dominates the gas mass in the primordial origin scenario. In Section 3.3, through a simplified model, we estimated an average number density of CO in the ring of order $(2\mbox{--}75)\times {10}^{-2}$ cm⁻³. In this model, a CO molecule sitting in the radial and vertical center of the ring will therefore have a CO column density of order ${10}^{12}\mbox{--}{10}^{14}$ cm⁻² surrounding it and, assuming a low CO/H₂ abundance ratio of 10⁻⁶ (similar to that found in the old TW Hya protoplanetary disk; Favre et al. 2013), a H₂ column density of order ${10}^{18}\mbox{--}{10}^{20}$ cm⁻². Since these CO and H₂ column densities are insufficient to shield CO over the system's lifetime (Visser et al. 2009), we conclude that the observed CO cannot have survived since the protoplanetary phase. We note that freeze-out onto grains is also negligible, due to the relatively low density of grains in the ring and their temperature being much above the CO freeze-out temperature (Matrà et al. 2015). The observed CO must therefore be of secondary origin, that is, recently produced through either continuous, steady-state replenishment or a recent stochastic event. We analyze both possibilities below.

For a stochastic collision, the requirement is recent production of the observed CO mass, at least $6.5\times {10}^{-8}\,{M}_{\oplus }$. Assuming a 10% CO+CO₂ mass fraction, this requires the destruction of a comet of total mass $6.5\times {10}^{-7}\,{M}_{\oplus }$, or about 300 Hale–Bopp masses (e.g., Weissman 2007). Given the observationally well constrained mass loss rate of small grains through the collisional cascade (Appendix B) and the known CO mass and photodestruction rate (previous section), we can estimate (see next section) the collisional mass loss rate of such large, CO+CO₂-rich bodies in the cascade, obtaining a range of $0.012\mbox{--}0.046\,{M}_{\oplus }$ Myr⁻¹. Then, we can estimate that the timescale for collisions between any such supercomets to take place is 14–54 years, meaning that the rate is 2.2–8.6 every 120 years (the survival timescale of CO). This would suggest that it is possible that the recent destruction of a large body within the belt, alone, produced all of the observed CO mass.

However, this conclusion is subject to the assumption that bodies of the required $6.5\times {10}^{-7}\,{M}_{\oplus }$ mass for CO release, which would be ∼100 km in size, lose mass at the same rate as other bodies participating in the collisional cascade. This implicitly assumes an extrapolation of the size distribution from small, observable grains up to bodies of this size. Although ∼100 km is consistent with observations of the Kuiper and asteroid belts in the solar system (Bottke et al. 2005; Fraser et al. 2014), as well as other planetesimal growth models (Johansen et al. 2015; Shannon et al. 2016), it remains to be determined whether such large bodies in the Fomalhaut belt exist and participate in the collisional cascade. If we extrapolate the power-law size distribution of Wyatt & Dent (2002) to large sizes ($n(D)\propto {D}^{2-3q}$ with $q=11/6$, Figure 8 in their paper), we find that there should be 7 × 10⁷ objects with this or larger mass within the Fomalhaut ring. Then, the collision timescale of one supercomet is $\sim 1.0\mbox{--}3.8\,\mathrm{Gyr}$, which is longer than the ∼440 Myr age of the system. This indicates supercomets would have seldom collided over the age of the system, and their size distribution would therefore have been set by growth processes during planet formation, where the latter is completely unconstrained observationally. In addition, if the extrapolation of the size distribution to these sizes were valid, the total mass of the belt (∼0.4 Jupiter masses) would be ∼4 times higher than the 29 M_⊕ that an initial minimum-mass solar nebula (MMSN) planetesimal disk would have had between 120 and 150 au (Kenyon & Bromley 2008). In general, these required high belt masses challenge the validity of our extrapolation and point to a likely steeper size distribution for the large primordial bodies (as is the case for Kuiper Belt objects, e.g., Schlichting et al. 2013).

Overall, this caveat would make our estimated high collision rate for large bodies an upper limit. To conclude, it is therefore possible that the destruction of a large icy body released all the CO observed in the Fomalhaut ring, although a reasonable likelihood for this event to happen requires the Fomalhaut belt to be massive, of order of a few times higher than the expectation from an MMSN-like disk.

On the other hand, we can consider the total gas release expected from the steady-state collisional cascade in the framework described in Matrà et al. (2015) and already applied to other debris disks hosting secondary gas (see Section 4.4). Regardless of the details of the ice-removal mechanism, this method can be used to estimate the CO+CO₂ ice mass fraction in Fomalhaut's exocomets that is required to produce the observed CO gas. In summary, we assume that a steady-state collisional cascade is in place within the ring, with large parent bodies grinding down to produce dust of sizes all the way down to the blow-out limit (e.g., Wyatt & Dent 2002).

Solid mass is being inputted through catastrophic collisions of the largest comets in the collisional cascade, and the CO+CO₂ mass fraction will be lost through gas release within the cascade. The condition of steady state imposes that the rate at which mass is being inputted by the largest bodies $({\dot{M}}_{{D}_{\max }})$ is equal to the sum of the rate at which mass is being lost through CO+CO₂ outgassing (${\dot{M}}_{\mathrm{CO}+{\mathrm{CO}}_{2}}$) and the rate at which mass is being lost through radiation forces at the bottom of the cascade, ${\dot{M}}_{{D}_{\min }}.$ Assuming that all of the CO+CO₂ ice is lost through the cascade before reaching the smallest-sized grains (see discussion in next section), we have ${\dot{M}}_{\mathrm{CO}+{\mathrm{CO}}_{2}}={f}_{\mathrm{CO}+{\mathrm{CO}}_{2}}{\dot{M}}_{{D}_{\max }},$ meaning that we can measure the CO+CO₂ ice mass fraction in exocomets ${f}_{\mathrm{CO}+{\mathrm{CO}}_{2}}$ through

$$\begin{eqnarray}&&{f}_{\mathrm{CO}+{\mathrm{CO}}_{2}}=\displaystyle \frac{1}{1+{\dot{M}}_{{D}_{\min }}/{\dot{M}}_{\mathrm{CO}+{\mathrm{CO}}_{2}}}.\end{eqnarray} \tag{ 1 }$$

In a steady state, the outgassing rate of CO+CO₂ molecules equals their destruction rate (where the latter is known) through ${\dot{M}}_{\mathrm{CO}+{\mathrm{CO}}_{2}}={M}_{{\mathrm{CO}}_{\mathrm{obs}}}/{t}_{\mathrm{phd}}$. Additionally, the mass loss rate ${\dot{M}}_{{D}_{\min }}$ of CO+CO₂-free particles at the small size end of the collisional cascade can be estimated by considering the collision timescale of particles just above the minimum size in the distribution (Appendix B). The latter is well constrained by observations and can be calculated as shown in Equation (21).

Overall, the CO+CO₂ ice mass fraction of exocomets can therefore be estimated through

$$\begin{eqnarray}&&{f}_{\mathrm{CO}+{\mathrm{CO}}_{2}}=\displaystyle \frac{1}{1+0.0012\ {R}^{1.5}{\rm{\Delta }}{R}^{-1}{f}^{2}{L}_{\star }{M}_{\star }^{-0.5}{t}_{\mathrm{phd}}{M}_{{\mathrm{CO}}_{\mathrm{obs}}}^{-1}},\end{eqnarray} \tag{ 2 }$$

where R and ${\rm{\Delta }}R$ are in au, L_⋆ and M_⋆ are in L_⊙ and ${M}_{\odot },{t}_{\mathrm{phd}}$ is in years, and ${M}_{{\mathrm{CO}}_{\mathrm{obs}}}$ in M_⊕. The CO+CO₂ ice mass fraction in Fomalhaut is therefore in the range 4.6%–76% (taking the observed belt parameters quoted in Appendix B).

Our compositional measurement in exocomets within the Fomalhaut ring adds to the other two measurements of the CO+CO₂ exocometary ice mass fraction obtained through ALMA observations of second-generation CO gas, β Pictoris (Matrà et al. 2017) and HD 181327 (Marino et al. 2016). For consistency, we also apply our updated, more accurate estimation of the ice mass fraction presented here to these systems, using the same belt and stellar parameters as in the original works. We find that the CO+CO₂ ice mass fraction becomes $\lt 32 \%$ for β Pictoris (using the 3σ upper limit on the CO mass) and in the range 0.8%–3.7% for HD 181327 (using the $\pm 1\sigma$ range of the CO mass). These two systems both belong to the β Pic moving group at an estimated age of 23 ± 3 Myr (Mamajek & Bell 2014) and are much younger than Fomalhaut, but they nonetheless present compositions that are within an order of magnitude of one another (Figure 6). For example, the weak detection in Fomalhaut compared to β Pic is easily explained by the fact that Fomalhaut is much less collisionally active. Indeed, the mass loss rate through the cascade is $\gtrsim 20$ times higher for the β Pictoris disk, and the higher amounts of CO also mean that some self-shielding takes place, increasing the photodissociation timescale by a factor of ∼2.5 (Matrà et al. 2017). Overall, for the same CO+CO₂ ice content, this indicates that the CO gas mass would be at least 50 times higher in β Pic with respect to Fomalhaut, in agreement with the observations. On the other hand, we note that our new Fomalhaut measurement is at least marginally inconsistent with that in the HD 181327 debris ring, potentially suggesting an intrinsic difference in composition or gas release mechanism between these two systems.

Zoom In Zoom Out Reset image size

Despite the rather large error bars on the measurements obtained so far, Figure 6 also shows that the CO+CO₂ mass fraction in exocomets is consistent with that of solar system comets. A caveat to keep in mind in this comparison is that mass fractions in solar system comets are derived from measurements of the relative abundance of CO and CO₂ compared to H₂O. On the other hand, their refractory-to-volatile mass ratio, required to derive the CO+CO₂ mass fraction, is only poorly (if at all) known, due to the difficulty in remote measurements of dust production rates. The only comet for which this is measured robustly in situ is comet 67P/Churyumov–Gerasimenko (67P), where the dust-to-ice mass ratio is 4 ± 2 in the coma (Rotundi et al. 2015) and consistent with (though not directly constrained by) density measurements of the interior of its nucleus (Pätzold et al. 2016). We therefore assumed all solar system comets to have an equal refractory-to-volatile mass ratio of 4, though note that this may be subject to scatter across comet families and individual objects.

Another underlying assumption is that we are using observation-based estimates of the gas production (outgassing) rates of CO and CO₂ (as well as H₂O for solar system comets) as a probe for their relative solid abundance; in other words, we are assuming that CO+CO₂ ice and refractories are removed at rates that reflect their internal composition. In solar system comets, linking outgassing rates in the comae and nuclei abundances depends not only on a comet's distance from the star (as clearly seen, e.g., for Hale–Bopp; Biver et al. 2002), but also on its detailed thermal and structural properties (Marboeuf et al. 2012). However, CO/H₂O and CO₂/H₂O outgassing rate ratios are generally measured near a comet's perihelion, where the process is dominated by water sublimation, which carries along CO and CO₂ trapped within it. Then, for a comet with a mostly homogeneous composition and dominated by clathrates (as shown to hold for 67P; Luspay-Kuti et al. 2016), outgassing ratios are expected to be a good representation of its ice abundance (see, e.g., Marboeuf & Schmitt 2014).

As with solar system comets, our exocometary abundance derivation also assumes that CO and CO₂ ice, as well as refractories, are released at rates that reflect their composition; in the context of our steady-state model, this is valid as long as two conditions apply:

• 1.  
  No CO or CO₂ is removed as ice by way of the smallest blow-out grains in the cascade, which is of order ∼10 μm in the Fomalhaut belt (though it can vary depending on grain composition). In other words, these ices are only removed as gas. This applies if ices are completely lost further up the collisional cascade, that is, if either the sublimation or photodesorption timescale is shorter than the collision timescale for grains larger than the blow-out size. Taking as an example a pure water ice grain 20 µm in size, the results of Grigorieva et al. (2007) can be used to show that its photodesorption timescale is ∼13,700 years. On the other hand, its collisional lifetime is ∼4 × 10⁵ years (Wyatt & Dent 2002), meaning that water ice cannot survive on the surface of grains at the bottom of the collisional cascade. A similar argument applies for a pure CO₂ ice grain, whose photodesorption timescale will be very similar to that of water, due to a similar photodesorption yield of ∼10⁻³ molecules photon⁻¹ (Grigorieva et al. 2007; Öberg et al. 2009).The next question is whether CO gas or CO₂ ice can remain trapped within refractories all the way down to the smallest sizes in the cascade. Due to its low sublimation temperature (∼20 K), CO is trapped already in gas form at the blackbody temperature of a dust grain within the Fomalhaut belt (∼48 K), and for small grains at the bottom of the cascade, it is likely to diffuse through the refractory layer. On the other hand, CO₂ has a higher sublimation temperature of ∼80 K (Collings et al. 2004). The hottest grains will be those near the blow-out limit (∼10 μm), which will have a temperature close to that observed in the spectral energy distribution (SED, ∼74 K; Kennedy & Wyatt 2014). Further heating is likely to take place through collisions themselves (e.g., with accelerated high-β grains, as proposed by Czechowski & Mann 2007), increasing the likelihood of CO₂ sublimation and subsequent release through diffusion. In conclusion, while further detailed modeling is necessary, we deem it very unlikely for any CO to be retained within grains down to the smallest sizes in the cascade, but on the other hand we cannot exclude that a fraction of CO₂ ice may instead be retained within a refractory mantle. The latter would cause an underestimate of the CO+CO₂ exocometary mass fractions presented here.
• 2.  
  CO and CO₂ production is not dominated by resurfacing collisions, which preferentially occur for large bodies at the top of the cascade (A. Bonsor et al. 2017, in preparation). These less energetic collisions, which are not taken into account in our model, can expose trapped volatiles as well as fresh surface ice that can then be rapidly lost through sublimation or photodesorption. This would produce more gas mass than predicted through catastrophic-only planetesimal collisions, meaning that our model would be overestimating the cometary CO+CO₂ ice mass fraction. If resurfacing collisions are the main driver for gas release, we expect most of this release to happen early in the lifetime of the belt, since it requires the big planetesimals (holding most of the CO mass) to only have been resurfaced rather than destroyed. However, we note that these large bodies may also be large enough to retain the released gas through an atmosphere, or they could lose dust as well as gas through drag (as observed in solar system comets). Overall, it remains to be established whether gas released through resurfacing collisions can dominate the released mass, though such a study is beyond the scope of this paper.

We find that CO+CO₂ mass fractions in exocomets are similar to each other and are of the same order (within about an order of magnitude) as observed for solar system comets (Le Roy et al. 2015 and references therein). As shown in Figure 6, this similarity appears to apply around host stars of a range of ages and spectral types. In terms of distance to the host star, exocometary gas is observed at 50–220 au around an 8.7 L_⊙ star (β Pic, Matrà et al. 2017), ∼81 au around a 3.3 L_⊙ star (HD 181327, Marino et al. 2016), and ∼135 au around a 16.6 L_⊙ star (this work, stellar luminosity from SED fit in Kennedy & Wyatt 2014). Assuming blackbody-like bodies, the equilibrium temperature in these belts would be equivalent to distances of 17–75, 45, and 33 au from the Sun in the solar system, meaning that the temperature of these exocomets should be of the same order as observed for the solar system's comet reservoir in the Kuiper Belt (30–50 au; e.g., Stern & Colwell 1997). Then, similar compositions and temperature environments for comets may indicate similar comet-formation conditions in younger protoplanetary disks, including the solar nebula.

Another aspect to consider is whether these cometary fractions of mass locked in CO and CO₂ are globally representative of the chemical heritage from the ISM (see Pontoppidan et al. 2014 for an extensive discussion of such inheritance). In the ISM, we know that the [CO/H₂] abundance ratio in the gas is of order 10⁻⁴, where H₂ dominates the gas mass, and the gas/dust ratio is of order ∼100. This yields a ${M}_{\mathrm{CO}\mathrm{gas},\mathrm{ISM}}/{M}_{\mathrm{total},\mathrm{solids},\mathrm{ISM}}\sim 14 \%$. On the other hand, the CO abundance in ISM ices is in the range [CO/H₂O]_ice,ISM ∼ 9%–36% compared to H₂O (Mumma & Charnley 2011, and references therein). The CO₂ content is dominated by its ice phase (e.g., van Dishoeck et al. 1996), with a [CO₂/H₂O]_ice,ISM abundance of ∼15%–44% (again, see Mumma & Charnley 2011 and references therein). This means that ISM material that is accreted in the outer regions of the protoplanetary disk will contain CO+CO₂ already in the ice form, with ${M}_{\mathrm{CO}+{\mathrm{CO}}_{2}\mathrm{ice},\mathrm{ISM}}/{M}_{{{\rm{H}}}_{2}{\rm{O}}\mathrm{ice},\mathrm{ISM}}\sim 0.51\mbox{--}1.64$, and CO in the gas form, with ${M}_{\mathrm{CO}\mathrm{gas},\mathrm{ISM}}/{M}_{\mathrm{total},\mathrm{solids},\mathrm{ISM}}\sim 0.14$.

Given their present location, we assume that cometary belts formed in protoplanetary disks outside the H₂O and CO₂ ice lines. This means that without significant vertical and outward radial mixing, H₂O and CO₂ remained locked on the grains with ISM abundances, producing cometary CO₂/H₂O abundances representative of the ISM. For the CO content, however, we consider two possible scenarios, where comets formed either outside or inside the CO ice line. Outside the CO ice line, the freeze-out of CO from ISM gas will enhance the ice-phase CO+CO₂ abundance compared to the ISM value. This simple scenario would yield a total CO+CO₂ mass fraction in the comets given by

$$\begin{eqnarray}&&{\left(\displaystyle \frac{{M}_{\mathrm{CO}+{\mathrm{CO}}_{2},\mathrm{ice}}}{{M}_{\mathrm{total},\mathrm{solids}}}\right)}_{\mathrm{comet}}=\displaystyle \frac{{({M}_{\mathrm{CO}+{\mathrm{CO}}_{2},\mathrm{ice}})}_{\mathrm{ISM}}+{({M}_{\mathrm{CO},\mathrm{gas}})}_{\mathrm{ISM}}}{{({M}_{\mathrm{total},\mathrm{solids}})}_{\mathrm{comet}}}\end{eqnarray} \tag{ 3 }$$

where

$$\begin{eqnarray}&&\displaystyle \frac{{({M}_{\mathrm{CO}+{\mathrm{CO}}_{2},\mathrm{ice}})}_{\mathrm{ISM}}}{{({M}_{\mathrm{total},\mathrm{solids}})}_{\mathrm{comet}}}=\displaystyle \frac{1}{1+{\left(\tfrac{{M}_{\mathrm{dust}}}{{M}_{\mathrm{ice}}}\right)}_{\mathrm{comet}}}\displaystyle \frac{1}{1+{\left(\tfrac{{M}_{\mathrm{CO}+{\mathrm{CO}}_{2},\mathrm{ice}}}{{M}_{{{\rm{H}}}_{2{\rm{O}}},\mathrm{ice}}}\right)}_{\mathrm{ISM}}^{-1}}\end{eqnarray} \tag{ 4 }$$

and

$$\begin{eqnarray}&&\displaystyle \frac{{({M}_{\mathrm{CO},\mathrm{gas}})}_{\mathrm{ISM}}}{{({M}_{\mathrm{total},\mathrm{solids}})}_{\mathrm{comet}}}=\displaystyle \frac{1}{1+{\left(\tfrac{{M}_{\mathrm{CO},\mathrm{gas}}}{{M}_{\mathrm{total},\mathrm{solids}}}\right)}_{\mathrm{ISM}}^{-1}}.\end{eqnarray} \tag{ 5 }$$

Given CO gas and CO+CO₂ ice abundances in the ISM, assuming a range of dust-to-ice ratios between the $\pm 1\sigma$ values measured in comet 67P, we obtain an expected CO+CO₂ cometary mass fraction of ∼17%–33% (darker shaded region in Figure 6). The growth of grains to comet-sized bodies allows CO ice originally on grain mantles to become trapped in other ices and refractories; this ensures that CO (and thus the CO+CO₂ mass fraction derived) can be retained once the protoplanetary disk is dispersed and the CO ice line moves farther out than the cometary belt location.

In the second formation scenario, where the cometary belt formed inside the CO ice line in the protoplanetary disk, no extra CO from ISM-like gas would be incorporated in the ice phase; this would imply a CO+CO₂ cometary mass fraction of ∼5%–21% (lighter shaded region in Figure 6). As well as assuming that all of the CO can remain trapped within other ices or refractories inside its ice line, this simplified evolutionary model ignores any grain surface chemistry, which is likely ongoing through the ISM, protoplanetary, and cometary phases of evolution. We expect that such chemistry will act to deplete the CO and CO₂ ice abundances, creating not only more complex volatiles (which are not dominant in cometary ice), but also organic refractories. Strictly speaking, we should therefore consider these ranges as upper limits to the CO+CO₂ cometary mass fractions derived in this ISM inheritance model.

While the large error bars in the solar system and exocometary measurements do not allow us to draw significant conclusions with regard to enhancement or depletion compared to ISM-inherited abundances, we show that such a comparison should be possible with increasingly accurate observations. For example, measuring the depletion of CO+CO₂ abundance with respect to ISM-inherited values would allow us to estimate the amount of CO or CO₂ that has been lost either due to sublimation during or immediately after dispersal of the protoplanetary disk (due to consequent outward movement of the CO ice line location) or due to grain surface chemistry and production of more complex organics. As well as achieving more accurate measurements over a larger sample of exo- and solar system comets, pinpointing the formation location of cometary belts within the protoplanetary disk with respect to the CO ice line will be crucial in allowing us to distinguish between the two possible ranges of ISM-inherited abundances discussed above.

In Section 3.5 we presented tentative evidence of an enhancement in CO emission near the planetesimal belt's pericenter. Here, we examine its possible physical origin in the framework of our steady-state model described in Section 4.3. Once again, we assume that the CO gas produced is a fraction of the solid mass lost as part of a steady-state collisional cascade, giving us access to the CO+CO₂ ice abundance. The total solid mass loss rate $\dot{M}(D)$ is calculated by multiplying the mass ${M}_{\mathrm{solid}}(D)$ in any given size bin with the collision rate ${R}_{\mathrm{col}}(D)$ of bodies in that size bin.

In Appendix A, we quantify how this mass loss rate within an eccentric cometary belt depends on the true anomaly f, for two possible regimes of grain sizes releasing CO. For the small grains, under the condition $D\lt {D}_{\min }($ 0.5 ${v}_{\mathrm{rel}}^{2}/{Q}_{{\rm{D}}}^{\star }{)}^{1/3}$, we find that the collision rate is independent of the true anomaly; this in turn means that the mass loss rate should be enhanced at the apocenter due to the expected enhancement in solid mass at this location (due to particles spending more time there). For larger grains ($D\gt {D}_{\min }($ 0.5 ${v}_{\mathrm{rel}}^{2}/{Q}_{{\rm{D}}}^{\star }{)}^{1/3}$), the azimuthal dependence of the mass loss rate is set by six parameters, namely the slope of the size distribution in the collisional cascade q, the mean proper eccentricity of material orbiting within the belt ${e}_{{\rm{p}}}$, the stellar mass ${M}_{\star }$, the belt semimajor axis a, the forced eccentricity of the belt ${e}_{\mathrm{frc}}$, and the specific incident energy required for a catastrophic collision, ${Q}_{{\rm{D}}}^{\star }$. For the case of the Fomalhaut belt, we have observational constraints on $q\sim 1.83$ (Ricci et al. 2012, where we assume it to be independent of true anomaly), ${M}_{\star }\sim 1.92$ M_⊙ (Mamajek 2012), $a\sim 136.3\,\mathrm{au}$, and ${e}_{\mathrm{frc}}\sim 0.12$ (MacGregor et al. 2017), meaning that our only free parameters are the mean proper eccentricity ${e}_{{\rm{p}}}$ and the specific strength ${Q}_{{\rm{D}}}^{\star }$ of the planetesimals. We can then quantify the fractional difference in mass loss rate at pericenter with respect to apocenter ($1-\dot{M}(f=180^\circ )/\dot{M}(f=0^\circ )$) using Equations (12) and (13) in Appendix A, shown for a wide range of ${e}_{{\rm{p}}}$ and ${Q}_{{\rm{D}}}^{\star }$ values in Figure 7.

Zoom In Zoom Out Reset image size

As expected, for low proper eccentricities and high planetesimal strengths, the effect of a higher mass at apocenter dominates over the effect of easier collisional disruption at pericenter. Vice versa, for high proper eccentricities and low planetesimal strengths, we see a pericenter enhancement in the mass loss rate due to the effect of easier collisional disruption at pericenter dominating over the effect of having a higher mass at apocenter. The upper and lower horizontal asymptotes of the curves correspond to the limits ${e}_{{\rm{p}}}\to 0$ and ${e}_{{\rm{p}}}\to \infty$ as calculated by setting ${e}_{\mathrm{frc}}=0.12$ in Equations (14) and (15) (Appendix A) and are independent of both ${e}_{{\rm{p}}}$ and ${Q}_{{\rm{D}}}^{\star }$. This implies a ∼15% maximum pericenter-to-apocenter enhancement in the mass loss rate of the Fomalhaut belt. For the $\pm 1\sigma$ range of proper eccentricities derived from continuum imaging (MacGregor et al. 2017) and ${Q}_{{\rm{D}}}^{\star }$ values in the range 10–10,000 J kg⁻¹, expected for weak ice grains of a wide range of sizes within the collisional cascade (Wyatt & Dent 2002; assuming the results of Benz & Asphaug 1999), we predict the fractional difference in mass loss rate between pericenter and apocenter in the range −9% to 12%. Then, given that the CO photodissociation rate is independent of true anomaly, and assuming that the CO+CO₂ ice fraction is too, the steady-state mass loss rate enhancement at pericenter directly translates into a CO mass enhancement.

Whether a CO mass enhancement translates into a CO flux enhancement at pericenter/apocenter, or in other words CO pericenter/apocenter glow, also depends on the excitation of the CO molecule (see Section 3.4), and particularly on the radial dependence of the gas kinetic temperature and electron density. This is because for example in Fomalhaut, CO emitted at the ring's pericenter will be closer to the star than the CO at apocenter, by a factor of $(1+e)/(1-e)\sim 1.27$. Making the simple assumption of a β Pic-like environment for CO excitation, with an electron density varying with radius as $\sim 300{(R/100\mathrm{au})}^{-1}$ cm⁻³ (Matrà et al. 2017), the expected flux at pericenter would be 2.5%–6.1% larger than that at apocenter, making the CO flux enhancement even larger than expected from a mass enhancement only. Therefore, we expect excitation effects to favor a CO flux enhancement at pericenter. We cautiously note that such a contribution from CO excitation will depend on the gaseous environment of the Fomalhaut belt; this may significantly differ from that of β Pic, though it may be characterized in the future using optically thin line ratio observations.

Finally, we need to take into account projection effects due to the viewing geometry; in particular, the observed pericenter-to-apocenter ratio will depend on the combined effect of the ring's vertical thickness, its inclination to the line of sight, and the on-sky angular distance between the ring's pericenter and the nearest ansa. This is because the azimuthal intensity distribution for a significantly inclined ring with a nonnegligible vertical thickness should present two enhancements at the location of the two ansae (see Marino et al. 2016 for the dependence of the azimuthal intensity distribution on a ring's scale height). To fully take this effect into account and demonstrate its impact on the measured pericenter or apocenter glow, we produce a simple model of CO in the Fomalhaut ring and produce a sky-projected model image using the RADMC-3D ¹⁵ radiative transfer code. In doing this, we first construct an eccentric CO ring model with radial mass distribution, inclination to the line of sight, position angle, forced eccentricity, and argument of pericenter equal to those of the dust ring. Then, we introduce the dependence of the CO surface density on the true anomaly as described above and derived in Appendix A. We assume the maximum possible pericenter versus apocenter CO mass enhancement for the best-fit ${e}_{{\rm{p}}}\sim 0.06$ from the continuum observations, corresponding to ∼12% (for a ${Q}_{{\rm{D}}}^{\star }$ of 10 J kg⁻¹). Furthermore, we make the simple assumptions that the gas kinetic temperature is equal to the blackbody temperature of the planetesimals and that the electron density follows the same radial dependence as found in β Pictoris.

In our model (see the spectrally integrated image in Figure 8), we assume a radially constant CO vertical aspect ratio equal to the average proper eccentricity of the planetesimals ($h\sim 0.06$), but we also consider extreme cases of a very small aspect ratio ($h\sim 0.0001$) and a larger one ($h\sim 0.14$, corresponding to the 3σ upper limit found for the HD 181327 ring; Marino et al. 2016). We find that in the infinitesimally vertically thin limit we recover the pericenter enhancement with respect to apocenter expected from our CO steady-state model from both mass (12%) and excitation (∼2.5%) effects. In the vertically thick cases, the value of this ratio is slightly reduced due to the additional "background" effect of the ansae enhancement. Their absolute flux difference, however, remains unchanged, due to the projection effect being axisymmetric with respect to the ring's geometric center.

Zoom In Zoom Out Reset image size

To conclude, the upper-limit pericenter versus apocenter enhancement predicted by our model (14.5%) is marginally consistent with our measurement from the current data set (88 ± 25%; see Section 3.5) at the 2.9σ level. If confirmed at high significance, this discrepancy would indicate the inability of our steady-state model to explain this asymmetry. In turn, this may favor a stochastic event such as the destruction of a large icy body near the belt's pericenter (as discussed in Section 4.2). Similarly, a recent impact was previously invoked to explain the observed CO asymmetry in the β Pictoris system (Dent et al. 2014), though this was recently ruled out through higher-resolution observations (Matrà et al. 2017).

Either way, a conclusive detection of this enhancement would allow us to determine the direction of orbital motion of the belt, and in turn whether its east or west side is nearest to Earth. This is because the finite lifetime of CO (∼120 years in such an optically thin environment) would cause either a tail in the direction of motion (in case of a stochastic event; see, e.g., Dent et al. 2014) or a ∼19° offset of the peak location (with respect to pericenter/apocenter, also in the direction of motion) due to CO pericenter/apocenter glow. Therefore, future, deeper ALMA observations of pericenter and apocenter are warranted to confirm this tentative evidence for asymmetry.