New Millimeter CO Observations of the Gas-rich Debris Disks 49 Cet and HD 32297

Figures 1(a) and (e) display the millimeter continuum images of our targets. Both sources are clearly detected, at a peak signal-to-noise ratio (S/N) of 14 and 31 for 49 Cet and HD 32297, respectively. The map of 49 Cet (Figure 1(a)) suggests that the disk is marginally resolved. To constrain the spatial distribution of emitting dust particles and to determine their total flux density at 1.33 mm we applied the uvmodelfit task in CASA that allows one to fit a single component source model to the uv data. Channels with CO emission were flagged and the data weights were recomputed using the statwt task.

Zoom In Zoom Out Reset image size

For 49 Cet an elliptical Gaussian model was fitted to the non-flagged visibility data. The derived disk geometry is in good agreement with the published results of Hughes et al. (2017) which were based on a much higher spatial resolution ALMA data set. Using uvmodelfit we measure a flux density of 5.3 ± 0.7 mJy, where the quoted uncertainty was computed as the quadratic sum of the measurement error and the typical ALMA calibration error (10% of the signal). By fitting the previously available (sub)millimeter photometry longward of 400 μm (at 0.45, 0.5, 0.85, and 9 mm, taken from Moór et al. 2015b; MacGregor et al. 2016; Holland et al. 2017; Hughes et al. 2017) with a straight line in the log–log space we can infer a flux density of 4.7 mJy at 1.33 mm. Our measurement is consistent with this prediction.

Utilizing an elliptical Gaussian model we found the disk around HD 32297 to be unresolved. Therefore in this case we applied a point-source model in the visibility domain fitting, resulting in a total flux density of 3.4 ± 0.4 mJy. The obtained center of the fitted point source is consistent with the stellar position. Previous millimeter continuum observations of HD 32297 with IRAM/MAMBO and with Submillimeter Array (SMA) at 1.2 and 1.3 mm yielded, in good agreement with our result, flux densities of 3.14 ± 0.82 mJy and 3.10 ± 0.74 mJy (Meeus et al. 2012). By analyzing their recent 1.3 mm high angular resolution ALMA observations of HD 32297, MacGregor et al. (2018) derived a total (as a combination of a planetesimal belt and a halo component) flux density of 3.7 ± 0.3 mJy. This result matches well the flux derived from our low-resolution observation, which indicates that no extended emission was filtered out by ALMA in the earlier much higher-resolution measurement.

Assuming isothermal, optically thin dust emission and using our obtained millimeter flux densities, we estimated the dust mass of the disks as (Hildebrand 1983)

$$\begin{eqnarray}&&{M}_{d}=\displaystyle \frac{{F}_{\nu }{d}^{2}}{{B}_{\nu }({T}_{\mathrm{dust}}){\kappa }_{\nu }},\end{eqnarray} \tag{ 1 }$$

where d is the distance, B_ν (T_dust) is the Planck function at 1.33 mm for a characteristic dust temperature ${T}_{\mathrm{dust}}$, and ${\kappa }_{\nu }$ is the dust opacity at the given frequency. Analysis of the SEDs implied that both disks have two components (Donaldson et al. 2013; Holland et al. 2017). The millimeter emission predominantly arises from the colder outer belts (Hughes et al. 2017; MacGregor et al. 2018). For HD 32297, MacGregor et al. (2018) derived a characteristic temperature of 47 K for the emitting large grains in the outer belt. 49 Cet harbors a very broad outer disk whose surface density peaks around 100 au and decreases toward larger radii (Hughes et al. 2017). Therefore to estimate a characteristic temperature for large particles we took the equilibrium temperature of blackbody grains located at a radial distance of 100 au, and obtained T_dust = 56 K. Following Beckwith et al. (1990), we adopted a value of 2.3 cm² g⁻¹ for the dust opacity. As a result we obtained dust masses of 0.15 ± 0.02 M_⊕ and 0.64 ± 0.07 M_⊕ for 49 Cet and HD 32297, respectively. The quoted uncertainties consider the errors of the measured flux density, of the distance of the object, and of the characteristic dust temperatures, but do not account for the probable factor of ∼2–3 uncertainty associated with the dust opacity (Miyake & Nakagawa 1993; Ossenkopf & Henning 1994).

By inspecting our CO data cubes, in addition to detecting ¹²CO emission at high peak S/N (>12), we discovered ¹³CO emission at peak S/Ns of >9 around the position and radial velocity of both stars. For HD 32297 even the C¹⁸O line was clearly detected. We started the analysis with the ¹²CO and ¹³CO data. We applied different-sized elliptical apertures with an axis ratio and a position angle identical to the appropriate beam to extract the line spectra from the data cubes and to determine the integrated line flux by summing all consecutive channels that show significant line emission, i.e., where the peak S/N is higher than 3 (within the aperture). By examining the obtained integrated line fluxes as a function of the applied apertures we then defined the minimum aperture size that includes all CO emission associated with the disk. We obtained identical aperture sizes for ¹²CO and ¹³CO data for both disks. In the case of C¹⁸O data the integration was performed with the same aperture radius and over the same velocity range as established for the more abundant species. To estimate the statistical error of the derived line fluxes we used the line-free regions of the spectra, and the final uncertainty was computed by adding a 10% flux calibration error to this. Zeroth moment maps and the acquired spectra for the examined transitions are displayed in Figure 1, and the obtained integrated line fluxes with their uncertainties are listed in Table 1. In the case of 49 Cet no significant C¹⁸O emission was detected at the position of the star but a source appears close to it with an offset of 57 (the positional uncertainty, estimated as the ratio of the beam size to achieved S/N, is ∼17). By centering the aperture on this source we obtained a spectrum whose shape resembles well that of the ¹²CO and ¹³CO lines of 49 Cet. Integration over the above specified velocity range results in an integrated line flux of 257 ± 58 mJy km s⁻¹. The similarity of the spectrum and the fact that there is no counterpart of this source in the ¹²CO and ¹³CO maps suggest that the observed emission comes from the disk. However, since we have no explanation for the observed large offset, the integrated line flux is considered as an upper limit in the following analysis.

