<figure><img src="image_1.png" /></figure>

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

Attila Moór¹, Quentin Kral², Péter Ábrahám¹ ID, Ágnes Kóspál¹,³ ID, Anne Dutrey⁴, Emmanuel Di Folco⁴, A. Meredith Hughes⁵,
Attila Juhász⁶, Ilaria Pascucci⁷ ID, and Nicole Pawellek¹,³ ID

¹ Konkoly Observatory, Research Centre for Astronomy and Earth Sciences, Hungarian Academy of Sciences, Konkoly-Thege Miklós út 15-17, 1121 Budapest, Hungary; moor@konkoly.hu

² LESIA, Observatoire de Paris, Université PSL, CNRS, Sorbonne Université, Univ. Paris Diderot, Sorbonne Paris Cité, 5 place Jules Janssen, F-92195 Meudon, France

³ Max Planck Institute for Astronomy, Königstuhl 17, D-69117 Heidelberg, Germany

⁴Laboratoire d'Astrophysique de Bordeaux, Univ. Bordeaux, CNRS, B18N, Allée Geoffroy Saint-Hilaire, F-33615 Pessac, France

⁵Department of Astronomy, Van Vleck Observatory, Wesleyan University, 96 Foss Hill Drive, Middletown, CT 06459, USA

⁶Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge CB3 0HA, UK

⁷Lunar and Planetary Laboratory, The University of Arizona, Tucson, AZ 85721, USA

Received 2019 January 16; revised 2019 September 2; accepted 2019 September 6; published 2019 October 16

Abstract

Previous observations revealed the existence of CO gas at nearly protoplanetary level in several dust-rich debris disks around young A-type stars. Here we used the Atacama Large Millimeter/submillimeter Array (ALMA) 7 m Array to measure $^{13}$CO and $^{18}$O emission toward two debris disks, 49 Cet and HD 32297, and detected similarly high CO content ($>0.01 M_{\oplus}$). These high CO masses imply a highly efficient shielding of CO molecules against stellar and interstellar ultraviolet photons. Adapting a recent secondary gas disk model that considers both shielding by carbon atoms and self-shielding of CO, we can explain the observed CO level in both systems. Based on the derived gas densities we suggest that, in the HD 32297 disk, dust and gas are coupled and the dynamics of small grains is affected by the gaseous component. For 49 Cet, the question of coupling remains undecided. We found that the main stellar and disk properties of 49 Cet and HD 32297 are very similar to those of previously identified debris disks with high CO content. These objects constitute together the first known representatives of shielded debris disks.

*Unified Astronomy Thesaurus concepts: Circumstellar disks (235); Debris disks (363); Circumstellar gas (238); A dwarf stars (3)*

**1. Introduction**

Circumstellar debris disks are mainly comprised of solids from micron- to planetesimal-size scales. Due to interactions with the stellar radiation field in such optically thin disks, dust grains are rapidly removed. Thus the smallest dust particles—that make these disks detectable for us at infrared (IR) wavelengths via their thermal emission—could not be leftovers from the primordial phase but should be second-generation grains stemming from collisional grinding of previously formed larger bodies (Wyatt 2008). These parent planetesimals are confined into a belt or belts (Hughes et al. 2018). If icy planetesimals and grains are present, then their collisions and sublimation may result in gas liberation as well (Zuckerman & Song 2012; Kral et al. 2017). Moreover, photodesorption of icy mantles can also contribute to gas release (Grigorieva et al. 2007). Similar to debris dust particles, gas produced in these ways is of second generation, rather than remnant gas from the protoplanetary phase.

Measuring the gas content of debris disks is crucial from
several different points of view. Observation of different gas
compounds offers a unique opportunity to constrain the volatile
composition of the parent icy planetesimals (Kral et al. 2016;
Matrà et al. 2017a, 2018b). The formation of rings, spirals, and
arcs in the distribution of small debris grains can be caused by
either the gravitational perturbation of an unseen planet or
interactions with a sufficient amount of gas in the system
(Takeuchi & Artymowicz 2001; Richert et al. 2018). To assess
which mechanism could be active, one needs to measure the
amount of gas. However, while the dust material of debris disks
is rather well studied, their possible gas content remains less

explored so far, especially due to observational difficulties.
Nevertheless, intensive investigations over the past decades
revealed gas emission in 19 debris disks. Though in a few
systems the gas component was detected through fine-structure
lines of neutral oxygen or ionized carbon (Riviere-Marichalar
et al. 2012, 2014), so far observation of low-level rotational
transitions of $^{12}$CO molecules has turned out to be the most
effective way to probe gas in debris disks. Thanks to
observations with single-dish radio telescopes and especially
with the Atacama Large Millimeter/submillimeter Array
(ALMA), as of today we know 17 CO-bearing debris disks.<fnref n="8" />
Thirteen disks surround A-type stars: 49 Cet (Zuckerman et al.
1995), HD 21997 (Moór et al. 2011), HD 32297 (Greaves et al.
2016), β Pic (Dent et al. 2014), HD 95086 (Booth et al. 2019),
HD 121191, HD 121617, HD 131488 (Moór et al. 2017),
HD 131835 (Moór et al. 2015a), HD 110058, HD 138813,
HD 156623 (Lieman-Sifry et al. 2016), and Fomalhaut (Matrà
et al. 2017a). Three disks are present in F-type systems,
HD 146897 (Lieman-Sifry et al. 2016), HD 181327 (Marino
et al. 2016), and η Crv (Marino et al. 2017), and recently Matrà
et al. (2019) revealed CO for the first time in a debris disk
hosted by an M-type star, TWA 7.

<footnote n="8">⁸ The CO-bearing disk around the young pre-main-sequence Herbig B9.5 Ve star, HD 141569, exhibits an IR excess at only a slightly higher level to that of the most dust-rich systems among the known gaseous debris disks, suggesting a relationship to them (e.g., Flaherty et al. 2016; Péricaud et al. 2016). However, in several aspects this system differs from genuine debris disks (Wyatt et al. 2015; White & Boley 2018). Because of its not fully clarified evolutionary state we will handle this object separately from the CO-bearing debris disk sample.</footnote>
The ¹²CO line J = 2–1 luminosities of the detected gaseous debris disks show a large spread of almost four orders of magnitude; the brightest eight of them (49 Cet, HD 21997, HD 32297, HD 121617, HD 131488, HD 131835, HD 138813, and HD 156623) exhibit line luminosities higher than ~10¹⁸ W, i.e., 6× higher than that of β Pic (Matrà et al. 2018b). Besides their strong ¹²CO emission, these eight disks share additional common properties. They all belong to the most dust-rich known debris systems with their dust fractional luminosities $f_d > 5 \times 10^{-4}$ (although there exist similarly dust-rich debris disks with no detectable CO, the most prominent example being HR 4796, Kennedy et al. 2018). They surround young (10–50 Myr) A-type stars, and the bulk of their detected gas material orbit at radial distances larger than ~20 au. Observations of the less abundant ¹³CO and C¹⁸O isotopologues in the disks around HD 21997, HD 121617, HD 131488, and HD 131835 brought a surprising result: the measured line ratios implied highly optically thick ¹²CO emission in all cases. The total CO gas masses of these systems, based on the analysis of the optically thin C¹⁸O line, are higher than 0.01 $M_⊕$ (Kóspál et al. 2013; Moór et al. 2017). The obtained CO masses are three to four orders of magnitude higher than those of the similarly young and dust-rich disks of β Pic and HD 181327 (Marino et al. 2016; Matrà et al. 2017b), and resemble more the CO gas quantity of protoplanetary disks around Herbig Ae stars (Moór et al. 2017; Kral et al. 2019).

