The Surprisingly Low Carbon Mass in the Debris Disk around HD 32297

No significant asymmetries are visible in either the line or the continuum image. Their appearance is consistent with a ring viewed edge-on. We measure the total C i and continuum emission using the images in Figure 1 by integrating all emission above 3σ. To estimate the error, we collect flux samples by shifting the integration region to parts of the image where no emission is seen. We adopt the standard deviation of the flux samples as our estimate of the statistical error. An additional 10% flux calibration error (Fomalont et al. 2014) is added in quadrature, and this turns out to be the dominant error source. We find a total 492 GHz C i flux of (4.0 ± 0.4) × 10⁻²⁰ W m⁻² and a total 487 GHz (616 μm) continuum flux of (2.2 ± 0.2) × 10⁻²⁸ W m⁻² Hz⁻¹ (22.0 ± 2 mJy). Including previous continuum measurements from Herschel/SPIRE²¹ at 350 and 500 μm and ALMA at 1.3 mm (MacGregor et al. 2018), we derive a submillimeter spectral index of 2.4. This value is comparable to the millimeter spectral indices determined for a sample of debris disks by MacGregor et al. (2016).

Figure 2 shows the PV diagram of the C i emission. We used the PA measured in Appendix A to rotate the data cube in order to align the midplane of the disk with the horizontal direction. The rotated data cube is then integrated from −06 to 07 in the vertical direction z, where z = 0 denotes the midplane.²² Analogous to the PV diagram for CO by MacGregor et al. (2018), we also plot the Keplerian velocity curve, i.e., the extremum radial velocity that can be reached at projected distance x from the star, assuming the same parameters as MacGregor et al. (2018): a stellar mass of 1.7 M_⊙ and an inclination of 836. In the absence of significant line broadening or resolution effects, one expects all emission to lie within these curves—as is observed. Hence, there is no evidence for sub-Keplerian rotation, contrary to the conclusion by MacGregor et al. (2018).

Zoom In Zoom Out Reset image size

We proceed to constrain the spatial distribution and mass of neutral carbon by fitting models to the C i data cube. For a given model, specified by the spatial distribution of the gas and the temperature profile, we first calculate the level population and gas emission in each cell of a three-dimensional grid by assuming local thermodynamic equilibrium (LTE). We verified that LTE is a reasonable approximation by calculating the full radiative transfer in non-LTE using the LIME code (Brinch & Hogerheijde 2010) for a model that fits the data (Figure 4, see below). For that model, the total flux calculated in LTE was only a factor 1.3 higher than the flux from LIME. We use atomic data from the Leiden Atomic and Molecular Database (LAMDA, Schöier et al. 2005). The line emission profile is assumed to be Gaussian with a broadening parameter of 1 km s⁻¹. The emission is red- or blueshifted according to the radial velocity (due to Keplerian rotation) of each grid cell. Then, the model is ray-traced as described in Cataldi et al. (2014) to account for optical depth. The resulting model cube is convolved in the spatial and spectral dimension to match the resolution of the data, and finally multiplied by the primary beam.

Short of a detailed heating–cooling analysis, we simply assume that the gas kinetic temperature behaves similarly to the local blackbody and varies with radial distance r as $T(r)={T}_{100}\sqrt{(100\,\mathrm{au})/r}$, where T₁₀₀ is a normalization constant obtained from the fit. The gas vertical scale height is given by (e.g., Armitage 2009)

$$\begin{eqnarray}&&H(r)=\sqrt{\displaystyle \frac{{{kTr}}^{3}}{\mu {m}_{p}{{GM}}_{* }}},\end{eqnarray} \tag{ 1 }$$

where k is the Boltzmann constant, μ the mean molecular weight, m_p the proton mass, G the gravitational constant, and M_* the stellar mass. We assume μ = 14 (gas mass dominated by carbon and oxygen; Kral et al. 2017), although μ might be higher because the CO mass in HD 32297 is large. The assumptions on the temperature and molecular weight should introduce only minor errors into our results.

We follow Kral et al. (2019) and consider a Gaussian surface density profile for neutral C:

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{{{\rm{C}}}^{0}}(r)={{\rm{\Sigma }}}_{0}\exp \left(-\displaystyle \frac{{\left(r-{r}_{0}\right)}^{2}}{2{\sigma }_{r}^{2}}\right),\end{eqnarray} \tag{ 2 }$$

where Σ₀ is the surface density at r = r₀ and σ_r describes the radial width of the ring-like mass distribution. The number density of neutral carbon is then given by