As Figures 1(i) and (j) demonstrate, the obtained CO spectra show typical double-peaked profiles indicating that the observed emission arises from rotating gas. Using the ALMA 12 m Array Hughes et al. (2017) obtained a high-S/N ¹²CO (3–2) line spectrum for 49 Cet (see their Figure 3). The shape and width of our ¹²CO (2–1) line are consistent with those of that 3–2 spectrum. Previous observation of the 2–1 line with SMA resulted in a line flux of 2.0 ± 0.3 Jy km s⁻¹ (Hughes et al. 2008), which is about two times lower than our value.

The first detection of CO emission toward HD 32297 was made by the 15 m James Clerk Maxwell Telescope (JCMT) radio telescope (Greaves et al. 2016). The obtained ¹²CO (2–1) line profile (Figure 3 in Greaves et al. 2016) is ∼18 km s⁻¹ wide. The line in our observation is only about half as broad and mostly corresponds to the blue wing of the JCMT spectrum in terms of its velocity range. Greaves et al. (2016) used two different baseline-fitting models and derived peak fluxes of ∼320 mJy and ∼440 mJy, i.e., >1.8 times higher than in our case. As a consequence of the mentioned differences their derived integrated line fluxes (2.7 Jy km s⁻¹ or 5.1 Jy km s⁻¹, depending on the applied baseline model) are also significantly higher than ours. Recently MacGregor et al. (2018) presented a ¹²CO (2–1) observation for HD 32297 with ALMA. The velocity interval of the CO emission (as shown in their position–velocity diagram, their Figure 3, left) as well as the derived integrated line flux of 1.02 Jy km s⁻¹ are consistent with our results. Therefore the somewhat deviating line observation of Greaves et al. (2016), derived from a lower-S/N measurement, will not be taken into account in the further analysis.

Using the obtained line fluxes, we derived line flux ratios of ${S}_{{}^{12}\mathrm{CO}}/{S}_{{}^{13}\mathrm{CO}}=2.3\pm 0.2$ and ${S}_{{}^{13}\mathrm{CO}}/{S}_{{{\rm{C}}}^{18}{\rm{O}}}\gt 6.5$ for 49 Cet and ${S}_{{}^{12}\mathrm{CO}}/{S}_{{}^{13}\mathrm{CO}}=2.1\pm 0.2$ and ${S}_{{}^{13}\mathrm{CO}}/{S}_{{{\rm{C}}}^{18}{\rm{O}}}=1.8\pm 0.3$ for HD 32297. Taking isotope ratios of [¹²C]/[¹³C] = 77 and [¹⁶O]/[¹⁸O] = 560, typical in the local interstellar matter (Wilson & Rood 1994), and assuming that all three isotopologues are optically thin and have identical excitation temperatures between 10 and 100 K, one would expect line ratios of ${S}_{{}^{12}\mathrm{CO}}/{S}_{{}^{13}\mathrm{CO}}\sim 85-91$ and ${S}_{{}^{13}\mathrm{CO}}/{S}_{{}^{18}\mathrm{CO}}\sim 7.3$. Thus the measured values imply highly optically thick ¹²CO emission in both disks, while the ¹³CO line is probably optically thick in HD 32297 and optically thin in 49 Cet.

We used the optically thin line observations to estimate the CO gas mass. As a first step we derived the mass of the specific ^yC^zO isotopologue as:

$$\begin{eqnarray}&&{M}_{{}^{y}{C}^{z}O}=4\pi {{md}}^{2}\displaystyle \frac{{S}_{21}}{{x}_{2}h{\nu }_{21}{A}_{21}},\end{eqnarray} \tag{ 2 }$$

where m is the mass of the given molecule, d is the distance of the target, and h is the Planck constant. S₂₁ is the integrated line flux, A₂₁ and ν₂₁ are the Einstein coefficient and the rest frequency of the specific rotational transition, while x₂ is the fractional population of the upper level.

We assumed that local thermodynamical equilibrium (LTE) holds and utilized the Boltzmann equation to compute the fractional populations. Our observations do not allow us to derive reliable gas temperatures for the targeted disks. Gas temperature estimates in CO-rich debris disks are uncertain, and previous studies typically resulted in low values. Based on our CO observations, we derived ≲10 K for the disk of HD 21997 (Kóspál et al. 2013). Using the measured peak flux value of the high angular resolution ¹²CO (3–2) map of 49 Cet (Hughes et al. 2017) we could derive a maximum brightness temperature of ∼26 K for that disk. For HD 141569, a system that resembles the most CO-rich debris disks in terms of dust/CO mass (Section 5.1), Flaherty et al. (2016) found a gas temperature of ${27}_{-4}^{+11}$ K at a radial distance of 150 au by modeling CO observations. In all these cases the derived temperatures are lower than the characteristic dust temperatures of the given systems.