Since the low dust content of debris disks does not provide efficient shielding against stellar and interstellar UV photons, gas molecules released from icy bodies are expected to be rapidly photodissociated. Unshielded second generation CO molecules are destroyed on a timescale of ~100 yr by high energy photons of the interstellar radiation field (Visser et al. 2009). This short lifetime is consistent with the low CO content of β Pic (Matrà et al. 2017b) or HD 181327 (Marino et al. 2016) but is in contradiction with the presence of the abovementioned substantially more CO-rich disks whose long-term maintenance would require an unrealistically high gas production rate. The existence of the latter disks can be understood if their CO molecules are shielded more effectively. The young age of these systems raises the possibility that the observed gas is not second generation but predominantly leftover from the primordial phase where remnant hydrogen molecules provide strong shielding to CO gas (Kóspál et al. 2013). Since the dust is thought to be second generation such disks would exhibit a hybrid nature (Kóspál et al. 2013; Péricaud et al. 2017). As an alternative solution Kral et al. (2019) recently demonstrated that shielding by neutral carbon gas (C⁰) produced mainly from photodissociation of CO and CO₂ in tandem with self-shielding by CO molecules can increase the lifetime of CO gas released from planetesimals significantly. They found that the gas quantity measured in all four abovementioned CO-rich debris disks can be explained by applying their shielded secondary gas disk model without the presence of any primordial material.

As part of our long-term project to study the properties and nature of all known young dust-rich debris disks exhibiting high ¹²CO line luminosity, in this work we investigated the CO content of two other representants, 49 Cet and HD 32297. In Section 2 we briefly summarize the currently known properties of our targets with special emphasis on previous studies regarding their gas material. Our new continuum and CO isotopologue observations carried out by the ALMA 7 m Array

are described in Section 3. Results from these observations including estimates for the dust and CO contents are presented in Section 4. Implications of these results are discussed in Section 5. Finally the outcomes of this investigation are summarized in Section 6.

## 2. Targets

49 Cet is an A1V-type star located at a distance of 57.0 ± 0.3 pc (Gaia DR2, Bailer-Jones et al. 2018; Lindegren et al. 2018; Gaia Collaboration et al. 2018). Based on its membership in the Argus moving group the estimated age of the star is 40–50 Myr (Zuckerman 2019). The debris disk of 49 Cet was discovered by Sadakane & Nishida (1986) through its strong far-infrared excess. The spectral energy distribution (SED) of the excess is best described by a two-component model (Roberge et al. 2013) with characteristic dust temperatures of 165 and 59 K and with a total fractional luminosity $f_d = 8.5 \times 10^{-4}$ (Holland et al. 2017). 49 Cet was the first debris system where millimeter CO emission was detected (Zuckerman et al. 1995). Since then several other gas species (e.g., C, C⁺, O) have been identified in the disk either in absorption (Roberge et al. 2014) or in emission (Roberge et al. 2013; Higuchi et al. 2017). Recently the spatial distribution of the ¹²CO gas and dust was investigated in great detail using high angular resolution (0″/4) submillimeter observations with the ALMA interferometer (Hughes et al. 2017). This study revealed a very extended, cold outer dust disk between 60 and 310 au. The CO emission was found to be axisymmetric and the bulk of the gas is located between 20 and 200 au radial distances. Since large dust grains emitting at millimeter wavelengths are not affected significantly by stellar radiation forces, the obtained dust map is thought to probe the location of the parent planetesimals. By analyzing near-infrared coronagraphic images, Choquet et al. (2017) showed that the disk's scattered light emission is also consistent with an axisymmetric model. Small dust particles explored by these measurements extend from ~65 to 250 au.

HD 32297, a young (<30 Myr old, Kalas 2005) A6-type main-sequence star, also harbors a well-known debris disk whose infrared excess was first identified by Silverstone (2000). Based on Gaia DR2 astrometry the star is located 132.3 ± 1.0 pc away (Bailer-Jones et al. 2018; Lindegren et al. 2018; Gaia Collaboration et al. 2018). Similar to 49 Cet, adequate fitting of the observed excess SED requires a two-temperature model. Donaldson et al. (2013) derived dust temperatures of 240 and 83 K for the warm and cold components. With its total fractional luminosity of $5.4 \times 10^{-3}$ (Kral et al. 2017), HD 32297 belongs to the most dust-rich known debris disks. The disk was resolved in scattered light at optical and near-IR wavelengths in several studies (Kalas 2005; Schneider et al. 2005; Debes et al. 2009; Boccaletti et al. 2012; Currie et al. 2012; Esposito et al. 2014; Rodigas et al. 2014), finding that it is nearly edge-on and extends far from the star (~1800 au,<fnref n="9" /> Schneider et al. 2014). These scattered light images also revealed a strong southwest-northeast brightness asymmetry. Using the JCMT radio telescope, Greaves et al. (2016) have recently detected CO (2–1) emission toward HD 32297. By analyzing their recent ALMA observations of the disk, MacGregor et al. (2018) found

<footnote n="9">⁹ Considering the new Gaia DR2 distance estimate instead of the Hipparcos-based 112 pc used by Schneider et al. (2014) in their paper.</footnote>
the CO gas to be colocated with large dust particles, whose spatial distribution was described by the combination of a planetesimal belt located between 78 and 122 au and a halo component that extends up to 440 au. CO is not the sole known gas component in the disk: HD 32297 is one of those rare debris disks where far-IR emission from ionized carbon was also detected with the Herschel Space Observatory (Donaldson et al. 2013). Prior to these observations, Redfield (2007) identified a stable circumstellar absorption component of NaI in the optical spectrum of the star.

We derived reddening values for the two stars by fitting their optical and near-IR photometry with ATLAS9 models (Castelli et al. 1997) using fixed $T_{\text{eff}}$, log g, and metallicity taken from Rebollido et al. (2018). The resulting $E(B-V)$ values are $0.024 \pm 0.009$ mag and $0.028 \pm 0.013$ mag for 49 Cet and HD 32297, respectively. Combining these data with the new Gaia DR2 distances we computed luminosities of $17.2 \pm 0.4 L_{\odot}$ for 49 Cet and $8.4 \pm 0.2 L_{\odot}$ for HD 32297.

### 3. Observations

49 Cet and HD 32297 were observed using the ALMA 7 m Array in stand-alone mode in the framework of our Cycle 4 project (2016.2.00200.S, PI: Á. Kóspál). We obtained Band 6 observations of the $J = 2-1$ rotational line of $^{12}$CO, $^{13}$CO, and $^{18}$O, as well as the 1.3 mm continuum. We used two spectral windows for the CO line studies. The $^{13}$CO and $^{18}$O lines were measured in a single baseband with a spectral resolution of 976.6 kHz and a bandwidth of 1875 MHz, while $^{12}$CO line was observed in another window with a 244.1 kHz spectral resolution and a bandwidth of 468.8 MHz. Two additional spectral windows with a bandwidth of 1875 MHz and centered at 217.0 and 233.5 GHz were defined to measure the continuum emission. Observations of 49 Cet were carried out in 2017 July with nine antennas providing projected baseline lengths between 9 and 45 m (6.7–34 kλ). The observation block began with pointing, bandpass and flux calibrations. Three quasars (J0522–3627, J0211+1051, J0116–1136) were used in the pointing calibrations, while quasar J0522–3627 and Uranus were targeted as bandpass and flux calibrators, respectively. Then, observations were alternated every ~6.5 minutes between 49 Cet and the quasar J0132–1654 that was used as a phase calibrator. Measurements of HD 32297 were taken in 2017 August using two observation blocks conducted on subsequent days (on August 18 and 19). Nine antennas were involved providing projected baseline lengths from 9 to 49 m (6.7–37 kλ). Observations were executed similar to that of 49 Cet, i.e., the measurement started with pointing, bandpass and flux calibrations followed by cycles between the target and a phase calibrator. J0522–3627 was used as bandpass calibrator in both blocks. In the first measurement Uranus was targeted as a flux calibrator and J0522–3627, J0211+1051, and J0457+0645 were used during the pointing sequences. In the second measurement the quasar J0510+1800 was the flux calibrator, while the pointing calibration was conducted using J0522–3627 and J0510+1800.