$$\begin{eqnarray}&&{n}_{{{\rm{C}}}^{0}}(r,z)=\displaystyle \frac{{{\rm{\Sigma }}}_{{{\rm{C}}}^{0}}(r)}{\sqrt{2\pi }{{Hm}}_{{\rm{C}}}}\exp \left(-\displaystyle \frac{{z}^{2}}{2{H}^{2}}\right),\end{eqnarray} \tag{ 3 }$$

where z is the distance from the midplane and m_C the mass of a C atom.

We use a Markov Chain Monte Carlo (MCMC) method implemented by the emcee package (version 2.2.1; Foreman-Mackey et al. 2013) to fit the model to our data (in the image space). The correlated noise in the ALMA data is handled as described by Booth et al. (2017). The following parameters are free to vary: the disk inclination i, the stellar mass M_*, the stellar radial velocity v_*, the ring location r₀, the ring width σ_r, the total mass of neutral carbon ${M}_{{{\rm{C}}}^{0}}$, and the temperature normalization T₁₀₀. In addition, we also fit for astrometric offsets of the disk center, Δx and Δz, which are parallel and perpendicular to the midplane, respectively. For each parameter θ, we assume uninformative priors, either with a location-invariant density (π(θ) = const) or a scale-invariant density (π(θ) ∝ θ⁻¹), over the range given in Table 1. The PA is fixed to 485 as determined in Appendix A. The MCMC is run for 1500 steps with 200 walkers. We discarded the first 200 samples of each walker. Table 1 shows derived parameters (50th percentile) together with error bars (16th and 84th percentiles). Figure 3 shows the posterior distribution of selected parameters. We note that the temperature is not well constrained by our data. This reflects the fact that the fractional population of the transition's upper level peaks at ∼30 K and then stays approximately constant for higher temperatures. The relation between the mass and the temperature is governed by the level population, as expected for an optically thin gas. The derived inclination of $77\buildrel{\circ}\over{.} {9}_{-1.5}^{+1.7}$ is lower than previously derived values from dust observations with ALMA ($83\buildrel{\circ}\over{.} {6}_{-0.4}^{+4.6}$; MacGregor et al. 2018) or from scattered light (∼88°; Boccaletti et al. 2012; Currie et al. 2012). The origin of this discrepancy is unclear, but one possibility is an asymmetry in the disk that our model does not account for. We run an additional MCMC simulation with a fixed inclination of 88° and verified that the results for the other parameters do not change. In Figure 4, we present the fit with the highest posterior probability in our chain. No strong residuals are left, suggesting that our model is a good representation of the true gas distribution. This model has a peak optical depth of 0.2 along the line of sight.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 1.  Fitted Parameters for the Gaussian Ring Model (Section 3.3) and the Double Power-law Model (Appendix B)

Parameter	Valuea		Unit	Priorb	Rangec
	Gaussian	DPLd
i	${77.9}_{-1.5}^{+1.7}$	${78.0}_{-1.5}^{+1.8}$	deg	loc	[65, 90]
r0	${109}_{-4}^{+3}$	${129}_{-11}^{+8}$	au	loc	[30, 200]
${\sigma }_{r}$	${35.4}_{-3.5}^{+3.8}$	⋯	au	sca	[0, 100]
${\beta }_{\mathrm{in}}$	⋯	${1.5}_{-0.6}^{+1.1}$	⋯	loc	[−10, 10]
${\beta }_{\mathrm{out}}$	⋯	$-{7.4}_{-1.6}^{+1.5}$	⋯	loc	[−10, 10]
T100	${117}_{-60}^{+98}$	${125}_{-64}^{+101}$	K	sca	[10, 300]
${M}_{{{\rm{C}}}^{0}}$	${3.5}_{-0.2}^{+0.2}$	${3.6}_{-0.2}^{+0.2}$	10−3 ${M}_{\oplus }$	sca	[0, 1000]
M*	${1.59}_{-0.05}^{+0.05}$	${1.59}_{-0.05}^{+0.05}$	${M}_{\odot }$	sca	[1, 3]
v*	${20.67}_{-0.04}^{+0.04}$	${20.67}_{-0.04}^{+0.04}$	km s−1	loc	[18, 22]
${\rm{\Delta }}x$	$-{0.06}_{-0.01}^{+0.01}$	$-{0.06}_{-0.01}^{+0.01}$	arcsec	loc	[−0.75, 0.75]
${\rm{\Delta }}z$	${0.04}_{-0.01}^{+0.01}$	${0.04}_{-0.01}^{+0.01}$	arcsec	loc	[−0.75, 0.75]