Considering these results, for the excitation temperature we adopted 20 K for both our disks. Based on Equation (2) we then obtained ${M}_{{}^{13}{{\rm{C}}}^{16}{\rm{O}}}$ = (1.5 ± 0.2) × 10⁻⁴ M_⊕ for 49 Cet and ${M}_{{}^{12}{{\rm{C}}}^{18}{\rm{O}}}$ = (1.4 ± 0.3) × 10⁻⁴ M_⊕ for HD 32297. The uncertainties account only for errors in the integrated line fluxes and distances of the disks. Different gas temperatures and/or non-LTE (NLTE) conditions would result in different x₂ values and thereby gas masses. Since with the adopted temperatures the fractional population is close to its maximum in LTE, these deviations can typically lead to higher gas masses (lower x₂). We note that if LTE holds, for gas temperatures between 7 and 74 K the gas mass estimates would not exceed our derived CO masses by more than a factor of two.

To estimate the total ¹²CO gas quantity, the obtained CO isotopologue masses were multiplied by the corresponding isotope ratios and took into account the difference in the molecular masses. These yielded CO gas masses ${M}_{{}^{12}{{\rm{C}}}^{16}{\rm{O}}}$ = (1.11 ± 0.13) × 10⁻² M_⊕ and (7.4 ± 1.3) × 10⁻² M_⊕ for 49 Cet and HD 32297, respectively. Different physical/chemical mechanisms, perhaps most importantly the isotope-selective photodissociation, can alter the abundance of the isotopologues, resulting in isotope ratios substantially different from the ones typical in the local interstellar medium (e.g., Visser et al. 2009). This can introduce serious uncertainties in the estimates of the total CO gas mass (see Section 5.2 and e.g., Miotello et al. 2014).

Our observations of rarer CO isotopologues in debris disks around 49 Cet and HD 32297 implied that their ¹²CO emission is optically thick and their total CO mass is higher than 0.01 M_⊕. These two objects do not stand alone with this characteristic: previous measurements revealed similarly CO-rich debris disks around young A-type stars HD 21997, HD 121617, HD 131488, and HD 131835. The estimated CO content of these systems is on average three orders of magnitude higher than that of the gaseous debris disk of β Pic (∼3.6 × 10⁻⁵ M_⊕, Matrà et al. 2018b) and is more similar to those of protoplanetary disks. To explore the latter point further, in Figure 2 we plotted the CO masses of these six debris disks (circles) along with disks around 10 nearby (≲160 pc) Herbig Ae stars (squares) as a function of age. Dust masses of the disks are also shown by colors of the symbols. These Herbig Ae disks can be considered as the precursors of young debris disks with A-type host stars. While their average dust mass is about two orders of magnitude higher than that of our debris disks, in terms of their CO gas quantity the difference is only one order of magnitude on average and there are some less massive Herbig Ae disks (e.g., V856 Sco) with CO content comparable to that of the debris sample. The disk around the well-known pre-main-sequence star, HD 141569, resembles CO-rich debris disks both in terms of dust and CO masses.

Zoom In Zoom Out Reset image size

Do the protoplanetary-like CO content and the high CO to dust mass ratios inevitably mean that 49 Cet and HD 32297 (as well as the other four debris disks) harbor primordial gas? An ultimate way to answer this question would be to measure the H₂ content of these disks. If the gas, similar to protoplanetary disks, is predominantly primordial then hydrogen molecules are the primary gaseous species. On the other hand, if the gas has a secondary origin then the abundance of H₂ molecules in the gas mixture is expected to be low (Kral et al. 2017, 2019), and the gas is rather composed of different molecules released from the icy bodies and their photodissociation products. This also means that the observed comparable CO masses do not neccessarily indicate similar total gas masses. Unfortunately, the detection of molecular hydrogen is notoriously difficult. Therefore here we use a different strategy and investigate whether the observed CO gas quantity of 49 Cet and HD 32297 could be explained within the framework of a secondary gas production scenario (see Kral et al. 2019).

High angular resolution millimeter continuum observations of 49 Cet and HD 32297 (Hughes et al. 2017; MacGregor et al. 2018) implied that they—similar to the other four CO-debris disks mentioned above—harbor massive, cold (<110 K) dust belts at radial distances larger than ∼30 au with which the observed CO gas is at least partly colocated. Supposing that the parent planetesimals are icy—which is permitted by the typical temperatures in the disks of 49 Cet and HD 32297—mutual collisions between these bodies produce not only smaller and smaller particles but can release gas as well (Moór et al. 2011; Zuckerman & Song 2012; Kral et al. 2017). Based on volatiles detected in spectroscopic surveys of solar system comets (Mumma & Charnley 2011), water, carbon monoxide, and carbon dioxide may be the most important products of this process. Due to stellar and interstellar UV photons, the released gas molecules are photodissociated. According to Kral et al. (2017), assuming a steady-state balance of gas production and loss, the CO mass of a disk can be estimated as

$$\begin{eqnarray}&&{M}_{\mathrm{CO}}={\dot{M}}_{\mathrm{CO}}{t}_{\mathrm{ph}}{\epsilon }_{\mathrm{CO}}={\dot{M}}_{\mathrm{loss}}\gamma {t}_{\mathrm{ph}}{\epsilon }_{\mathrm{CO}},\end{eqnarray} \tag{ 3 }$$

where ${\dot{M}}_{\mathrm{loss}}$ is the mass-loss rate in the collisional cascade, γ is the CO+CO₂ ice mass fraction of planetesimals, t_ph is the photodissociation timescale of unshielded CO gas, while _CO is a factor that shows how well shielded the molecules are. Here we consider that photodissociation of CO₂ molecules also contribute to the observed CO gas. If we know t_ph, γ, and ${\dot{M}}_{\mathrm{loss}}$ then the minimum necessary shielding factor in the 49 Cet and HD 32297 systems can be estimated. As a lower threshold for the UV flux, we consider only the influence of interstellar UV photons and following Visser et al. (2009) we set a dissociation timescale of ∼120 yr for unshielded CO molecules. In solar system comets the ice mass fraction is not higher than 0.27 (Mumma & Charnley 2011; Matrà et al. 2017b), thus we adopt this value for the γ parameter. Assuming a steady-state collisional cascade in the planetesimal belt, we estimated the mass-loss rate using the formula proposed by Matrà et al. (2017a, see their Equation (21)):