Both data sets were calibrated and flagged with the ALMA reduction tool Common Astronomy Software Applications (CASA v4.7.2; McMullin et al. 2007) using the pipeline script delivered with the data. For 49 Cet, data from one of the antennas were flagged entirely by this calibration pipeline. We used the CASA task uvcontsub to fit a first-order polynomial to line-free channels in the uv space and subtract

<table><caption>**Table 1**
Observational Data</caption><thead><tr><th>Parameters</th><th>49 Cet</th><th>HD 32297</th></tr></thead><tbody><tr><td colspan="3"><strong>Continuum</strong></td></tr><tr><td>Beam size (")</td><td>6"/6 &times; 5"/1</td><td>7"/8 &times; 4"/6</td></tr><tr><td>Beam PA (°)</td><td>-81°/6</td><td>+68°/5</td></tr><tr><td>rms (mJy beam<sup>-1</sup>)</td><td>0.25</td><td>0.11</td></tr><tr><td colspan="3"><strong><sup>12</sup>CO (2&ndash;1)</strong></td></tr><tr><td>Beam size (")</td><td>6"/4 &times; 4"/9</td><td>7"/6 &times; 4"/4</td></tr><tr><td>Beam PA (°)</td><td>-81°/5</td><td>+69°/1</td></tr><tr><td>rms (mJy beam<sup>-1</sup> chan.<sup>-1</sup>)</td><td>28.5</td><td>15.2</td></tr><tr><td colspan="3"><strong><sup>13</sup>CO (2&ndash;1)</strong></td></tr><tr><td>Beam size (")</td><td>6"/8 &times; 5"/2</td><td>7"/9 &times; 4"/6</td></tr><tr><td>Beam PA (°)</td><td>-80°/6</td><td>+68°/9</td></tr><tr><td>rms (mJy beam<sup>-1</sup> chan.<sup>-1</sup>)</td><td>17.7</td><td>9.0</td></tr><tr><td colspan="3"><strong>C<sup>18</sup>O (2&ndash;1)</strong></td></tr><tr><td>Beam size (")</td><td>6"/8 &times; 5"/2</td><td>8"/0 &times; 4"/6</td></tr><tr><td>Beam PA (°)</td><td>-80°/6</td><td>+68°/9</td></tr><tr><td>rms (mJy beam<sup>-1</sup> chan.<sup>-1</sup>)</td><td>11.9</td><td>7.7</td></tr></tbody></table>

<table><caption>Continuum flux densities and CO integrated line fluxes<sup>a</sup></caption><tr><td>$F_{\nu}$ at 1.33 mm (mJy)</td><td>5.3 ± 0.7(0.4)</td><td>3.4 ± 0.4(0.1)</td></tr><tr><td>$S_{12\text{CO}}$ (Jy km s<sup>-1</sup>)</td><td>3.87 ± 0.41(0.14)</td><td>1.05 ± 0.12(0.05)</td></tr><tr><td>$S_{13\text{CO}}$ (Jy km s<sup>-1</sup>)</td><td>1.68 ± 0.19(0.09)</td><td>0.50 ± 0.07(0.04)</td></tr><tr><td>$S_{\text{C}^{18}\text{O}}$ (Jy km s<sup>-1</sup>)</td><td>&lt;0.26</td><td>0.27 ± 0.05(0.04)</td></tr></table>

**Note.**

the continuum from the obtained CO emission. Then the tclean task was used to create spectral cubes of the three CO lines from the continuum-subtracted calibrated visibilities adopting natural weighting. The channel width was set to 0.7 km s<sup>-1</sup> for $^{12}$CO and to 1.4 km s<sup>-1</sup> for $^{13}$CO and $^{18}$O lines. To construct the continuum images we combined data from the two dedicated continuum spectral windows with line-free channels of the two other spectral windows. The continuum maps were also cleaned using natural weighting. The achieved sensitivities as well as the size of the resulting synthesized beams are summarized in Table 1.

### 4. Results

#### 4.1. Continuum Emission

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.

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

<footnote><sup>a</sup> The quoted uncertainties include the calibration uncertainties; the numbers in parentheses show the measurement errors.</footnote>
<figure><img src="image_2.png" /><figcaption>Figure 1. Continuum and zeroth moment maps for three different CO isotopologues obtained with the ALMA 7 m Array in Band 6 for 49 Cet (first row) and HD 32297 (second row) as well as the spatially integrated spectra of the $^{12}$CO, $^{13}$CO, and C$^{18}$O emission (third row).</figcaption></figure>

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)

$$M_d = \frac{F_\nu d^2}{B_\nu(T_{\text{dust}}) \kappa_\nu}, \quad (1)$$

where d is the distance, $B_\nu(T_{\text{dust}})$ is the Planck function at 1.33 mm for a characteristic dust temperature $T_{\text{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_{\text{dust}} = 56$ K. Following Beckwith et al. (1990), we adopted a value of 2.3 cm$^2$ g$^{-1}$ for the dust opacity. As a result we obtained dust masses of 0.15 ± 0.02 $M_\oplus$ and 0.64 ± 0.07 $M_\oplus$ 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).

## 4.2. CO Line Emission

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 5.7 (the positional uncertainty, estimated as the ratio of the beam size to achieved S/N, is ~1.7). 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¹²CO/S¹³CO = 2.3 ± 0.2 and S¹³CO/S¹⁸O > 6.5 for 49 Cet and S¹²CO/S¹³CO = 2.1 ± 0.2 and S¹³CO/S¹⁸O = 1.8 ± 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¹²CO/S¹³CO ~ 85 – 91 and S¹³CO/S¹⁸CO ~ 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 γC²O isotopologue as:

$$M_{\gamma C^{\prime \prime} O} = 4\pi m d^2 \frac{S_{21}}{x_2 h \nu_{21} A_{21}}, \quad (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+11/-4 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¹³C¹⁶O = (1.5 ± 0.2) × 10⁻⁴ M⊕ for 49 Cet and M¹²C¹⁸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\text{C}^{16}\text{O}} = (1.11 \pm 0.13) \times 10^{-2} M_{\oplus}$ and $(7.4 \pm 1.3) \times 10^{-2} M_{\oplus}$ 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).

# 5. Discussion

## 5.1. Young Debris Disks with Large CO Content

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_{\oplus}$. 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_{\oplus}$, 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.

## 5.2. Shielding of CO Gas in 49 Cet and HD 32297

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 necessarily 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