Notes. The parameters β_in and β_out are restricted to the double power-law model.

^a50th percentile with error bars corresponding to the 16th and 84th percentile of the posterior distribution. ^bLocation-invariant (loc) or scale-invariant (sca) prior. ^cParameter range explored by the MCMC fit. ^dDouble power law.

Download table as:  ASCIITypeset image

We note that the calibration uncertainty, not included in the MCMC fit, might introduce some additional uncertainty in the fitted parameters.

In summary, the C⁰ gas distribution can be fit with a ring centered at ∼110 au with an FWHM of ∼80 au. In Appendix B, we also explore fitting the data using power-law models, to accommodate the possibility that the gas has spread into a broad disk. Our results there confirm that the gas is indeed distributed in a ring-like geometry. In Appendix C, we use a simple ionization calculation to confirm that our model is consistent with the C ii emission detected by Donaldson et al. (2013).

Here, we present a simple model to capture a few main processes that affect the time evolution of CO and carbon. This includes UV shielding of CO by carbon, CO self-shielding, and carbon removal (e.g., via recapture by dust grains). The shielding aspect of our model is similar in spirit to that in Kral et al. (2019), Moór et al. (2019), and Marino et al. (2020), but differs (substantially) in the treatment of shielding.

We consider only one radial zone of well-mixed gases and discard the effect of gas removal by viscous spreading (as is justified here; see Section 5.2). The CO mass evolves as

$$\begin{eqnarray}&&\displaystyle \frac{{{dM}}_{\mathrm{CO}}}{{dt}}={P}_{\mathrm{CO}}-{D}_{\mathrm{CO}},\end{eqnarray} \tag{ 5 }$$

where P_CO is the CO production rate by the debris ring and D_CO the CO photodestruction rate. CO is produced in a collisional cascade among the debris.²⁴ Intuitively, we might expect some time dependence of P_CO such as a gradual decrease due to the collisional evolution of the disk (Marino et al. 2020). However, for simplicity, we assume a constant CO production rate. The carbon mass evolves with time as

$$\begin{eqnarray}&&\displaystyle \frac{{{dM}}_{{\rm{C}}}}{{dt}}=\displaystyle \frac{12}{28}{D}_{\mathrm{CO}}-\displaystyle \frac{{M}_{{\rm{C}}}}{{\tau }_{{\rm{C}}}}.\end{eqnarray} \tag{ 6 }$$

The second term on the right-hand side models a carbon sink (e.g., capture onto dust grains). To convert M_C into masses of neutral and ionized carbon, a constant ionization fraction of 0.2 is assumed (Appendix C).

In contrast to Kral et al. (2019), we obtain the photodissociation rate (D_CO) using photon counting. The probability that a given ISRF photon either ionizes C⁰ or dissociates CO is given by

$$\begin{eqnarray}&&{f}_{1}=1-\exp (-{\tau }_{{{\rm{C}}}^{0}}-{\tau }_{\mathrm{CO}})=1-\exp (-{\tau }_{\mathrm{tot}})\end{eqnarray} \tag{ 7 }$$

where the optical depth ${\tau }_{x}(\lambda )={N}_{x}{\sigma }_{x}(\lambda )$ with N_x the vertical column density and σ_x the wavelength-dependent cross section. Under our assumption of perfect mixing, the fraction of these photons dissociating CO is

$$\begin{eqnarray}&&{f}_{2}=\displaystyle \frac{{n}_{\mathrm{CO}}{\sigma }_{\mathrm{CO}}}{{n}_{{{\rm{C}}}^{0}}{\sigma }_{{{\rm{C}}}^{0}}+{n}_{\mathrm{CO}}{\sigma }_{\mathrm{CO}}}=\displaystyle \frac{{\tau }_{\mathrm{CO}}}{{\tau }_{\mathrm{tot}}}\end{eqnarray} \tag{ 8 }$$

where n denotes the number density. Let ϕ_ν be the number of ISRF photons hitting the disk per unit time, per unit frequency, per unit area (taken from Draine 1978). The photodestruction mass rate for CO is

$$\begin{eqnarray}&&{D}_{\mathrm{CO}}={m}_{\mathrm{CO}}\int {\phi }_{\nu }{{Af}}_{1}{f}_{2}\,d\nu \end{eqnarray} \tag{ 9 }$$