$$\begin{eqnarray}&&{\dot{M}}_{\mathrm{loss}}({M}_{\oplus }\,{\mathrm{Myr}}^{-1})=1200{\left(\displaystyle \frac{R}{\mathrm{au}}\right)}^{1.5}{\left(\displaystyle \frac{{\rm{\Delta }}R}{\mathrm{au}}\right)}^{-1}{f}_{d}^{2}\left(\displaystyle \frac{{L}_{* }}{{L}_{\odot }}\right){\left(\displaystyle \frac{{M}_{* }}{{M}_{\odot }}\right)}^{-0.5}.\end{eqnarray} \tag{ 4 }$$

Based on MacGregor et al. (2018), for the outer belt of HD 32297 we adopted a radius (R) and width (Δ R) of 100 au and 40 au, respectively. In the case of the very broad disk of 49 Cet we used the double power-law surface density profile model proposed by Hughes et al. (2017) and the inner (71 au) and outer radii (153 au) of the disk was set at half of the surface density maximum. This resulted in R = 112 au and Δ R = 82 au. Fractional dust luminosities and stellar parameters needed for this calculation were taken from Section 2 and from the literature (Hughes et al. 2017; MacGregor et al. 2018). Using Equation (4) then we obtained mass-loss rates of 0.1 and 5.2 M_⊕ Myr⁻¹ for 49 Cet and HD 32297, respectively. Applying these values in Equation (3) yielded shielding factors of _CO ∼ 3600 and ∼400 for 49 Cet and HD 32297, respectively. Since the ice mass fraction is likely lower than 0.27 and—especially in the case of 49 Cet—stellar UV photons can affect the dissociation of CO molecules in the inner part of the disk these factors should be considered as lower limits. We note that the lower _CO value of HD 32297 is due to its prominently high fractional luminosity (the dust production rate being proportional to ${f}_{d}^{2}$).

What processes can contribute to the shielding in these systems? While in protoplanetary disks continuum shielding by dust can substantially reduce the flux of incoming UV photons, in debris disks their absorption by optically thin dust is insignificant. Since photodissociation of CO molecules occurs primarily through discrete absorptions into bound electronic states, they are subject to self-shielding. According to Visser et al. (2009) at ¹²CO column densities higher than 10¹⁵ cm⁻² self-shielding can reduce significantly the photodissociation rate of ¹²CO molecules located in deeper regions. By interpolating in their Table 5 we found that to achieve the above shielding factors the vertical ¹²CO column density in the disks around 49 Cet and HD 32297 must exceed ∼1.7 × 10¹⁹ cm⁻² and ∼1.5 × 10¹⁸ cm⁻², respectively. Though MacGregor et al. (2018) did not provide a detailed model for the distribution of CO gas in HD 32297 they argued that it is colocated with dust. Assuming that the gas follows the same spatial distribution as the dust component and taking the total CO mass derived in Section 4.2, we derived a maximum vertical column density of ∼10¹⁸ cm⁻² for this disk. This indicates that at least in the densest regions the self-shielding mechanism can play an important role in the long-term survival of ¹²CO molecules. In the case of the less CO-rich 49 Cet, the typical vertical CO column densities fall short of the value required for effective shielding. In Figure 3 we display the necessary shielding factors for our two targets together with the other four high CO mass debris disks and the disk of β Pic. Stellar and disk parameters needed for calculations of ${\dot{M}}_{\mathrm{loss}}$ were taken from the literature (Lieman-Sifry et al. 2016; Hughes et al. 2017; Moór et al. 2017; Kral et al. 2019; MacGregor et al. 2018; Matrà et al. 2018a). The obtained ${\dot{M}}_{\mathrm{loss}}$ values range between 0.04 and 4.2 M_⊕ Myr⁻¹. The CO mass estimates are from Figure 2, except for β Pic where it was taken from Matrà et al. (2018b). As the plot shows, while the CO content of β Pic can be explained with an unshielded or weakly shielded model, all CO-rich disks require substantial shielding with _CO > 100, i.e., following the proposal of Kral et al. (2019) the latter objects can be rightfully called shielded debris disks.

Zoom In Zoom Out Reset image size

CO column densities needed for efficient self-shielding of these systems are also shown in Figure 3. By comparing these values with the vertical CO column densities of the given disks we found that self-shielding can play an important role only in the HD 32297 and HD 121617 systems. Even in these two disks, however, more efficient shielding is needed for explaining the observed CO quantity. Note that all of our previous calculations on shielding dealt only with ¹²CO molecules whose total mass was estimated from optically thin C¹⁸O or ¹³CO observations assuming abundance ratios to be equal to the [¹⁶O]/[¹⁸O] and [¹²C]/[¹³C] isotope ratios measured in the local interstellar matter (see Section 4.2 and Moór et al. 2017). UV absorption lines of ¹³CO and C¹⁸O are only partly overlapping with those of ¹²CO making these lower abundance isotopologues more vulnerable to photodissociation; for example, based on Visser et al. (2009), in a disk where only self-shielding is taken into account the photodissociation rate of C¹⁸O molecules—depending on the ¹²CO column densities—could be up to ∼15 times higher than that of ¹²CO. By leading to higher ¹²CO/¹³CO and ¹²CO/C¹⁸O ratios this isotope-selective photodissociation would result in higher ¹²CO masses than we assumed in our analysis thereby further increasing the need for more efficient shielding mechanisms. This also means that to avoid the uncertainties associated with the abundance ratios and construct a more self-consistent model we rather have to reproduce the results of the optically thin ¹³CO or C¹⁸O line observations.