<figure><img src="image_3.png" /><figcaption>**Figure 2.** CO masses of circumstellar disks as a function of age. Dust masses are also shown by the colors of the symbols. Debris disks are plotted by circles, while disks around Herbig Ae stars are displayed by squares. Data for 49 Cet and HD 32297 are taken from Sections 4.1 and 4.2. For the other four CO-rich debris disks we used data from Moór et al. (2017), but considering their new Gaia DR2-based distances from Bailer-Jones et al. (2018). To estimate the CO masses of protoplanetary disks we followed the outline described in Section 4.2. Apart from the case of HD 141569, where a ¹³CO observation was used (Péricaud 2016), all our calculations are based on C¹⁸O line fluxes that are taken from the literature (Favre et al. 2013; Rosenfeld et al. 2013; Ansdell et al. 2016; Fedele et al. 2017; Hales et al. 2018; Kastner et al. 2018; van der Plas et al. 2019) or derived by processing data from the ALMA archive. Dust masses are derived from (sub)millimeter continuum measurements of the targets (Henning et al. 1994; Sylvester et al. 1996; Raman et al. 2006; Isella et al. 2007; Chapillon et al. 2008; Andrews et al. 2012; Meeus et al. 2012; Ansdell et al. 2016; Hales et al. 2018; Kastner et al. 2018), assuming optically thin emission. For Herbig Ae stars we adopted gas and dust temperatures of 50 K in these calculations.</figcaption></figure>

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

$$M_{\text{CO}} = \dot{M}_{\text{CO}}t_{\text{ph}}\epsilon_{\text{CO}} = \dot{M}_{\text{loss}}\gamma t_{\text{ph}}\epsilon_{\text{CO}}, \quad (3)$$