where $A=2\pi {r}_{0}{\rm{\Delta }}r$ is the surface of the disk. An alternative derivation of this equation is given in Appendix E.

Equation (9) can be used to derive an upper limit on the CO destruction rate. We assume f₁ = f₂ = 1 for any frequency ν, where σ_CO(ν) > 0, i.e., no photon that could potentially destroy a CO molecule is lost. This leads to ${D}_{\mathrm{CO}}^{\max }\sim 50$ M_⊕ Myr⁻¹ (with r₀ = 110 au and Δr = 80 au). In the case that ${P}_{\mathrm{CO}}\gt {D}_{\mathrm{CO}}^{\max }$, no steady state is possible.

We extend our model to also include the evolution of ¹³CO and C¹⁸O. This is easily achieved by replacing τ_CO in Equation (7) by ${\tau }_{{}^{12}\mathrm{CO}}+{\tau }_{{}^{13}\mathrm{CO}}+{\tau }_{{{\rm{C}}}^{18}{\rm{O}}}$. Then, the masses of the isotopologues are evolved the same way as ¹²CO. The wavelength-dependent cross sections of the different isotopologues only overlap partially. Therefore, the lower abundance isotopologues are more readily photodissociated because their self-shielding is weaker. This effect is called isotopologue selective photodissociation. It increases the ratios of ¹²CO/¹³CO and ¹²CO/C¹⁸O compared to the ISM ratios and is naturally taken into account by our model.

The cross sections for CO photodissociation were provided by A. Heays (2019, private communication) and are based on Visser et al. (2009), while the ionization cross section σ_C is from the database of "Photodissociation and photoionization of astrophysically relevant molecules"²⁵ (hereafter PPARM; Heays et al. 2017); σ_C is essentially constant (∼1.6 × 10⁻¹⁷ cm²) over the wavelength range where CO photodissocation happens (Rollins & Rawlings 2012).

We note that the Kral et al. (2019) model has a different approach to the calculation of the CO destruction rate, and we argue that they underestimate it in the shielded regime. In their model,

$$\begin{eqnarray}&&{D}_{\mathrm{CO}}={M}_{\mathrm{CO}}\times {\chi }_{0}\times S({N}_{\mathrm{CO}})\times \exp (-{N}_{{{\rm{C}}}^{0}}{\sigma }_{{{\rm{C}}}^{0}})\end{eqnarray} \tag{ 10 }$$

where χ₀ is the unshielded ISRF photodissociation rate and S(N_CO) is the reduction on χ₀ due to CO self-shielding under a column density N_CO (Visser et al. 2009). However, Equation (10) is the dissociation rate for a CO molecule that sits behind the full column densities of C and CO. First, such a model ignores the CO dissociation in the shielding layer and therefore overestimates the effect of shielding. In contrast, our model predicts higher CO destruction rates because even if the midplane is strongly self-shielded, ISRF photons are still used to dissociate CO in the higher layers of the disk. Second, the Kral et al. (2019) model also implicitly assumes that C forms a separate layer at the surface of the disk, as expressed by the exponential term in Equation (10). This maximizes the effect of C shielding. In contrast, in our model, C and CO are well mixed. Observations of C i that resolve the vertical structure are needed to decide which prescription is a better representation of a real disk.

Next, we consider carbon condensation onto dust grains. We estimate the total dust cross section by adopting a dust mass of 0.57 M_⊕ (MacGregor et al. 2018), a minimum grain diameter equal to the blowout limit (∼10 μm), a maximum grain diameter of 1 mm, and a power-law size distribution with index −3.5. Adopting T = 100 K and n_C = 370 cm⁻³ (from our Gaussian ring model fit), we find that the mean free time of a C atom before colliding with a dust grain is ∼5000 yr. Grigorieva et al. (2007) also estimated the recondensation of water onto dust grains. By adapting their Equation (17)²⁶ to C, we find that all C is recondensed on a timescale of ∼8000 yr. So, this process might be important. Although collisions between C atoms and grains have nearly unity sticking coefficient (e.g., Hollenbach & Salpeter 1970; Leitch-Devlin & Williams 1985), thermal desorption might remove them again quickly. In a later discussion, we consider this more in detail and also discuss whether the accreted atoms can form fresh CO (Section 5.1).