We note that the currently known sample of CO-bearing debris disks may contain three additional young shielded debris disks. HD 138813 and HD 156623 exhibit higher ¹²CO line luminosities than 49 Cet and HD 21997. Though, because of the lack of ¹³CO and/or C¹⁸O observations, their CO mass is less constrained, the lower mass limits derived from the ¹²CO measurements imply shielding factors of >10 and >100 for HD 156623 and HD 138813, respectively. These results suggest that these systems likely harbor CO-rich shielded debris disks too. A third system, HD 121191, harbors a CO mass of 2.5 × 10⁻³ M_⊕ (taking into account the new Gaia DR2 distance; Moór et al. 2017). This remains below the average CO mass of the abovementioned six systems but almost two orders of magnitude higher than that of β Pic. Though the paucity of information on the spatial distribution of gas and dust and their relationship in this system makes the modeling less reliable, this disk might also require strong shielding (Moór et al. 2017; Kral et al. 2019). Similar to the previously discussed CO-rich systems all these three disks surround A-type star with ages <50 Myr.

In the literature two models are proposed to explain the existence of shielded debris disks. The hybrid disk model (Kóspál et al. 2013; Péricaud et al. 2017) is motivated by the fact that all known representatives of this class are very young. According to this scenario, the gas in these systems is mostly primordial and dominated by hydrogen molecules that provide effective shielding for CO. The observed CO molecules could partly be primordial and collisions of icy bodies and grains—that are thought to be responsible for the continuous replenishment of the dust component—may also produce gas. Recently Kral et al. (2019) proposed an alternative model in which both the dust and gas components are second generation. According to this scenario second-generation molecules, including CO, are produced from volatile-rich solid bodies located in a planetesimal belt. Should the CO production rate be sufficiently high and the viscous evolution of the gas slow enough, neutral atomic carbon gas, the photodissociation product of CO, can accumulate so considerably that it becomes optically thick to UV radiation. Since the photoionization of C⁰ occurs at the same UV wavelengths as the photodissociation of the CO molecules, the carbon component could shield CO efficiently. This leads to the formation of a CO-rich shielded debris disk. Kral et al. (2019) found that the disks around HD 21997, HD 121191, HD 121617, HD 131488, and HD 131835 can be explained with this model.

5.3.1. Upgraded Shielded Secondary Gas Disk Model

5.3.2. Modelling Results

In the course of modeling we assumed that the collisional cascade and thereby the production of secondary CO gas is initiated 5 Myr after the birth of the systems and then our simulations are run for 40 Myr and 25 Myr for 49 Cet and HD 32297, respectively. For HD 32297, we adopted a planetesimal belt between 80 and 120 au, while in the case of 49 Cet we used a disk model with inner and outer radii of 71 au and 153 au, respectively (see Section 5.2). We took a uniform gas temperature of 20 K (Section 4.2) in both disks. In both cases we hypothesized that collisions between planetesimals result in fragmentation in the whole disk throughout our simulation, i.e., the disk is prestirred (Wyatt 2008). When releasing the gas, we initially assume a constant mass input rate with radial distance to not favor any gas-release mechanism over others. To find a best-fit model to the ¹³CO and C¹⁸O masses derived from observations for 49 Cet and HD 32297, we run a large grid of 100 models to explore the parameter space in ${\dot{M}}_{\mathrm{CO}}$ and α. The grid is logarithmic for both parameters and comprised of 10 elements going from 10⁻³ to 10⁻¹ M_⊕ Myr⁻¹ for ${\dot{M}}_{\mathrm{CO}}$ and from 10⁻⁵ to 10⁻² for α. For each model, we compute the total ¹³CO (for 49 Cet) or C¹⁸O (for HD 32297) masses and compare them to observed values. Masses found by the model for each part of the parameter space are shown in Figure 4. Solid black contours correspond to the observed ¹³CO and C¹⁸O masses, while the black dashed contours mark the 3σ uncertainties (Table 1).

Zoom In Zoom Out Reset image size

For 49 Cet (left panel), we found that models with ${\dot{M}}_{\mathrm{CO}}$ ranging from 0.004 to 0.01 M_⊕ Myr⁻¹ and with any α can explain the measured ¹³CO mass. To further constrain the acceptable parameter space we need to consider other observational properties. By combining ALMA CO (3–2) observations with previous SMA CO (2–1) data, Hughes et al. (2017) derived an inner-disk radius of 20 au as a best fit. Detecting C i ³P₀–³P₁ emission (at a rest frequency of 492.161 GHz) toward 49 Cet with the single-dish ASTE radio telescope, Higuchi et al. (2017) derived 25.8 ± 2.3 Jy km s⁻¹ for the integrated line flux.¹⁰ Models with different α and ${\dot{M}}_{\mathrm{CO}}$ result in quite different inner CO disk radii (${R}_{\mathrm{in},\mathrm{CO}}$) and C i emissions. The solid blue line in Figure 4 (left) corresponds to models where ${R}_{\mathrm{in},\mathrm{CO}}=20\,\mathrm{au}$, the best-fit value suggested by Hughes et al. (2017), while the solid green marks models which reproduce the measured C i integrated line flux. In the computation of C i model fluxes we used the LIME radiation transfer code (Brinch & Hogerheijde 2010) assuming LTE conditions. The dashed blue and green lines denote the models are consistent at a 3σ level with the observed parameters. We found that there exists a region in the parameter space, at 2 × 10⁻⁵ ≲ α ≲ 9 × 10⁻⁵ and $0.004\lesssim {\dot{M}}_{\mathrm{CO}}\lesssim 0.006$ M_⊕ Myr⁻¹, where all three observational constraints, represented by the blue (inner CO disk radius), green (C i emission), and black (¹³CO mass) lines, are simultaneously satisfied within their 3σ uncertainties (dashed lines). As a representative grid point for the subsequent simulation, we adopted ${\dot{M}}_{\mathrm{CO}}\sim 4.6\times {10}^{-3}$ M_⊕ Myr⁻¹ and α ∼ 4.6 × 10⁻⁵.