where $\dot{M}_{\text{loss}}$ is the mass-loss rate in the collisional cascade, $\gamma$ is the CO+CO₂ ice mass fraction of planetesimals, $t_{\text{ph}}$ is the photodissociation timescale of unshielded CO gas, while $\epsilon_{\text{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_{\text{ph}}$, $\gamma$, and $\dot{M}_{\text{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 $\gamma$ 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)):

$$\dot{M}_{\text{loss}}(M_{\oplus} \text{ Myr}^{-1}) = 1200 \left(\frac{R}{\text{au}}\right)^{1.5} \left(\frac{\Delta R}{\text{au}}\right)^{-1} f_d^2 \left(\frac{L_{\odot}}{L_{\odot}}\right) \left(\frac{M_{\oplus}}{M_{\odot}}\right)^{-0.5} . \quad (4)$$

Based on MacGregor et al. (2018), for the outer belt of HD 32297 we adopted a radius ($R$) and width ($\Delta 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 $\Delta 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_{\oplus} \text{ Myr}^{-1}$ for 49 Cet and HD 32297, respectively. Applying these values in Equation (3) yielded shielding factors of $\epsilon_{\text{CO}} \sim 3600$ and $\sim 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 $\epsilon_{\text{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 $^{12}\text{CO}$ column densities higher than $10^{15} \text{ cm}^{-2}$ self-shielding can reduce significantly the photodissociation rate of $^{12}\text{CO}$ molecules located in deeper regions. By interpolating in their Table 5 we found that to achieve the above shielding factors the vertical $^{12}\text{CO}$ column density in the

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Shielding factors necessary to explain the CO content of the six CO-rich debris disks and the disk of $\beta$ Pic (Section 5.2). 49 Cet and HD 32297, the targets of our study, are plotted with blue symbols. CO column densities corresponding to the given shielding factors are also indicated on the right-hand side. These column densities were computed by interpolating in Table 5 of Visser et al. (2009). We note that in the case of HD 21997, where the required shielding factor is extremely high, we had to use extrapolation to obtain the relevant CO column density.</figcaption></figure>

disks around 49 Cet and HD 32297 must exceed $\sim 1.7 \times 10^{19} \text{ cm}^{-2}$ and $\sim 1.5 \times 10^{18} \text{ cm}^{-2}$, 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 $\sim 10^{18} \text{ cm}^{-2}$ 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 $^{12}\text{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 $\beta$ Pic. Stellar and disk parameters needed for calculations of $\dot{M}_{\text{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}_{\text{loss}}$ values range between 0.04 and $4.2 M_{\oplus} \text{ Myr}^{-1}$. The CO mass estimates are from Figure 2, except for $\beta$ Pic where it was taken from Matrà et al. (2018b). As the plot shows, while the CO content of $\beta$ Pic can be explained with an unshielded or weakly shielded model, all CO-rich disks require substantial shielding with $\epsilon_{\text{CO}} > 100$, i.e., following the proposal of Kral et al. (2019) the latter objects can be rightfully called *shielded debris disks*.

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 $^{12}\text{CO}$ molecules whose total mass was estimated from optically thin $C^{18}\text{O}$ or $^{13}\text{CO}$ observations assuming abundance ratios to be
equal to the [$^{16}$O]/[$^{18}$O] and [ $^{12}$C]/[ $^{13}$C] isotope ratios measured in the local interstellar matter (see Section 4.2 and Moór et al. 2017). UV absorption lines of $^{13}$CO and C$^{18}$O are only partly overlapping with those of $^{12}$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$^{18}$O molecules—depending on the $^{12}$CO column densities— could be up to ~15 times higher than that of $^{12}$CO. By leading to higher $^{12}$CO/$^{13}$CO and $^{12}$CO/C$^{18}$O ratios this isotope-selective photodissociation would result in higher $^{12}$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 $^{13}$CO or C$^{18}$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 $^{12}$CO line luminosities than 49 Cet and HD 21997. Though, because of the lack of $^{13}$CO and/or C$^{18}$O observations, their CO mass is less constrained, the lower mass limits derived from the $^{12}$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 \times 10^{-3} M_{\oplus}$ (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 $\beta$ 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.

### 5.3. Origin of Gas in 49 Cet and HD 32297

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$^0$ 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

In the following, we evaluate whether this shielded secondary gas disk model can reproduce the observed quantities of 49 Cet and HD 32297 as well. The model presented in Kral et al. (2019) can be used to model gas released from exometary bodies and it is now able to follow the evolution of massive CO disks thanks to the inclusion of a new ingredient in their modeling: C$^0$ shielding of CO. It is also the first model to follow the evolution of CO together with C$^0$. CO is originated from the cometary bodies at the parent belt location at a rate $M_{CO}$. CO is then initially photodissociated in ~100 yr when no shielding occurs, which forms a gas disk made of atomic carbon and oxygen. In this new model, Kral et al. (2019) follow the evolution of carbon taking into account that at each time step there is some new carbon coming from CO destruction and that carbon depletes at a given radial location R because it viscously spreads. The viscosity $\nu$ is parameterized with an $\alpha$ parameter such that $\nu = \alpha c_s H$, where $c_s$ is the sound speed and H the scale height of the disk. In turn, when and if enough carbon accumulates, the model includes a back reaction on the CO photodissociation timescale that becomes longer owing to carbon shielding that can be very efficient. When the CO photodissociation timescale becomes longer than the viscous timescale (roughly equal to $R^2/\nu$), the CO gas disk has also time to viscously spread, which is included in the model.

We note that in debris disks $\alpha$ values could be very different from those in protoplanetary disks for several reasons. Kral & Latter (2016) found that the magnetorotational instability could work in a debris disk environment and could even lead to high $\alpha$ values because of the high ionization fraction that can be reached in unshielded disks (mostly made of ionized carbon, such as in $\beta$ Pic, Cataldi et al. 2018). In HD 32297 and 49 Cet, the ionization fraction could be much lower because of the strong shielding leading to drastically lower $\alpha$ values.

Stellar radiation pressure can have a substantial effect on gas particles, leading to the blowout of certain species (Fernández et al. 2006). Especially around the more luminous 49 Cet, carbon atoms can also be subject to radiation forces. However, as Kral et al. (2017) demonstrated, when a disk becomes optically thick to UV radiation in the radial direction, the efficiency of the star's radiation pressure decreases dramatically. They found that in disks with a CO production rate of $>10^{-4} M_{\oplus} \text{Myr}^{-1}$—which is significantly lower than in our targets (see Section 5.3.2)—carbon can be protected from expulsion. Therefore this effect was not considered in our modeling.

The model we use in this paper is based on the Kral et al. (2019) model described above but it has been improved. Now, instead of just following the evolution at the parent belt location, we follow the evolution at all radial locations. This new 1D model (which will be presented in detail in a forthcoming paper) thus solves in parallel the viscous diffusion equations (see Lynden-Bell & Pringle 1974, and also how we proceeded in Kral et al. 2016) for CO and C$^0$ together with the abovementioned source and sink terms. We also follow the viscous evolution of $^{13}$CO and C$^{18}$O (which we assume to be released at the ISM proportion compared to $^{12}$CO) together with the fact that self-shielding differs for these different
<figure><img src="image_5.png" /><figcaption>**Figure 4.** $^{13}$CO (49 Cet, left) and C$^{18}$O (HD 32297, right) gas masses derived in our secondary gas disk models for different CO production rates $\dot{M}_{\text{CO}}$ ($M_{\oplus}$ Myr$^{-1}$) and different $\alpha$ viscosities. Solid black contours correspond to the observed $^{13}$CO and C$^{18}$O masses, while the black dashed contours mark the $3\sigma$ uncertainties. Blue lines denote models with different inner CO disk radii ($R_{\text{in,CO}}$). Solid and dashed green lines in the left panel outline models producing the measured C I line flux and its higher $3\sigma$ uncertainty. The light gray lines in the right panel show models reproducing the observed $^{12}$CO line flux assuming disk inclinations of 84° (solid line) and 88° (dashed). In the case of 49 Cet the acceptable region in the parameter space is defined by the $3\sigma$ upper limit of the C I line flux (green dashed line at 33 Jy km s$^{-1}$) and by the $3\sigma$ upper limit of $R_{\text{in,CO}}$ (blue dashed line at 50 au) and by the $\pm 3\sigma$ uncertainties of the $^{13}$CO mass (dashed black lines). In the case of HD 32297 the acceptable region is defined by the $3\sigma$ lower limit of the $R_{\text{in,CO}}$ (blue dashed line at 55 au) and by the $\pm 3\sigma$ uncertainties of the C$^{18}$O mass (dashed black lines).</figcaption></figure>

isotopes (Visser et al. 2009). We therefore take into account the isotope-selective photodissociation for the first time in a debris disk model. Another upgrade is that this new 1D code calculates the ionization fraction in the disk at every radial location taking into account the UV radiation from the interstellar medium as well as the star. This code is thus ideally suited to model disks where spatial information on the gas and dust locations have previously been obtained from resolved ALMA observations. We will now describe how we used it for the specific systems around 49 Cet and HD 32297.

### 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 $^{13}$CO and C$^{18}$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}_{\text{CO}}$ and $\alpha$. The grid is logarithmic for both parameters and comprised of 10 elements going from $10^{-3}$ to $10^{-1} M_{\oplus}$ Myr$^{-1}$ for $\dot{M}_{\text{CO}}$ and from $10^{-5}$ to $10^{-2}$ for $\alpha$. For each model, we compute the total $^{13}$CO (for 49 Cet) or C$^{18}$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 $^{13}$CO and C$^{18}$O masses, while the black dashed contours mark the $3\sigma$ uncertainties (Table 1).

For 49 Cet (left panel), we found that models with $\dot{M}_{\text{CO}}$ ranging from 0.004 to 0.01 $M_{\oplus}$ Myr$^{-1}$ and with any $\alpha$ can explain the measured $^{13}$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 $^{3}$P$_{0}$–$^{3}$P$_{1}$ 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 \pm 2.3$ Jy km s$^{-1}$ for the integrated line flux.<fnref n="10" /> Models with different $\alpha$ and $\dot{M}_{\text{CO}}$ result in quite different inner CO disk radii ($R_{\text{in,CO}}$) and C I emissions. The solid blue line in Figure 4 (left) corresponds to models where $R_{\text{in,CO}} = 20$ 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\sigma$ level with the observed parameters. We found that there exists a region in the parameter space, at $2 \times 10^{-5} \lesssim \alpha \lesssim 9 \times 10^{-5}$ and $0.004 \lesssim \dot{M}_{\text{CO}} \lesssim 0.006 M_{\oplus}$ Myr$^{-1}$, where all three observational constraints, represented by the blue (inner CO disk radius), green (C I emission), and black (C $^{13}$CO mass) lines, are simultaneously satisfied within their $3\sigma$ uncertainties (dashed lines). As a representative

<footnote n="10">¹⁰ Higuchi et al. (2017) quoted the line flux in unit of Kkms$^{-1}$ (0.45 ± 0.04 Kkms$^{-1}$). To convert K to Jy we followed the description at https://alma.mtk.nao.ac.jp/aste/cfp2012/note.html adopting a half-power beamwidth of 17" based on Higuchi et al. (2017).</footnote>
<figure><img src="image_6.png" /><figcaption>**Figure 5.** The mass evolution of ¹²CO, ¹³CO, C¹⁸O gas and C⁰ gas components in representative secondary gas disk models of 49 Cet (left) and HD 32297 (right). Variation of the *M*₁₃CO/*M*₁₂CO and *M*₄¹⁸O/*M*₁₂CO mass ratios during the disk evolution are also shown. The models predict substantial variations in the mass ratios. There is an initial (very short for HD 32297) period when the mass ratios decrease. The ¹²CO line becomes gradually optically thick leading to stronger self-shielding of the ¹²CO molecules. At the same time the optically thin rarer isotopologues remain less efficiently shielded (isotope-selective photodissociation). This leads to their decreasing mass ratios with respect to ¹²CO. Later, the shielding of ¹³CO and C¹⁸O also becomes more efficient due to the increasing carbon and CO column densities and the mass ratios start to increase.</figcaption></figure>

grid point for the subsequent simulation, we adopted
$\dot{M}_{\text{CO}} \sim 4.6 \times 10^{-3} M_{\oplus} \text{Myr}^{-1}$ and $\alpha \sim 4.6 \times 10^{-5}$.

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ᵢₘₐ = 70 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ᵢₘₐ = 55 au models. The grid point located closest to these crossings has $\dot{M}_{\text{CO}} \sim 3.6 \times 10^{-2} M_{\oplus} \text{Myr}^{-1}$ and $\alpha \sim 1.5 \times 10^{-5}$ 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
~100($\frac{T}{20\,\text{K}})^{-1/2}$ ($\frac{\Sigma_n}{10^{-5}\,\text{g cm}^{-2}}$) ($\frac{f}{10^{-2}}$), where T is the gas temper-
ature, $\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^{-2}$ (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<sub>⊕</sub> Myr<sup>-1</sup> and ~0.035 M<sub>⊕</sub> Myr<sup>-1</sup> for 49 Cet
and HD 32297, respectively. Considering the estimated mass-
loss rates, 0.1 M<sub>⊕</sub> Myr<sup>-1</sup> for 49 Cet and 5.2 M<sub>⊕</sub> Myr<sup>-1</sup> for
HD 32297 (see Section 5.2), these M<sub>CO</sub> rates can be reproduced
with CO+CO<sub>2</sub> ice mass fractions of γ ~ 5% and ~0.7%,
respectively. Using solar system comets as analogs, where this
fraction ranges between 2% and 27% (Mumma & Charn-
ley 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<sub>⊕</sub> and ~0.9 M<sub>⊕</sub>,
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<sub>⊕</sub> 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<sub>13CO</sub>/M<sub>12CO</sub> and M<sub>48O</sub>/M<sub>12CO</sub> mass ratios show
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Radial surface density distributions of $^{12}$CO, $^{13}$CO, and C$^{18}$O molecules and C$^0$ atoms at three different evolutionary phases in the best-fit secondary gas disk models of 49 Cet (left) and HD 32297 (right). The C$^0$ analytics line (orange) is the level C$^0$ would reach at steady state for the given M$_{CO}$ and α values without shielding (assuming an ionization fraction of zero), see Equation (1) in Kral et al. (2019).</figcaption></figure>

<figure><img src="image_8.png" /></figure>

a significant variation during the studied evolutionary period. These changes are related to the isotope-selective photodissociation. As the $^{12}$CO line becomes more and more optically thick and the column density of neutral carbon increases, the $^{12}$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, $^{13}$CO, and C$^{18}$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 $^{12}$CO masses we obtain in these simulations (0.1 M$_{\oplus}$ for 49 Cet and 0.9 M$_{\oplus}$ 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$^0$ masses of 0.03 and 0.04 M$_{\oplus}$ for HD 32297 and 49 Cet, respectively. These values are similar to the C$^0$ mass found around HD 131835 (Kral et al. 2019).

### 5.3.4. Evolution of Gas Surface Density Distribution

Figure 6 shows the actual radial surface density distributions of the $^{12}$CO, $^{13}$CO, and C$^{18}$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$^0$ 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$^0$ 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).

### 5.3.5. A Possible Primordial Origin of the Gas Component

We found that the shielded secondary gas disk model offers a viable explanation for the prominent CO content of the young debris disks around 49 Cet and HD 32297. Although this finding implies that there is no need for additional large amounts of primordial H$_2$ molecules to explain the shielding of CO gas, we cannot rule out their presence based on our current results. To clarify this question, further observations are needed. One possible way to constrain the overall composition of the gas is to derive the scale height of the disk that depends on the gas temperature and the mean molecular weight (Hughes et al. 2017; Kral et al. 2019). By measuring the kinetic gas temperature and the vertical extent of the disk we can thus estimate the mean molecular weight that is expected to be very different in the primordial/secondary gas scenarios due to the different gas composition. Hughes et al. (2017) has already applied this method for 49 Cet. Though the vertical size of the disk remained unresolved in their observation, the obtained upper limit favors a higher mean molecular weight, i.e., a secondary gas composition. Due to its nearly edge-on orientation and brightness in CO lines the disk around HD 32297 is ideal for such observation. Though the detection of cold H$_2$ gas in emission is very difficult, its electron transitions can be observed in absorption at far-ultraviolet wavelengths offering an opportunity to measure the H$_2$ column density in disks seen close to edge-on. By analyzing spectroscopic observations obtained with the *Far Ultraviolet Spectroscopic Explorer* satellite, Lecavelier des Etangs et al. (2001) were able to set an upper limit of $<10^{18}$ cm$^{-2}$ on the column density of H$_2$ molecules in the gaseous edge-on disk of β Pic. Comparing this result to the column density of CO molecules detected in absorption resulted in a CO to H$_2$ ratio of $>6 \times 10^{-4}$ indicating—in good accordance with conclusions from other observations (Fernández et al. 2006; Kral et al. 2016; Matrà et al. 2017b; Cataldi et al. 2018)—that the gas is of second generation in β Pic. Later, Martin-Zaïdi et al. (2008) reported an even stricter upper limit of N$_{H_2}$ < $2.6 \times 10^{17}$ cm$^{-2}$ for β Pic. No similar data are available for any of the shielded
debris disks. Among them, HD 32297 harbors the most edge-on disk ($83.6^\circ < i < 88^\circ$) offering a special opportunity to directly measure the CO/H$_2$ ratio using the same method in the future. For outlines of additional possible observational methods that could be used to further constrain the origin of gas in shielded debris disks we direct the reader to Kral et al. (2019).

### 5.4. Gas-dust Interaction

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:

$$St \approx \frac{1}{\beta} \left( \frac{\Sigma_g}{9 \times 10^{-5} \text{ g cm}^{-2}} \right)^{-1} \left( \frac{L_*}{L_{\odot}} \right) \left( \frac{M_*}{M_{\odot}} \right)^{-1}, \quad (5)$$

where $\beta$ is the ratio of radiation to gravitational force, $\Sigma_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, $\beta$ is equal to 0.5 for the smallest bound grains (Krivov 2010). Concerning $\Sigma_g$, we took the highest total (C$^0$+CO) gas surface density of our secondary gas disk models (Figure 6). Then Equation (5) yields $St \sim 20$ and $\sim 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$_2$, and H$_2$O, have not been taken into account in this calculation. In solar system comets the H$_2$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 \lesssim 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$_2$O/CO ratios. In the 49 Cet system more H$_2$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 $\beta$ parameter can be expressed as $\beta = 0.574 \left(\frac{L_*}{L_{\odot}}\right) \left(\frac{M_*}{M_{\odot}}\right)^{-1} \left(\frac{1 \mu\text{m}}{s}\right) \left(\frac{1 \text{ g cm}^{-3}}{\rho}\right)$, where s and $\rho$ are the size and density of the grains (Burns et al. 1979). Taking a grain density of 2 g cm$^{-3}$, the grain size $\beta = 0.5$ particles (blowout grain size) in the 49 Cet and HD 32297 systems are $\sim 5 \mu\text{m}$ and $\sim 3 \mu\text{m}$.

The presence of gas has an influence on the dynamics of unbound grains with $\beta > 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 $\mu\text{m}$ and 17.9 $\mu\text{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 $\sim 0.1 \mu\text{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_{bb}(au) = \left(\frac{L_*}{L_{\odot}}\right)^{0.5} \left(\frac{278.3 \text{ K}}{T_{\text{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 $\sim 100$ au (average of the best-fit inner and outer belt radii). From this radius we obtain a $\Gamma = R_{\text{disk}}/R_{\text{bb}} = T_{\text{disk}}^2/T_{\text{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 $\Gamma$ ratios for other resolved debris disks around stars with similar luminosity ($6-10 L_{\odot}$), 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 $\mu\text{m}$ observation, Fitzgerald et al. (2007) argued that warm dust grains emitting at this wavelength are situated beyond a radial separation of $\gtrsim 0.65$ ($\gtrsim 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 $\mu\text{m}$ and 18.3 $\mu\text{m}$ images of the source, Moerchen et al. (2007) also found that the emitting grains are mostly concentrated in the outer regions ($\gtrsim 94$ au, considering the new Gaia DR2 distance) and derived a typical dust temperature of $\sim 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 \lesssim 1$ they saw several of the abovementioned features. We assume that in the case of HD 32297, where the Stokes number for $\beta = 0.5$ dust particles is also low ($St \lesssim 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_2$ molecules. Adopting the canonical $H_2$/CO abundance ratio of $10^4$ (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.

## 6. Summary

We used the ALMA 7 m Array to measure the $^{13}$CO and $^{18}$O isotopologues toward young debris disks around 49 Cet and HD 32297. Our observations revealed unexpectedly high CO gas contents comparable to the levels in protoplanetary disks. Though collisions and evaporation of icy bodies in a debris disk can naturally lead to gas production, the presence of the observed huge amount of CO gas in these systems needs a special environment in which the released CO molecules are shielded against UV photons from the stellar and interstellar radiation field.

Previous secondary gas models were not able to account for debris systems with such high CO masses putting forward the concept of hybrid disks where the gas component is dominated by residual primordial $H_2$ molecules. To examine whether the

disks of 49 Cet and HD 32297 could be fully secondary we utilized the shielded secondary gas model proposed recently by Kral et al. (2019). We found that by assuming a CO production rate of ~0.005–0.03 $M_\oplus$ Myr$^{-1}$ and slow viscous evolution with $\alpha$ of <10$^{-4}$, shielding by carbon atoms and self-shielding by CO together could be so effective that the obtained $^{13}$CO (for 49 Cet) and $^{18}$O (for HD 32297) masses can be explained in the framework of this model. In the case of 49 Cet the line flux of the predicted C I component was also found to be consistent with the observed one. The necessary CO production rates in the two systems are consistent with their dust production rates derived from observations, if we adopt ice mass fractions typical in solar system comets. Moreover, the total CO mass released from the planetesimals over the lifetime of the disks does not exceed the CO content of known nearby protoplanetary disks around Herbig Ae stars.

In the presence of sufficiently dense gas, dust particles and gas can interact with each other in various ways. Based on gas surface densities derived from our new shielded secondary gas models we argue that small grains with sizes close to the blowout limit in the disk of HD 32297 are at least moderately coupled to the gas and likely subject to such interactions. In accordance with this, several observational results indicate that the dust material of the gaseous outer ring of HD 32297 is unusually warm suggesting an overabundance of small solids in this region. The high gas mass of HD 32297 makes this source one of the most promising targets to study the possible outcomes of gas–dust interactions in optically thin dust disks. Though, based on our secondary gas model, the gas–dust coupling might be weaker in the 49 Cet system, there are also indications for the presence of a considerable amount of particles with sizes smaller than the blowout limit in this system. The retention of these particles suggests that gas–dust interaction can happen at some level in this system, too. We note that if the gas has a primordial origin then even large grains with sizes >100 µm could be entrained by the gas.

49 Cet and HD 32297 are in many ways similar to HD 21997, HD 121617, HD 131488, and HD 131835. The central stars in these systems are younger than 50 Myr, have A spectral type, and have no known stellar companion. Their debris disks have dust-rich ($f_d > 5 \times 10^{-4}$) outer belts at radii larger than 30 au. These belts are at least partly colocated with a gas component, having a protoplanetary disk–like CO mass of $M_{CO} > 0.01 M_\oplus$. Though the exceptional gas content of all six disks can be explained in the framework of the shielded secondary disk scenario, we cannot rule out that in some of them the very effective shielding of CO is due to leftover primordial $H_2$ gas. Future observations are needed for a final clarification of this question.

We thank the anonymous referee for providing insightful comments that helped us to improve the quality of this paper. We also thank Alycia Weinberger, Meredith MacGregor, and Gerrit van der Plas for sharing their results with us prior to publication. This paper makes use of the following ALMA data: ADS/JAO.ALMA#2017.2.00200.S. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada) and NSC and ASIAA (Taiwan) and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ. This work has made use of data from the European Space Agency (ESA)
mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/ consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. Our work was supported by the Hungarian OTKA grants K119993 and KH130526.

Facility: ALMA.

Software: CASA (v4.7.2 McMullin et al. 2007), LIME (Brinch & Hogerheijde 2010).

## ORCID iDs

Péter Ábrahám https://orcid.org/0000-0001-6015-646X

Ágnes Kóspál https://orcid.org/0000-0001-7157-6275

Ilaria Pascucci https://orcid.org/0000-0001-7962-1683

Nicole Pawellek https://orcid.org/0000-0002-9385-9820

## References

Andrews, S. M., Wilner, D. J., Hughes, A. M., et al. 2012, ApJ, 744, 162

Ansdell, M., Williams, J. P., van der Marel, N., et al. 2016, ApJ, 828, 46

Backman, D. E., & Paresce, F. 1993, in Protostars and Planets III, ed. E. H. Levy & J. I. Lunine (Tucson, AZ: Univ. Arizona Press), 1253

Bailer-Jones, C. A. L., Rybizki, J., Fouesneau, M., Mantelet, G., & Andrae, R. 2018, AJ, 156, 58

Beckwith, S. V. W., Sargent, A. I., Chini, R. S., & Guesten, R. 1990, AJ, 99, 924

Besla, G., & Wu, Y. 2007, ApJ, 655, 528

Bhowmik, T., Boccaletti, A., Thébault, P., et al. 2019, A&A, 630, 85

Boccaletti, A., Augereau, J.-C., Lagrange, A.-M., et al. 2012, A&A, 544, A85

Booth, M., Matrà, L., Su, K. Y. L., et al. 2019, MNRAS, 482, 3443

Brinch, C., & Hogerheijde, M. R. 2010, A&A, 523, A25

Burns, J. A., Lamy, P. L., & Soter, S. 1979, Icar, 40, 1

Castelli, F., Gratton, R. G., & Kurucz, R. L. 1997, A&A, 318, 841

Cataldi, G., Brandeker, A., Wu, Y., et al. 2018, ApJ, 861, 72

Chapillon, E., Guilloteau, S., Dutrey, A., & Piétu, V. 2008, A&A, 488, 565

Choquet, É., Milli, J., Wahhaj, Z., et al. 2017, ApJL, 834, L12

Currie, T., Rodigas, T. J., Debes, J., et al. 2012, ApJ, 757, 28

Debes, J. H., Weinberger, A. J., & Kuchner, M. J. 2009, ApJ, 702, 318

Dent, W. R. F., Wyatt, M. C., Roberge, A., et al. 2014, Sci, 343, 1490

Donaldson, J. K., Lebreton, J., Roberge, A., Augereau, J.-C., & Krivov, A. V. 2013, ApJ, 772, 17

Esposito, T. M., Fitzgerald, M. P., Graham, J. R., & Kalas, P. 2014, ApJ, 780, 25

Favre, C., Cleeves, L. I., Bergin, E. A., Qi, C., & Blake, G. A. 2013, ApJL, 776, L38

Fedele, D., Carney, M., Hogerheijde, M. R., et al. 2017, A&A, 600, A72

Feldt, M., Olofsson, J., Boccaletti, A., et al. 2017, A&A, 601, A7

Fernández, R., Brandeker, A., & Wu, Y. 2006, ApJ, 643, 509

Fitzgerald, M. P., Kalas, P. G., & Graham, J. R. 2007, ApJ, 670, 557

Flaherty, K. M., Hughes, A. M., Andrews, S. M., et al. 2016, ApJ, 818, 97

Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2018, A&A, 616, A1

Greaves, J. S., Holland, W. S., Matthews, B. C., et al. 2016, MNRAS, 461, 3910

Grigorieva, A., Thébault, P., Artymowicz, P., & Brandeker, A. 2007, A&A, 475, 755

Hales, A. S., Pérez, S., Saito, M., et al. 2018, ApJ, 859, 111

Henning, T., Launhardt, R., Steinacker, J., & Thamm, E. 1994, A&A, 291, 546

Higuchi, A. E., Sato, A., Tsukagoshi, T., et al. 2017, ApJL, 839, L14

Hildebrand, R. H. 1983, QJRAS, 24, 267

Holland, W. S., Matthews, B. C., Kennedy, G. M., et al. 2017, MNRAS, 470, 3606

Hughes, A. M., Duchêne, G., & Matthews, B. C. 2018, ARA&A, 56, 541

Hughes, A. M., Lieman-Sifry, J., Flaherty, K. M., et al. 2017, ApJ, 839, 86

Hughes, A. M., Wilner, D. J., Kamp, I., & Hogerheijde, M. R. 2008, ApJ, 681, 626

Hung, L.-W., Fitzgerald, M. P., Chen, C. H., et al. 2015, ApJ, 802, 138

Isella, A., Testi, L., Natta, A., et al. 2007, A&A, 469, 213

Kalas, P. 2005, ApJL, 635, L169

Kastner, J. H., Qi, C., Dickson-Vandervelde, D. A., et al. 2018, ApJ, 863, 106

Kennedy, G. M., Marino, S., Matrà, L., et al. 2018, MNRAS, 475, 4924

Kóspál, Á., Moór, A., Juhász, A., et al. 2013, ApJ, 776, 77

Kral, Q., & Latter, H. 2016, MNRAS, 461, 1614

Kral, Q., Marino, S., Wyatt, M. C., Kama, M., & Matra, L. 2019, MNRAS, 489, 3670

Kral, Q., Matrà, L., Wyatt, M. C., & Kennedy, G. M. 2017, MNRAS, 469, 521

Kral, Q., Wyatt, M., Carswell, R. F., et al. 2016, MNRAS, 461, 845

Krivov, A. V. 2010, RAA, 10, 383

Lecavelier des Etangs, A., Vidal-Madjar, A., Roberge, A., et al. 2001, Natur, 412, 706

Lieman-Sifry, J., Hughes, A. M., Carpenter, J. M., et al. 2016, ApJ, 828, 25

Lindegren, L., Hernández, J., Bombrun, A., et al. 2018, A&A, 616, A2

Lynden-Bell, D., & Pringle, J. E. 1974, MNRAS, 168, 603

MacGregor, M. A., Weinberger, A. J., Hughes, A. M., et al. 2018, ApJ, 869, 75

MacGregor, M. A., Wilner, D. J., Chandler, C., et al. 2016, ApJ, 823, 79

Marino, S., Matrà, L., Stark, C., et al. 2016, MNRAS, 460, 2933

Marino, S., Wyatt, M. C., Panić, O., et al. 2017, MNRAS, 465, 2595

Martin-Zaïdi, C., Deleuil, M., Le Bourlot, J., et al. 2008, A&A, 484, 225

Matrà, L., Dent, W. R. F., Wyatt, M. C., et al. 2017a, MNRAS, 464, 1415

Matrà, L., MacGregor, M. A., Kalas, P., et al. 2017b, ApJ, 842, 9

Matrà, L., Marino, S., Kennedy, G. M., et al. 2018a, ApJ, 859, 72

Matrà, L., Öberg, K. I., Wilner, D. J., Olofsson, J., & Bayo, A. 2019, AJ, 157, 117

Matrà, L., Wilner, D. J., Öberg, K. I., et al. 2018b, ApJ, 853, 147

McMullin, J. P., Waters, B., Schiebel, D., Young, W., & Golap, K. 2007, in ASP Conf. Ser. 376, Astronomical Data Analysis Software and Systems XVI, ed. R. A. Shaw, F. Hill, & D. J. Bell (San Francisco, CA: ASP), 127

Meeus, G., Montesinos, B., Mendigutía, I., et al. 2012, A&A, 544, A78

Miotello, A., Bruderer, S., & van Dishoeck, E. F. 2014, A&A, 572, A96

Miyake, K., & Nakagawa, Y. 1993, Icar, 106, 20

Moenchen, M. M., Telesco, C. M., De Buizer, J. M., Packham, C., & Radomski, J. T. 2007, ApJL, 666, L109

Moór, A., Ábrahám, P., Juhász, A., et al. 2011, ApJL, 740, L7

Moór, A., Curé, M., Kóspál, Á., et al. 2017, ApJ, 849, 123

Moór, A., Henning, T., Juhász, A., et al. 2015a, ApJ, 814, 42

Moór, A., Kóspál, Á., Ábrahám, P., et al. 2015b, MNRAS, 447, 577

Morales, F. Y., Bryden, G., Werner, M. W., & Stapelfeldt, K. R. 2016, ApJ, 831, 97

Mumma, M. J., & Charnley, S. B. 2011, ARA&A, 49, 471

Ossenkopf, V., & Henning, T. 1994, A&A, 291, 943

Pawellek, N., Krivov, A. V., Marshall, J. P., et al. 2014, ApJ, 792, 65

Péricaud, J. 2016, PhD thesis, L'université de Bordeaux, Bordeaux

Péricaud, J., di Folco, E., Dutrey, A., et al. 2016, in IAU Symp. 314, Young Stars and Planets Near the Sun, ed. J. H. Kastner, B. Stelzer, & S. A. Metchev (Cambridge: Cambridge Univ. Press), 201

Péricaud, J., Di Folco, E., Dutrey, A., Guilloteau, S., & Piétu, V. 2017, A&A, 600, A62

Raman, A., Lisanti, M., Wilner, D. J., Qi, C., & Hogerheijde, M. 2006, AJ, 131, 2290

Rebollido, I., Eiroa, C., Montesinos, B., et al. 2018, A&A, 614, A3

Redfield, S. 2007, ApJL, 656, L97

Richert, A. J. W., Lyra, W., & Kuchner, M. J. 2018, ApJ, 856, 41

Riviere-Marichalar, P., Barrado, D., Augereau, J.-C., et al. 2012, A&A, 546, L8

Riviere-Marichalar, P., Barrado, D., Montesinos, B., et al. 2014, A&A, 565, A68

Roberge, A., Kamp, I., Montesinos, B., et al. 2013, ApJ, 771, 69

Roberge, A., Welsh, B. Y., Kamp, I., Weinberger, A. J., & Grady, C. A. 2014, ApJL, 796, L11

Rodigas, T. J., Debes, J. H., Hinz, P. M., et al. 2014, ApJ, 783, 21

Rosenfeld, K. A., Andrews, S. M., Hughes, A. M., Wilner, D. J., & Qi, C. 2013, ApJ, 774, 16

Sadakane, K., & Nishida, M. 1986, PASP, 98, 685

Schneider, G., Grady, C. A., Hines, D. C., et al. 2014, AJ, 148, 59

Schneider, G., Silverstone, M. D., & Hines, D. C. 2005, ApJL, 629, L117

Silverstone, M. D. 2000, PhD thesis, Univ. California, Los Angeles

Sylvester, R. J., Skinner, C. J., Barlow, M. J., & Mannings, V. 1996, MNRAS, 279, 915

Takeuchi, T., & Artymowicz, P. 2001, ApJ, 557, 990

van der Plas, G., Ménard, F., Gonzalez, J.-F., et al. 2019, A&A, 624, A33

Visser, R., van Dishoeck, E. F., & Black, J. H. 2009, A&A, 503, 323

Wahhaj, Z., Koerner, D. W., & Sargent, A. I. 2007, ApJ, 661, 36

White, J. A., & Boley, A. C. 2018, ApJ, 859, 103

Wilson, T. L., & Rood, R. 1994, ARA&A, 32, 191

Wyatt, M. C. 2008, ARA&A, 46, 339

Wyatt, M. C. 2018, in Handbook of Exoplanets, ed. H. J. Deeg & J. A. Belmonte (Berlin: Springer), 146

Wyatt, M. C., Panić, O., Kennedy, G. M., & Matrà, L. 2015, Ap&SS, 357, 103
Zuckerman, B. 2019, ApJ, 870, 27

Zuckerman, B., Forveille, T., & Kastner, J. H. 1995, Natur, 373, 494

Zuckerman, B., & Song, I. 2012, ApJ, 758, 77