In the following, we leave the removal timescale τ_C as a free parameter in our model. Figure 5 illustrates three examples of our model for different rates of CO production without and with removal of carbon (taking ${\tau }_{{\rm{C}}}=\infty$ and 5000 yr, respectively). For these models, we fixed r₀ = 110 au and Δr = 80 au. We find that without C removal, once either C⁰ shielding or CO self-shielding is significant, both C and CO rise linearly and the C⁰/CO ratio remains constant. When the CO production rate is low, C shielding is necessary before CO can survive rapid destruction. This explains the initial delay in its mass rise (left panel). Later on, both masses grow linearly, with a fixed fraction of the freshly made CO being turned into C. The ratio of C⁰/CO is above unity. In comparison, in the models with a high CO production rate (middle and right panels), even from the very beginning, CO self-shielding is already significant. We still observe a linear growth in both masses, but the C⁰/CO mass ratio is consistently below unity. On the other hand, removal of C gas allows the system to settle to a steady state in which the removal of C is balanced by its production from CO dissociation. We note that the C⁰/CO ratio decreases with increasing CO production rate, regardless of whether or not C is removed. For HD 32297, the observed C⁰/CO ≪ 1 suggests that a large CO production rate will be required to explain the observations.

Zoom In Zoom Out Reset image size

Here, we apply our above model to the specific case of HD 32297. The planetesimal belt is observed to have a width of ∼40 au if we ignore the "halo component" (MacGregor et al. 2018). Thus, the C gas with its FWHM of ∼80 au is observed to be mildly more spread out. For simplicity, we ignore this difference and fix the width of the gas disk to 80 au.

Our model contains two unknowns: the event time (t₀) at which CO (and presumably dust) production commenced, and the CO production rate (P_CO). We aim to reproduce two observables: the C⁰ mass (3.5 × 10⁻³ M_⊕) and the mass of C¹⁸O (1.4 × 10⁻⁴ M_⊕, Moór et al. 2019). We use the C¹⁸O mass because the main component of the CO gas, ¹²CO, is optically thick and its mass cannot be determined accurately. We also ignore ¹³CO as it may be optically thick as well.

4.2.1. Steady-state Solution

4.2.2. Full Calculation

We now do not assume steady state and instead evolve a grid of models with different event times (t₀) and CO production rates (P_CO) with the goal of reproducing the observed masses of C⁰ and C¹⁸O. We fix the isotopologue ratios in the gas-producing solids (i.e., the CO production rate ratios) to ISM ratios.

Our model results are presented in Figure 6. In the case of no C removal (${\tau }_{{\rm{C}}}=\infty$, left panel), we find that only the combination of an extremely high CO production rate and a recent event can explain the data: P_CO ≈ 100 M_⊕ Myr⁻¹ and t₀ ≈ 1000 yr. As can be seen from the figure, for lower CO production rates, that model exceeds the observed C mass long before it has built up the observed C¹⁸O mass.

Zoom In Zoom Out Reset image size

Figure 6 (right) shows a scenario were τ_C is adjusted such that the observed gas masses can be explained with a steady-state solution. We find τ_C ∼ 700 yr is necessary, i.e., C gas needs to be removed efficiently in order to reproduce the small C mass we observe: for larger values of τ_C, the results are qualitatively similar to the ${\tau }_{{\rm{C}}}=\infty$ case, while for smaller values, no solution is found. Even in this scenario, the required CO production rate remains high at ∼15 M_⊕ Myr⁻¹ (as mentioned in Section 4, based on the dust mass-loss rate, we would expect a CO production rate of ∼0.5 M_⊕ Myr⁻¹). Steady state is reached after a few times 10⁵ yr. The steady-state CO mass is 1.3 M_⊕. These results are consistent with the steady-state solution described in Section 4.2.1.

In Figure 7, we show the two models that can formally explain the observed gas masses, with parameters as identified from Figure 6. As mentioned above, the first model (no removal of carbon) requires a very recent event and an extremely large CO production rate. By the time the model fits the observed C⁰ and C¹⁸O masses, the total CO mass released is ∼0.1 M_⊕. This model can thus be interpreted as the recent destruction of a ∼1 M_⊕ body with a CO fraction of 10% releasing all of its CO within 1000 yr. In the following, we focus on discussing the second model. In this latter case, the model does not give us a unique answer regarding the event time, other than that it must be greater than ∼10⁵ yr. For comparison, the system lifetime is larger than 15 Myr. The CO production rate of 15 M_⊕ Myr⁻¹ is in contrast with the result from Moór et al. (2019), where they found that P_CO ≈ 3. 6 × 10⁻² M_⊕ Myr⁻¹. However, the Moór et al. (2019) model overpredicts the C⁰ mass by a factor of 10 compared to our data, despite underestimating the CO breakdown rate due to the different treatment of shielding (Section 4.1).