As Figure 4 (right) demonstrates, in the case of HD 32297 the observed C¹⁸O mass can be reproduced with CO production rates between 0.01 and 0.04 M_⊕ Myr⁻¹. Although MacGregor et al. (2018) presented no detailed model for the spatial distribution of ¹²CO, they argued that the gas and dust components in this disk are colocated. For the inner edge of the dust disk they inferred a radius of 78.5 ± 8.1 au. Figure 4 (right) shows models with 70 au and 55 au inner CO disk radii represented by solid and dashed blue lines, respectively. The latter value can be considered as a lower limit for the inner radius. These models suggest a constraint of α ≲ 2 × 10⁻⁵. Models with ${R}_{\mathrm{in},\mathrm{CO}}=70\,\mathrm{au}$ have no intersection with the black line in the explored parameter space. Furthermore we also displayed the loci of models (gray lines) that produce ¹²CO line fluxes corresponding to the observed value (1.05 Jy km s⁻¹, Table 1). In the calculations we assumed two different inclination values, 88° (dashed line) and 84° (solid), the best-fit parameters inferred from scattered light (Boccaletti et al. 2012) and continuum millimeter (MacGregor et al. 2018) observations of the disk. Model line-flux computations were performed with the LIME tool using an LTE approach. Intersection of the latter loci with the black lines is broadly consistent with the ${R}_{\mathrm{in},\mathrm{CO}}=55\,\mathrm{au}$ models. The grid point located closest to these crossings has ${\dot{M}}_{\mathrm{CO}}\sim 3.6\times {10}^{-2}\,{M}_{\oplus }\,{\mathrm{Myr}}^{-1}$ and α ∼ 1.5 × 10⁻⁵ for HD 32297.

We note that in the case of HD 32297 the agreement between the model and data is more marginal than for 49 Cet and requires very low α values. As we mentioned earlier, because of strong shielding, the ionization fraction of this disk could be significantly lower than in the more tenuous gas disk of β Pic. The lower ionization fraction will also favor ambipolar diffusion that could become important in these low gas density disks. Using Equation (4) in Kral & Latter (2016), we calculate that the ambipolar parameter for carbon is $\sim 100{\left(\tfrac{T}{20{\rm{K}}}\right)}^{-1/2}\left(\tfrac{{{\rm{\Sigma }}}_{n}}{{10}^{-5}\,{\rm{g}}\,{\mathrm{cm}}^{-2}}\right)\left(\tfrac{f}{{10}^{-2}}\right)$, where T is the gas temperature, ${{\rm{\Sigma }}}_{n}$ the neutral carbon density and f the ionization fraction of carbon. We see that indeed, in HD 32297, for an ionization fraction lower than about 10⁻² (which is typical for shielded disks as shown in Kral et al. 2019), the ambipolar parameter becomes lower than 100 and thus the ions and neutrals are not well coupled, which provides very low α values in MRI simulations (see Kral & Latter 2016, and references therein).

Our best-fitting secondary gas models predict CO production rates of ∼0.005 M_⊕ Myr⁻¹ and ∼0.035 M_⊕ Myr⁻¹ for 49 Cet and HD 32297, respectively. Considering the estimated mass-loss rates, 0.1 M_⊕ Myr⁻¹ for 49 Cet and 5.2 M_⊕ Myr⁻¹ for HD 32297 (see Section 5.2), these ${\dot{M}}_{\mathrm{CO}}$ rates can be reproduced with CO+CO₂ ice mass fractions of γ ∼ 5% and ∼0.7%, respectively. Using solar system comets as analogs, where this fraction ranges between 2% and 27% (Mumma & Charnley 2011; Matrà et al. 2017b), the required γ values are not unrealistically high. Considering the obtained CO mass production rates and the adopted evolutionary times, the total CO gas mass release of the disks around 49 Cet and HD 32297 throughout their evolution is ∼0.2 M_⊕ and ∼0.9 M_⊕, respectively.

5.3.3. Mass Evolution of the Gas Components