Zoom In Zoom Out Reset image size

Here, we consider the fate of carbon atoms that collide with dust grains. We also discuss the possibility that accreted C atoms will react with O to re-form CO. Although O i 63 μm emission was not detected by Herschel (Donaldson et al. 2013), we expect O to be at least as abundant as C, thanks to photodissociation of CO, CO₂, and H₂O, all of which should be present in the cometary parent bodies.

We first consider if atoms can remain on the grain or are quickly rereleased to the gas phase by thermal desorption. The rate of thermal desorption R_td depends on the binding energy E_i (Hollenbach & Salpeter 1970; Tielens & Hagen 1982):

$$\begin{eqnarray}&&{R}_{\mathrm{td}}\approx {\nu }_{i}{e}^{-{E}_{d}/{{kT}}_{\mathrm{gr}}}\end{eqnarray} \tag{ 11 }$$

where ν_i ∼ 10¹³ s⁻¹ (Tielens & Hagen 1982) is the characteristic oscillation frequency for the atom in the lattice potential. The value of the binding energy depends on the type of interaction between the adsorbate and the surface (e.g., Tielens 2005). Binding by van der Waals forces (physiosorption) is relatively weak (0.01–0.2 eV; Tielens 2005), while valence bonds (chemisorption) are stronger (∼1 eV). Whether adsorbed atoms are physiosorbed or chemisorbed depends on the grain properties (e.g., chemical composition, temperature; Y. Aikawa 2020, private communication). For instance, let us assume a grain temperature of 50 K. If C is physiosorbed with a binding energy of 0.1 eV, it would be thermally desorbed in a fraction of a second. In such a case, dust grains cannot act as a C sink. On the other hand, chemisorbed C can be accumulated on the grains. More studies (e.g., quantum chemical computations such as in Shimonishi et al. 2018) are needed to determine the binding energy of C (and O) on dust grains.

Next, we consider the possibility that accreted C and O atoms re-form CO on the dust grains. Detailed modeling to determine if and under which conditions C+O→CO can proceed on the dust grains²⁷ is beyond the scope of this paper. However, if the reaction proceeds, it has interesting consequences. CO is then continuously recycled: it is released from grains, photodissociated by the ISRF, C and O reaccreted onto grains, CO re-formed rapidly, and rereleased into space. In this picture, as long as a certain amount of CO is injected as gas at some point, there is no need for fresh CO production. The CO production rate that we infer would simply be a result of CO re-formation. Consider a C mass of 3 × 10⁻³ M_⊕ and a removal time of τ_C = 700 yr. This yields the rate of CO re-formation, ∼10 M_⊕ Myr⁻¹, as we obtain above. It is not necessary for the CO production rate to be limited by the dust grind-down rate. CO re-formation would also explain why the carbon gas is concentrated into a narrow radial belt—all gas components have been released so recently that they have not had time to spread viscously.

We experimented with a modified version of our model with an additional CO source term equal to $28/12{M}_{{\rm{C}}}/{\tau }_{{\rm{C}}}$, i.e., each removed C atom is immediately converted to a CO molecule. Then, the observed C and CO masses can indeed be reproduced with much smaller CO production rates (P_CO ≲ 1 M_⊕ Myr⁻¹) compared to our standard model.

While consideration of dust capture yields a carbon removal time of τ_C ∼ 5000 yr (Section 4.1), our model requires a much shorter τ_c ∼ 700 yr, in order to explain the low C⁰ mass we infer from the ALMA observations. We consider a number of alternative carbon removal mechanisms, but do not find a satisfactory solution to this problem at the moment.

• 1.  
  Gas accretion into the star appears ineffective in HD 32297, as the gas ring appears narrow, not radially spread into an accretion disk. Indeed, our modeling of the C i emission suggests that the gas is distributed in a ring that is only mildly broader than the dust ring.
• 2.  
  Can one or several planets inward of the dust belt prevent the formation of an accretion disk by accreting gas themselves? So far, no giant planets have been found in direct imaging searches around HD 32297 (Boccaletti et al. 2012; Rodigas et al. 2014; Bhowmik et al. 2019). In the most recent study, Bhowmik et al. (2019) exclude the presence of planets with 2.5–6 Jupiter masses at projected distances beyond ∼80 au. However, because of the edge-on geometry of the disk, planets might be hiding close to the line of sight toward the star. Furthermore, Marino et al. (2020) found that smaller, currently undetectable Neptune-class planets can be sufficient to accrete gas.
• 3.  
  Marino et al. (2020) recently suggested that C might become undetectable as it viscously spreads outwards to regions at larger radii with lower surface density. More detailed modeling is needed to test this proposition. Naively, one might expect this to be less of a problem for an edge-on disk like HD 32297 where the column density along the line of sight is high compared to a face-on disk. This scenario would also imply efficient spreading and thus would again require a mechanism to prevent the formation of an inward accretion disk (see the previous point).
• 4.  
  What about carbon removal by gas-phase chemistry? Consulting the chemical network and reaction rates relevant for protoplanetary disks, as compiled by Glover et al. (2010), we find that the most probable pathway for carbon removal is the following reaction:²⁸

  $$\begin{eqnarray}&&\mathrm{OH}+{\rm{C}}\to \mathrm{CO}+{\rm{H}},\end{eqnarray} \tag{ 12 }$$

  with a reaction rate of $k={10}^{-10}$ cm³ s⁻¹ (similarly, OH + C⁺ → CO⁺ + H has k = 7.7 × 10⁻¹⁰ cm³ s⁻¹). Thus, to match our required C removal timescale of τ_C = 700 yr, we would require ${n}_{\mathrm{OH}}={({\tau }_{{\rm{C}}}k)}^{-1}=0.5$ cm⁻³. The OH molecules are likely the results of water photodissociation, and they are highly vulnerable to UV photodissociation. So, for this channel to be effective in removing carbon gas, one requires a copious amount of water production. We encourage more detailed modeling to assess how chemistry influences the evolution of debris disk gas.
• 5.  
  Can carbon gas spontaneously condense into solid form, without the assistance of dust grains? The low-temperature phase diagram of carbon (see, e.g., Jaworski et al. 2016) shows that the vapor pressure for condensation of gaseous atomic C into solid state at the relevant temperature range is orders of magnitude below the gas pressure in the disk. However, typical condensation requires the presence of seed nuclei. These nuclei ought to be scarce as they are quickly removed by the radiation pressure from the A star (Weingartner & Draine 2001), although sufficient amounts of gas could lengthen their residence time.

Another possible solution might be to assume a CO production that is enhanced in isotopologues. For instance, we experimented with a model where the production rate of ¹³CO and C¹⁸O, relative to ¹²CO, is enhanced by a factor of 10 from the ISM ratio. At steady state, we require a CO production rate of 5 M_⊕ Myr⁻¹, and a carbon removal time of τ_c ∼ 2000 yr, more similar to the 5000 yr inferred in Section 4.1.

The conclusions we draw from our modeling rely heavily on our main observational result: the low mass of neutral carbon. How robust is this mass measurement? Could it be that the small error bars (see Table 1) derived from the MCMC modeling are underestimating the real uncertainty? One scenario that could increase the C⁰ mass estimate is a low kinetic temperature. This could lead to a high optical depth, allowing for an arbitrarily high C⁰ mass, in principle. We discuss this possibility in Appendix D. Overall, we believe that our mass error bars are reasonable.

Regarding UV shielding, we can exclude the possibility of another shielding agent. Other than carbon and CO, there are no promising molecules or atoms in the PPARM database that have the required high cross section at the relevant wavelengths to provide shielding.

Next, the strength of the ISRF may be somewhat uncertain. Very few direct observations exist (e.g., Henry 2002), and there may also be some spatial variations caused by proximity to massive stars (e.g., van Dishoeck 1994). A weaker ISRF than the canonical value assumed here would reduce the CO destruction rate, and correspondingly, the inferred CO production rate.

Finally, our model assumes that CO and C are well mixed. When instead assuming that all C is concentrated at the surface of the disk (this configuration maximizes C shielding), the derived CO production rates and C removal timescales only change by a factor of ∼2. This reflects the modest shielding by the small measured C mass. For larger C masses, the difference can be more dramatic.

When was the dusty disk around HD 32297 produced? Is it as old as the star, or is it the result of a recent catastrophic event? What is the total mass of the underlying planetesimal belt? Our model without C removal requires a catastrophic event 1000 yr ago, i.e., the CO gas was released within an orbital timescale. In such a scenario, the symmetric appearance of the disk in CO and C i is challenging to explain. On the other hand, our steady-state model cannot give a unique answer to the event time, so we look to other evidence for guidance.