Figure 5 displays the mass evolution of the ¹²CO, ¹³CO, C¹⁸O gas, and carbon components in the best-fit secondary gas disk models of 49 Cet and HD 32297. According to this, the disk of HD 32297 reaches a CO mass of ∼0.01 M_⊕ after ∼0.5 million years of evolution and then appears as a CO-rich debris system for nearly all of its lifetime. In the 49 Cet system it takes ∼12 Myr to reach the same level of CO content. As these plots illustrate, the ${M}_{{}^{13}\mathrm{CO}}/{M}_{{}^{12}\mathrm{CO}}$ and ${M}_{{{\rm{C}}}^{18}{\rm{O}}}/{M}_{{}^{12}\mathrm{CO}}$ mass ratios show a significant variation during the studied evolutionary period. These changes are related to the isotope-selective photodissociation. As the ¹²CO line becomes more and more optically thick and the column density of neutral carbon increases, the ¹²CO line acquires strong shielding while the photodissociation timescale of rarer isotopologues, that remain optically thin, remains short. This leads to decreasing mass ratios. Then, as the column densities of neutral carbon, ¹³CO, and C¹⁸O increase further, and thus the shielding of less abundant CO isotopologues become more efficient, the mass ratios start to increase but do not reach the typical values found for the local interstellar matter (that we initially assumed for our icy planetesimals) even at the end of the studied period. As a result, the ¹²CO masses we obtain in these simulations (0.1 M_⊕ for 49 Cet and 0.9 M_⊕ for HD 32297) are about an order of magnitude higher than those we found in Section 4.2. In the case of HD 32297 the evolution is so rapid that the initial decreasing trend is less clearly outlined and significantly shorter than for 49 Cet. While in HD 32297 the CO mass exceeds that in 49 Cet by a large factor, our model predicts very similar C⁰ masses of 0.03 and 0.04 M_⊕ for HD 32297 and 49 Cet, respectively. These values are similar to the C⁰ mass found around HD 131835 (Kral et al. 2019).

Zoom In Zoom Out Reset image size

5.3.4. Evolution of Gas Surface Density Distribution

Figure 6 shows the actual radial surface density distributions of the ¹²CO, ¹³CO, and C¹⁸O molecules as well as of carbon atoms for three different evolutionary times using the best-fit model of the specific target. We predict that in both systems the radial distribution of C⁰ is more extended than those of the different CO molecules. This is because when the carbon gas spreads viscously away from the gas production site, the most displaced regions will become too thin to shield CO molecules. Drastic depletion of CO through photodissociation in these zones then results in a larger C⁰ to CO width. Future high angular resolution mapping of C i emission will make it possible to check these predictions and give further constraints on the models (see Kral et al. 2019).

Zoom In Zoom Out Reset image size

5.3.5. A Possible Primordial Origin of the Gas Component

Gas, if presents in sufficient amount, can have a significant influence on the dynamics of small debris particles (e.g., Takeuchi & Artymowicz 2001; Besla & Wu 2007; Richert et al. 2018). With their protoplanetary level of CO gas, the debris disks of 49 Cet and HD 32297 are among the best candidates where gas–dust interactions can really manifest. Kral et al. (2019) pointed out that such interactions could be important in the CO-rich disk around HD 131835 (see also Feldt et al. 2017). In order to assess how well the gas and grains are coupled in 49 Cet and HD 32297, we estimated the dimensionless stopping time that shows the timescale (in orbital units) needed for a dust particle to be entrained by the gas. According to Richert et al. (2018) the dimensionless stopping time (a.k.a. the Stokes number, St) can be computed as:

$$\begin{eqnarray}&&{St}\approx \displaystyle \frac{1}{\beta }{\left(\displaystyle \frac{{{\rm{\Sigma }}}_{{\rm{g}}}}{9\times {10}^{-5}{\rm{g}}{\mathrm{cm}}^{-2}}\right)}^{-1}\left(\displaystyle \frac{{L}_{* }}{{L}_{\odot }}\right){\left(\displaystyle \frac{{M}_{* }}{{M}_{\odot }}\right)}^{-1},\end{eqnarray} \tag{ 5 }$$

where β is the ratio of radiation to gravitational force, Σ_g is the gas surface density, while L_* and M_* are the luminosity and mass of the central star. A smaller Stokes number means stronger coupling between dust and gas. Assuming parent bodies on a circular orbit, β is equal to 0.5 for the smallest bound grains (Krivov 2010). Concerning Σ_g, we took the highest total (C⁰+CO) gas surface density of our secondary gas disk models (Figure 6). Then Equation (5) yields St ∼ 20 and ∼1 for the smallest bound grains in 49 Cet and HD 32297, respectively. These values, however, should be considered as upper limits since additional daughter products, H and O, from photodissociation of CO, CO₂, and H₂O, have not been taken into account in this calculation. In solar system comets the H₂O/CO mass ratio ranges between 3 and 200 (Mumma & Charnley 2011). Considering these additional products, like H and O, the St numbers could be at least a few times lower. This certainly results in St ≲ 1 for HD 32297, indicating that the smallest bound grains are at least marginally coupled to gas via drag even if the icy parent bodies have very low H₂O/CO ratios. In the 49 Cet system more H₂O daughter products, thus more water-rich planetesimals are needed to reach the same level of gas–dust coupling for such grains. Assuming compact spherical grains with radiation pressure efficiency of unity, the β parameter can be expressed as $\beta =0.574\left(\tfrac{{L}_{* }}{{L}_{\odot }}\right){\left(\tfrac{{M}_{* }}{{M}_{\odot }}\right)}^{-1}$ $\left(\tfrac{1\,\mu {\rm{m}}}{s}\right)\left(\tfrac{1\,{\rm{g}}\ {\mathrm{cm}}^{-3}}{\rho }\right)$, where s and ρ are the size and density of the grains (Burns et al. 1979). Taking a grain density of 2 g cm⁻³, the grain size β = 0.5 particles (blowout grain size) in the 49 Cet and HD 32297 systems are ∼5 μm and ∼3 μm.

The presence of gas has an influence on the dynamics of unbound grains with β > 0.5 as well. While they would normally escape from the system on the local dynamical timescale, their interaction with the gas, e.g., gas drag, can result in longer lifetimes leading to an accumulation of very small particles compared to a gas-free situation (Lieman-Sifry et al. 2016; Kral et al. 2019; Richert et al. 2018; Wyatt 2018). In both of our targets there are indications that trapping of small grains is really happening. Analyzing 12.5 μm and 17.9 μm Keck images of 49 Cet, Wahhaj et al. (2007) argued that the bulk of the observed mid-IR emission comes from grains with a size of ∼0.1 μm located between 30 and 60 au from the star. This overlaps with the region where gas was also found (Hughes et al. 2017). Interactions with this gas material can explain how such small particles, well below the blowout limit, can be retained in the disk.

Based on SED analysis Donaldson et al. (2013) derived a characteristic dust temperature of 83 K for the cold outer component of HD 32297. Supposing that dust grains act like blackbodies, this temperature would correspond to a dust ring located at a radius of 32 au from the star (${R}_{\mathrm{bb}}(\mathrm{au})={\left(\tfrac{{L}_{* }}{{L}_{\odot }}\right)}^{0.5}{\left(\tfrac{278.3{\rm{K}}}{{T}_{\mathrm{dust}}}\right)}^{2}$, Backman & Paresce 1993). High angular resolution millimeter imaging of the source, however, implies a substantially larger radius for the outer belt that was modeled by a planetesimal belt encompassed by a very extended halo (MacGregor et al. 2018). The radius of this planetesimal belt component is ∼100 au (average of the best-fit inner and outer belt radii). From this radius we obtain a ${\rm{\Gamma }}={R}_{\mathrm{disk}}/{R}_{\mathrm{bb}}={T}_{\mathrm{disk}}^{2}/{T}_{\mathrm{bb}}^{2}$ ratio of 3.1 that, because of the presence of the halo, should be considered as a lower limit. Pawellek et al. (2014) and Morales et al. (2016) measured lower Γ ratios for other resolved debris disks around stars with similar luminosity (6–10 L_☉), indicating that the outer disk of HD 32297 is unusually warm. The larger average dust temperature might be attributed to the presence of warm micron-sized particles that may have longer lifetimes than in gas-poor debris disks around similarly luminous host stars.

High angular resolution mid-IR thermal emission and near-IR scattered light images of HD 32297 provide further evidence for the presence of small dust. Based on their 11.2 μm observation, Fitzgerald et al. (2007) argued that warm dust grains emitting at this wavelength are situated beyond a radial separation of ≳065 (≳86 au), thus are colocated with the outer gas and dust belt mapped in the millimeter regime (MacGregor et al. 2018). By analyzing their 11.7 μm and 18.3 μm images of the source, Moerchen et al. (2007) also found that the emitting grains are mostly concentrated in the outer regions (≳94 au, considering the new Gaia DR2 distance) and derived a typical dust temperature of ∼180–190 K. Both works concluded that small submicron-sized particles were present at large distances from the star. It is worth noting this also means that the warm component inferred from the SED analysis (Donaldson et al. 2013) cannot be exclusively linked to an inner dust belt. By analyzing their multiband scattered light data on HD 32297 Bhowmik et al. (2019) measured a blue color in the near-IR spectrum that was interpreted as a signature of a significant presence of grains smaller than the blowout limit.

Several other CO-bearing debris disks show signatures of considerable amounts of small dust particles. By modeling resolved debris disks in their survey Lieman-Sifry et al. (2016) found evidence for the presence of grains smaller than the blowout size in two gaseous debris disks, HD 138813 and HD 156623. Moór et al. (2017) discovered that the disks of HD 121617 and HD 131488—similar to HD 32297—exhibit outstandingly high Γ factors, indicative of warm, small solids. The spatially resolved mid-IR image of HD 131835 also revealed very hot, supposedly submicron-sized particles via their strong emission in the outer disk (Hung et al. 2015). Thus, the existence of small grains in CO-rich debris disks seems to be a common phenomenon.

Due to gas drag, grains can move radially inward or outward in the disk. The direction, speed, and degree of this radial migration depend on many different factors, such as the density and temperature distribution of the gas component, the size of the particle, the strength of the radiative forces, and the collisional timescales. Grains at the lower end of the size distribution typically migrate outward (Takeuchi & Artymowicz 2001) potentially leading to a detachment between the ring of parent planetesimals and the small grain population. In a recent model of gas–dust interactions in optically thin disks, Richert et al. (2018) considered dust–gas drag with back reaction, photoelectric heating of gas by dust, and stellar radiation pressure of dust. In their simulations, they detected various structures like concentric rings, arcs, and spirals arms in the dust distribution, attributed to photoelectric instability. Though their results cannot be directly applied to our cases (e.g., they adopted a G-type star as a central object resulting in a substantially weaker radiation field), in their models with St ≲ 1 they saw several of the abovementioned features. We assume that in the case of HD 32297, where the Stokes number for β = 0.5 dust particles is also low (St ≲ 1), the same physical processes like gas drag–driven detachment and photoelectric instability could be operational. The resulting structures, if any, can be best observed in the distribution of small grains outlined by high spatial resolution scattered light images as pointed out in Kral et al. (2019).

So far we supposed that our two targets harbor secondary gas disks. As we already noted in Section 5.3, it cannot be excluded that the detected gas component is the residual of the primordial disk whose mass is dominated by H₂ molecules. Adopting the canonical H₂/CO abundance ratio of 10⁴ (i.e., a mass ratio of ∼700), the gas surface densities in our two disks could be several hundred times higher, meaning that even the >100 μm sized dust particles are well coupled to the gas. From the analysis of the 850 μm continuum ALMA image of 49 Cet, Hughes et al. (2017) found a single power-law disk model with a ring-like enhancement at 110 au in their best-fit model. If the gas component is primordial, this dust ring could be attributed to gas–dust interactions.