If the CO gas has been steadily produced over a long timescale, it and its descendants (carbon and oxygen gas) ought to have viscously spread out and eventually accreted onto the star. For a gas belt of width Δr and an α parameter of 0.01, as likely appropriate for gaseous debris disks (Kral & Latter 2016), the time required to double the radial width is of order

$$\begin{eqnarray}&&t\sim \displaystyle \frac{{\left({\rm{\Delta }}r\right)}^{2}}{\alpha {c}_{s}H}\sim {10}^{6}\,\mathrm{yr}{\left(\displaystyle \frac{{\rm{\Delta }}r}{40\mathrm{au}}\right)}^{2}.\end{eqnarray} \tag{ 13 }$$

However, our modeling of the C i emission (Section 3.3) suggests that the gas is distributed in a ring that is only mildly broader than the dust ring, indicating that there has been insignificant radial spreading around HD 32297 since gas production started. Could this imply that the event time is shorter than the system lifetime (≥15 Myr)? The answer is no. Any newly outgassed CO is necessarily concentrated near the grains' orbits. It is photodissociated within a time short compared to the above spreading time (mean lifetime in our model is 0.1 Myr). And its descendant, the C gas, is rapidly removed (τ_C ∼ 700 yr). So no radial spreading is allowed in the steady-state model.

The high CO production rate that we obtain (∼15 M_⊕ Myr⁻¹) can in principle place a constraint on the event time, as it is not sustainable over a long period of time. However, this is not true if CO recycling occurs. In the latter case, even if gas production started right after the dispersal of the protoplanetary disk, the minimum CO required is still only of order M_⊕ (the current CO gas mass), roughly the amount of CO contained by a comet belt of some 10 M_⊕. This total mass is consistent with the estimates based on the gradual luminosity evolution of debris disks (e.g., Shannon & Wu 2011).

Previously, Cataldi et al. (2018) concluded that the debris disk around β Pic had to have been produced recently. In that system, the ¹²CO emission is optically thin, and the CO mass can be directly measured (Matrà et al. 2017a). It is at least three orders of magnitude lower than that in HD 32297, while the carbon masses are comparable. Moreover, both CO and carbon are distributed in a highly asymmetric, radially confined ring (Cataldi et al. 2018). The low C and CO masses imply little shielding, and the CO lifetime is as short as ∼50 yr, or a CO destruction rate of 0.6 M_⊕ Myr⁻¹. Assuming that the low amount of C is the integrated result of this destruction rate, Cataldi et al. (2018) concluded that carbon can be accumulated in a time as short as ≲10⁴ yr. However, this picture of a short event time would drastically be changed if C is removed. In such a case, the low carbon mass no longer places a constraint on the event time directly.

Besides β Pic and HD 32297, there are currently two other debris disks with resolved CO and C i observations: HD 131835 and 49 Ceti. For HD 131835, Kral et al. (2019) derived r₀ = 90 au, Δr = 70 au, and a C⁰ mass between 2.7 × 10⁻³ M_⊕ and 1.2 × 10⁻² M_⊕. The CO mass is estimated at 0.04 M_⊕ (Moór et al. 2017; Kral et al. 2019). The spatial distribution of the gas is uncertain, and there remains the possibility that an accretion disk has formed (Kral et al. 2019). For 49 Ceti, the dust belt is concentrated around r₀ ∼ 100 au, with a width Δr ∼ 80 au (the double power-law model of Hughes et al. 2017), while the carbon mass is ${M}_{{{\rm{C}}}^{0}}=3\,\times {10}^{-3}$ M_⊕ (assuming optically thin emission; A. Higuchi 2019, private communication) and the CO mass (if ISM ratios apply) (1.11 ± 0.13) × 10⁻² M_⊕ (Moór et al. 2019). Hughes et al. (2017) found that the CO observations are well described by a model where the surface density increases with radius between ∼20 and ∼220 au, suggesting that an accretion disk is unlikely. Both disks might be qualitatively explained by a scenario similar to HD 32297 (Figure 7 right).

Recently, Marino et al. (2020) presented a population synthesis study by coupling the secondary gas production model by Kral et al. (2019) with a dust collisional evolution model (i.e., the dust grind-down rate, and thus CO production, decays with time). While they find agreement with measured CO masses, their model tends to overpredict the masses of neutral C, although they get better agreement by assuming that C beyond the belt location is not detectable in observations. This result seems analogous to our finding that the small C mass around HD 32297 is difficult to reproduce with a CO production rate derived from the dust mass-loss rate.