Primordial or Secondary? Testing Models of Debris Disk Gas with ALMA^*

The goal of the present study is to get an overview of the C and CO gas in debris disks. Therefore, our sample includes any disks listed in the "Debris" category of the Catalog of Circumstellar Disks ¹⁴ that had observations of the C i ³P₁–³P₀ line and at least one CO line in the ALMA archive as of 2021 September. We used the ALminer software to query the ALMA archive. We identified 14 disks meeting our criteria. For 10 of these disks, the ALMA observations of CI have not been published previously. Table 1 gives an overview of our sample and lists the stellar and disk parameters we used.

Table 1. Stellar and Disk Parameters Used in This Study

Star	SpT	d	vsys	M*	Teff	$\mathrm{log}g$	[M/H]	L*	Age	f	i	PA	rin	rout	r Reference
		(pc)	(km s−1)	(M⊙)	(K)	(log(cm s−2))	(⋯)	(L⊙)	(Myr)	(×10−3)	(deg)	(deg)	(au)	(au)
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)	(10)	(11)	(12)	(13)	(14)	(15)	(16)
49 Ceti	A1V (1)	57.1	11.92 (2)	2.1 (2)	9120 (1)	4.32 (1)	0 (1)	17.2 (3)	40 (4)	1.1 (5)	79.5 (2)	−72 (2)	19.3 a	153.8 a	CO (2)
β Pictoris	A6V (6)	19.4	20.5 (7)	1.75 (8)	8052 (6)	4.15 (6)	0.05 (6)	8.1 (9)	23 (10)	2.1 (9)	86.6 (11)		50	160	CO (12)
HD 21997	A3IV/V (1)	69.6	17.92 (13)	1.8 (13)	8520 (1)	4.27 (1)	0 (1)	9.29 (9) a	42 (14)	0.6 (9)	32.6 (13)	202.6 (13)	25.1 a	133.5 a	CO (13)
HD 32297	A6V (15)	132.8	20.67 (16)	1.59 (16)	7980 (1)	3.77 (1)	−0.5 (1)	8.46 (3) a	<30 (17)	4.4 (18)	77.9 (16)		67.6 a	151.2 a	C i (16)
HD 48370	G8V (19)	36.1	23.6 (20,21)	0.94 (22)	5600 (22)	4.5 (22)	−0.03 (23)	0.76 (22) a	42 (14,22)	0.6 (22)	69.8 (24)	67.3 (24)	73.6	145.7	Millimeter dust (24)
HD 61005	G8Vk (6)	36.5	22.5 (25)	0.9 (9)	5500 (25)	4.5 (25)	0.01 (25)	0.75(9) a	40 (9)	2.3 (9)	85.6 (26)	70.3 (26)	41.9	78	Millimeter dust (26)
HD 95086	A8III (27)	86.4	17 (28)	1.6 (29)	7750 (28)	4 (28)	−0.25 (28)	6.4 (30) a	13 (31)	1.7 (5)	31 (32)	278 (32)	110.4 a	337.3 a	Millimeter dust (32)
HD 110058	A6/7V (33)	130.1	12.26 (33)	1.84 (33)	7839 (34)	4 (34)	−0.5 (1)	8.65 (30) a	15 (35,1)	1.4 (18)	85.5 (33)	335.1 (33)	7.4	21.2	CO (33)
HD 121191	A5IV/V (1)	132.1	10.1 (36)	1.6 (36)	7970 (1)	4.38 (1)	0 (1)	7.75 (30) a	16 (35,1)	4.7 (37)	28 (36)	25 (36)	14	30	CO (36)
HD 121617	A1V (1)	116.9	7.8 (38)	1.9 (9)	9160 (1)	4.13 (1)	0 (1)	14.13 (30) a	16 (35,1)	4.8 (39)	37 (30)	43 (30)	49.2 a	101.2 a	Millimeter dust (30)
HD 131835	A2IV (1)	133.7	3.45 (40)	1.77 (41)	8610 (1)	4.24 (1)	0 (1)	11 (40)	16 (35,1)	3 (41)	68.6 (40)	57.7 (40)	22.8 a	174.9 a	Millimeter dust (42)
HD 146897	F2/F3V (1)	131.5	−3.1 (1)	1.5 (9)	6700 (1)	4.3 (1)	0 (1)	3.25 (9) a	10 (9)	8.2 (9)	84.6 (43)	114.6 (44)	75.9 a	97.3 a	Millimeter dust (42)
HD 181327	F5/F6V (1)	48.2	0.1 (45)	1.36 (46)	6360 (1)	4.09 (1)	−0.05 (1)	2.87 (47)	23 (10,1)	4.1 (47)	30 (45)	98.8 (45)	69.2 a	90.8 a	Millimeter dust (45)
HR 4796	A0V (48)	71.9	7.5 (1)	2.18 (49)	10060 (1)	4.44 (1)	−0.5 (1)	22.83 (50) a	8 (18)	3 (51)	76.6 (52)	26.7 (52)	70.5 a	84.4 a	Millimeter dust (52)

Notes. References are given in parenthesis. (1) star name; (2) spectral type; (3) distance from Gaia Collaboration et al. (2018), except for β Pictoris where we use the more precise value by van Leeuwen (2007); (4) systemic velocity in the barycentric frame; (5) stellar mass; (6) stellar effective temperature; (7) surface gravity; (8) stellar metallicity; (9) stellar luminosity; (10) system age; (11) disk fractional luminosity; (12) disk inclination; (13) disk position angle (missing values indicate that it was not necessary to adopt a PA because no flux measurements were carried out for that target); (14) adopted inner edge of the disk; (15) adopted outer edge of the disk; (16) type of observations that were used to define r_in and r_out.

^a Indicates scaling to updated GAIA distance measures.

References. (1) Rebollido et al. (2020), (2) Hughes et al. (2017), (3) Moór et al. (2019), (4) Zuckerman & Song (2012), (5) Moór et al. (2015a), (6) Gray et al. (2006), (7) Brandeker (2011), (8) Crifo et al. (1997), (9) Matrà et al. (2018a), (10) Mamajek & Bell (2014), (11) Matrà et al. (2019b), (12) Dent et al. (2014), (13) Kóspál et al. (2013), (14) Bell et al. (2015), (15) Rodigas et al. (2014), (16) Cataldi et al. (2020), (17) Kalas (2005), (18) Chen et al. (2014), (19) Torres et al. (2008), (20) Gaia Collaboration et al. (2016), (21) Gaia Collaboration et al. (2018), (22) Moór et al. (2016), (23) Gáspár et al. (2016), (24) A. Moor (2023, in preparation), (25) Desidera et al. (2011), (26) MacGregor et al. (2018), (27) Houk & Cowley (1975), (28) Moór et al. (2013a), (29) Booth et al. (2019), (30) Moór et al. (2017), (31) Booth et al. (2021), (32) Su et al. (2017), (33) Hales et al. (2022), (34) Saffe et al. (2021), (35) Pecaut & Mamajek (2016), (36) Kral et al. (2020b), (37) Vican et al. (2016), (38) Rebollido et al. (2018), (39) Moór et al. (2011), (40) Hales et al. (2019), (41) Moór et al. (2015b), (42) Lieman-Sifry et al. (2016), (43) Engler et al. (2017), (44) Goebel et al. (2018), (45) Marino et al. (2016), (46) Lebreton et al. (2012), (47) Pawellek et al. (2021), (48) Houk (1982), (49) Gerbaldi et al. (1999), (50) Kral et al. (2017), (51) Meeus et al. (2012), (52) Kennedy et al. (2018).

Download table as:  ASCIITypeset image

For each disk, we define an inner radius r_in and an outer radius r_out based on a dust or gas model presented in the literature, as indicated in the last row of Table 1. If the literature model is a Gaussian ring (HD 121191, HD 121617, HD 32297, HD 181327, HD 48370), we set r_in = r₀ − FWHM/2 and r_out = r₀ + FWHM/2, where r₀ is the center of the Gaussian. For the disk models of 49 Ceti by Hughes et al. (2017) and HD 110058 by Hales et al. (2022), no outer radius is defined. Thus, we adopt the radius where the model surface density profile is at 50% of the peak value as r_out. The same procedure is employed for HD 61005. While the power-law disk model by MacGregor et al. (2018) for HD 61005 does specify an outer cutoff radius (at 188 au), the power-law index is very steep, meaning that our adopted r_out = 78 au better represents the location of the bulk disk material. In all other cases, r_in and r_out are adopted directly from the literature. Where necessary, r_in and r_out (as well as the stellar luminosity) were scaled to take into account updated distance measurements.

Besides CO and C i observations by ALMA, we also used C ii ²P_3/2–²P_1/2 observations from Herschel/PACS where available. Table 2 provides an overview of the emission line data we used for each disk. Note the following data sets available in the ALMA archive that we ignored. For C i toward HD 121191, we found that when combining the ALMA data of program 2019.1.01175.S with the ALMA Compact Array (ACA) data of program 2018.1.00633.S, the signal-to-noise ratio (S/N) of the disk-integrated flux worsens. Thus, we ignored the noisier ACA data. For ¹²CO 3–2 toward HD 181327, we ignored data from program 2013.1.00025.S because the total integration time is only 60 s. For ¹²CO 2–1 toward HD 61005, we follow Olofsson et al. (2016) and only use the data acquired on 2014 March 20 for which the amount of water vapor in the atmosphere was smallest.

Table 2 also lists the disk-integrated line fluxes we used. For emission lines without published fluxes, we derived disk-integrated fluxes as described in Section 3. Otherwise, we used fluxes from the literature, with the following exceptions. For HD 21997, fluxes of ¹²CO 2–1 and 3–2, ¹³CO 2–1 and 3–2, and C¹⁸O 2–1 from program 2011.0.00780.S were published by Kóspál et al. (2013). However, these data were recently reanalyzed ¹⁵ and supplemented by new observations of the 3–2 transition of ¹³CO and C¹⁸O (program 2017.1.01575.S). We use naturally weighted data cubes derived from the reanalyzed 2011.0.00780.S data and the 2017.1.01575.S data. These will be presented in more detail in a forthcoming publication (L. Matrà et al. 2023, in preparation). For ¹³CO 3–2, which was observed by both programs, we measure the disk-integrated fluxes from both data cubes (as described in Section 3.1.2) and adopt the weighted average as follows. The weights are ${w}_{i}={\sigma }_{i}^{-2}$, with σ_i the errors on the fluxes. The error on the weighted mean is then calculated as

$$\begin{eqnarray}&&\sigma =\sqrt{{\left(\displaystyle \sum _{i=1}^{2}{\sigma }_{i}^{-2}\right)}^{-1}}.\end{eqnarray} \tag{ 1 }$$

Our results are consistent with those of Kóspál et al. (2013), except for C¹⁸O 2–1 where we find a significantly higher flux.

For HD 95086, Booth et al. (2019) presented fluxes of ¹²CO 1–0, 2–1, and 3–2, where the 2–1 transition was tentatively detected (see also Zapata et al. 2018). However, Booth et al. (2021) recently proposed that HD 95086 belongs to the Carina young association, implying a different systemic velocity than assumed by Booth et al. (2019). Thus, the ¹²CO 2–1 signal is likely spurious and we re-derived the fluxes using the updated systemic velocity.

For HD 110058, Hales et al. (2022) published CO fluxes after we already had conducted our flux measurements. Furthermore, for ¹²CO 2–1, they used data from program 2018.1.00500.S only, while we included additional data from program 2012.1.00688.S. We stick with the CO fluxes derived in this work, which are consistent with the results by Hales et al. (2022).

Finally, we also re-derived the C ii ²P_3/2–²P_1/2 flux toward β Pictoris (β Pic), as described in Section 3.2. The same data were already analyzed by Brandeker et al. (2016).

We note that after 2021 September, more disks that meet our selection criteria as well as more line observations for the already selected disks have become available in the ALMA archive. Analyzing an extended sample including these new data will be the subject of a future publication.

Here we describe the calibration and imaging of the ALMA (and ACA) data that were used to derive disk-integrated fluxes. We start with the pipeline-calibrated visibility data downloaded from the ALMA archive that we process and image as described below. The one exception is the CO data of HD 21997 where we use image cubes that will be presented in more detail in a forthcoming publication (L. Matrà et al. 2023, in preparation, see Section 2.1).

2.2.1. Calibration

2.2.2. Imaging

When available, we use published C ii ²P_3/2–²P_1/2 fluxes measured with Herschel (Table 2). However, for HD 21997 the Herschel data have not been published yet. We also decided to reanalyze the data for β Pic. For both targets, we used observations by Herschel/PACS (Pilbratt et al. 2010; Poglitsch et al. 2010). The PACS spectrometer array consists of 5 × 5 spatial pixels (spaxels) covering a 47'' × 47'' field of view. We directly use the data downloaded from the Herschel Science Archive.

HD 21997 was observed in the "Mapping" observing mode, meaning that a raster map was constructed using several pointings. We use the spatially resampled and mosaicked projected cube (HPS3DPR). The beam FWHM is approximately 115; ¹⁷ thus, we do not expect the target to be resolved. We extract a spectrum by spatially integrating the data cube over a circular aperture with a diameter equal to the beam FWHM, centered on the proper motion-corrected position of the star. We only consider the region between 157.1 and 158.5 μm because the noise increases toward the edges of the spectral window. The FWHM of an unresolved line is Δλ = 0.126 μm (239 km s⁻¹), ¹⁸ implying that the C ii emission line will not be resolved. We subtract the continuum by fitting a linear polynomial to the spectrum, excluding the region λ₀ ± (1.5Δλ) where λ₀ is the rest wavelength of the line.

β Pic was observed in the "Pointed" observing mode, meaning that only a single pointing was carried out. Thus, instead of a raster map as in the case of HD 21997, we only get a single spectrum for each of the 25 spaxels of PACS. We use the spectrally rebinned cube (HPS3DRR). Given the beam size of 115, the spaxel size of 94 × 94 and a disk diameter of roughly 10'', we can expect that emission is restricted to the central nine spaxels (see also Figure 3 of Brandeker et al. 2016). We extract a single spectrum by summing over the central nine spaxels. We then subtract the continuum in the same way as for HD 21997.

3.1.1. Moment 0 Maps

To measure CO and C i fluxes, we first produce moment 0 maps. For each emission line, we first define a velocity range v_sys ± Δv where v_sys is the systemic velocity and Δv is the projected Keplerian velocity at the inner edge of the disk r_in (see Table 1 for the adopted values of r_in). Thus,

$$\begin{eqnarray}&&{\rm{\Delta }}v=\sqrt{\displaystyle \frac{{{GM}}_{* }}{{r}_{\mathrm{in}}}}\sin i,\end{eqnarray} \tag{ 2 }$$

with the stellar mass M_∗ and the disk inclination i as listed in Table 1. For low inclinations, this tends to zero, but our sample does not include very low inclination disks and thus the smallest Δv is 1.9 km s⁻¹. The adopted velocity range always spans several channels. We then extract a spectrum by spatially integrating the image cube over an elliptical region defined by the disk inclination, position angle, and disk outer edge r_out. We use this spectrum to verify by eye that our velocity range encompasses all emission (obviously, this is only possible in the cases where emission is detected). This is indeed the case, with one exception: for HD 121617, we needed to manually increase the velocity range from 6.7 to 10 km s⁻¹. Since for this disk, r_in is based on continuum observations, this could indicate that the gas extends further inwards than the dust. We then create a moment 0 map using the bettermoments package by integrating the image cube over the adopted velocity range. The moment 0 maps for the disks with new C i detections are shown in Figures 1 (HD 21997), 2 (HD 121191), and 3 (HD 121617). For each disk, we show all moment 0 maps and continuum maps we created. The complete set of moment 0 and continuum maps we produced and used for flux measurements is available in the online journal.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

View figure set (12 panels)

Tarball of figure set images

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

3.1.2. Measurement of Disk-integrated Fluxes

We now measure disk-integrated CO and C i fluxes from the moment 0 maps by using apertures. We consider two kinds of apertures: (1) apertures constructed using the known geometrical parameters of the disk (geometrical aperture), and (2) empirical apertures based on the observed emission in the moment 0 map (empirical aperture). The latter can only be applied if the S/N is sufficiently high.

Let us first consider the geometrical apertures. We consider a disk model with inner and outer cutoffs given by r_in and r_out (Table 1) and constant scale height

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

with k the Boltzmann constant, T the temperature (fixed to 50 K for simplicity), μ the mean molecular weight of the gas (for simplicity fixed to 14, Kral et al. 2017), and m_p the proton mass. We define the mean radius r_mean as

$$\begin{eqnarray}&&{r}_{\mathrm{mean}}^{2}=\displaystyle \frac{{r}_{\mathrm{in}}^{2}+{r}_{\mathrm{out}}^{2}}{2},\end{eqnarray} \tag{ 4 }$$

such that the areas (and thus gas masses, given the constant scale height) inside and outside of r_mean are equal. The disk is placed in a 3D grid taking into account its inclination and position angle. We set any point within our grid to 1 if it satisfies r_in ≤ r ≤ r_out (with r the cylindrical coordinate) and is within ±2H of the disk midplane. All other points are set to 0. We then sum the grid along the line of sight. In the resulting 2D map, any point larger than 0 is set to 1. Finally, we convolve the map with the beam of the observations and normalize it to a peak of 1. Our aperture then consists of all points larger than 0.2.

In some cases (most significantly for the CO lines of HD 21997, see Figure 1), there was a visible offset between the expected and observed disk center. We then shifted the geometrical aperture by eye so that all flux is included.

We then measure the noise ξ in a region outside of the disk. Note that we use the non-primary beam-corrected images to assure that the noise is uniform throughout the image. This is justified because the disks are small compared to the primary beam. Indeed, considering the normalized primary beam response (peak 1), no aperture extends further out than the 0.73 level of the primary beam. In 95% of the cases, the aperture stays within the 0.86 level of the primary beam. We also verified that disk-integrated fluxes measured from primary beam-corrected images are not significantly different. To this end, we measured fluxes from the primary beam-corrected images using the same apertures as for the uncorrected images. The difference was always smaller than 0.4σ, and for the vast majority of the images, smaller than 0.1σ, where σ denotes the error on the disk-integrated flux.

Next, we measure the peak S/N within the geometrical aperture. If the peak S/N within the geometrical aperture is below 8, the flux obtained by integrating over this aperture is used as our final flux measurement. If the S/N is above 8, we measured the flux from an empirical aperture instead. The empirical aperture is constructed as follows. We consider all pixels of the image with a flux F > 2.5ξ, where ξ is the noise in the moment 0 map. This results in a number of disjoint regions. The region corresponding to the disk can be easily identified as the one containing the highest flux. This region is adopted as our empirical aperture. We found that a cutoff at 2.5ξ provides a good compromise in including most of the real emission and excluding noise.

The moment 0 and continuum maps of Figures 1–3 show the apertures we used to measure the fluxes. The complete set of moment 0 and continuum maps is available in the online journal.

To estimate the error σ on the disk-integrated flux, we place the (geometrical or empirical) aperture at various positions outside the disk (without overlap) and take the standard deviation of the collected flux samples. However, sometimes it is not possible to collect a sufficient number of non-overlapping flux samples because the aperture is too big and/or the image is too small. Thus, if less than 20 flux samples can be collected, we instead calculate σ analytically as

$$\begin{eqnarray}&&\sigma =\xi {{\rm{\Omega }}}_{p}\sqrt{{N}_{p}}f=\xi {{\rm{\Omega }}}_{{\rm{b}}}\sqrt{{N}_{b}},\end{eqnarray} \tag{ 5 }$$

where Ω_p and Ω_b are the solid angles (in units of steradian) of a pixel and the beam, respectively, N_p is the number of pixels in the aperture, $f=\sqrt{{{\rm{\Omega }}}_{b}/{{\rm{\Omega }}}_{p}}$ is the noise correlation ratio (e.g., Booth et al. 2017) and N_b is the number of beams in the aperture. The noise ξ is in units of W m⁻² sr⁻¹.

For the lines for which enough flux samples can be collected, we compared the two error estimation methods. We find that the ratio σ_sample/σ_analytical = 0.7 ± 0.1, i.e., the analytical method tends to give larger errors, but the two methods are in reasonable agreement. Finally, we conservatively add a 10% calibration error in quadrature. ¹⁹ Our results are listed in Table 2. Our work results in three new detections of C i toward HD 21997, HD 121191, and HD 121617 (see Figures 1, 2 and 3).

We also measured continuum fluxes with the same procedure as for line emission. For HD 95086 we excluded the background source identified by Zapata et al. (2018) from our aperture, which implies that our measurements will tend to underestimate the total disk flux. Our continuum measurements are presented in Table 11 of Appendix A.

As discussed in Section 2.2.2, we used image cubes without the JvM correction applied. As a test, we computed disk-integrated line fluxes using JvM-corrected image cubes. We find that for all disks and lines, the fluxes agree with the values given in Table 2 within 0.7σ. For the continuum, 84% of the flux measurements are within 1σ, with the maximum difference of 2.2σ occurring for HD 181327. These differences are partially due to the fact that the aperture chosen to measure the flux depends on the noise in the image, and the JvM correction changes the noise. We also find that the derived flux errors agree well. Overall, this indicates that our images have a dirty beam that is well approximated by a Gaussian. For naturally weighted images, there is a larger proportion of cases where the corrected and uncorrected line fluxes differ: 84% are within 1σ, with the maximum difference of 3.0σ being identified for ¹³CO 3–2 toward HD 110058.

We measure fluxes from the C ii spectra extracted from the Herschel/PACS data (see Section 2.3). First, we measure the standard deviation ς of the flux data points, excluding the region λ₀ ± (1.5Δλ) (where again λ₀ and Δλ are the rest wavelength and FWHM of the line, respectively). We define a model for the emission line as

$$\begin{eqnarray}&&F(\lambda ,{F}_{0})=\displaystyle \frac{{F}_{0}}{\sqrt{2\pi }\sigma }\exp \left(-\displaystyle \frac{{\left(\lambda -{\lambda }_{0}\right)}^{2}}{2{\sigma }^{2}}\right),\end{eqnarray} \tag{ 6 }$$

where $\sigma ={\rm{\Delta }}\lambda /(2\sqrt{2\mathrm{ln}2})$ and F₀ is the total line flux. We then consider the likelihood

$$\begin{eqnarray}&&{ \mathcal L }({F}_{0})=\displaystyle \prod _{i}\displaystyle \frac{1}{\sqrt{2\pi }f\varsigma }\exp \left(-\displaystyle \frac{{\left(F({\lambda }_{i},{F}_{0})-{D}_{i}\right)}^{2}}{2{\left(f\varsigma \right)}^{2}}\right),\end{eqnarray} \tag{ 7 }$$

where the product goes over all data points i of the spectrum, D_i are the measured flux values, and $f=\sqrt{{\rm{\Delta }}\lambda /\delta \lambda }$ (with δ λ the channel width) is the noise correlation ratio (e.g., Booth et al. 2017). We verified that the likelihood is very well approximated by a Gaussian. Using ${ \mathcal L }$ we can then determine the maximum likelihood estimate (e.g., Barlow 1989), and error bars corresponding to the fluxes where $\mathrm{ln}{ \mathcal L }$ has decreased by 0.5 with respect to its peak. This yields C ii fluxes of (0.6 ± 2.0) × 10⁻¹⁹ W m⁻² for HD 21997 and (3.3 ± 0.5) × 10⁻¹⁷ W m⁻² for β Pic. For HD 21997, a 3σ upper limit of 6.6 × 10⁻¹⁹ W m⁻² is derived by determining the flux where $\mathrm{ln}{ \mathcal L }$ has decreased by 4.5 with respect to its peak. These maximum likelihood estimates of the C ii fluxes are reported in Table 2 and used in our Markov Chain Monte Carlo (MCMC) modeling. However, we may go one step further and derive a posterior probability distribution by multiplying the likelihood with a prior that is zero for F₀ < 0, and flat for F₀ ≥ 0. For HD 21997, this yields ${1.6}_{-1.1}^{+1.6}\times {10}^{-19}$ W m⁻² (50th percentile with 15.9th and 84.1th percentiles as error bars), or a 99% upper limit of 5.6 × 10⁻¹⁹ W m⁻². For β Pic, the maximum likelihood and the posterior distribution yield identical results because the line is clearly detected.

Brandeker et al. (2016) also analyzed the C ii data of β Pic, but they only report the fluxes from the four spaxels for which they detect emission. Summing their values (with errors summed in quadrature) gives (3.2 ± 0.3) × 10⁻¹⁷ W m⁻², consistent with our result. Note that C ii toward β Pic was also observed by Herschel/HIFI, but the HIFI beam barely covers the disk (Cataldi et al. 2014). Therefore, the data from Herschel/PACS are better suited to measure the total flux.

We first present a model where local thermodynamic equilibrium (LTE) is assumed. The free parameters of the model are the total CO mass M(CO), the total C⁰ mass M(C⁰), the total C⁺ mass M(C⁺) (if C ii was observed), the CO temperature T_CO, and the carbon temperature T_C; see Table 3. Thus, we assume that C⁰ and C⁺ share the same temperature. This assumption is made because half of our disks do not have C ii data. In Appendix B.1, for the disks with C ii data, we run additional fits where the temperatures of C⁰ and C⁺ are separated. No significant changes in the derived gas masses are observed when separating the C⁰ and C⁺ temperatures.

The mass ratios of C⁰/¹³C⁰, ¹²CO/¹³CO, and ¹²CO/C¹⁸O are fixed to interstellar medium (ISM) abundances (¹²C/¹³C = 69, ¹⁶O/¹⁸O = 557, Wilson 1999). With this assumption we explicitly ignore the possibility of CO isotopologue-selective photodissociation which would tend to increase the ¹²CO/¹³CO and ¹²CO/C¹⁸O ratios (e.g., Visser et al. 2009). However, the assumption is necessary to constrain the CO mass in the cases where the ¹²CO emission is optically thick, although it might lead to an underestimation of the CO mass if isotopologue-selective photodissociation is active.

We assume a disk with a uniform density for r_in ≤ r ≤ r_out (see Table 1 for the adopted values), and within ±1.5H of the disk midplane, where H is the scale height given by Equation (3). For simplicity, the scale height is fixed, that is, it does not vary with radius nor with the temperature: r = r_mean and T = 50 K are fixed to calculate H. Choosing T = 50 K is arbitrary, but given the weak dependence of H on temperature, we do not expect this choice to affect our results. Indeed, compared to H(T = 50 K), H would vary by a factor ∼2 at most for the temperature range (10–200 K) we consider. We also note that in reality, the vertical structure is unlikely to be the same for each species as it depends on a number of factors and the interplay between photochemistry and vertical diffusion (e.g., Marino et al. 2022).

For a given set of parameter values, we can then compute the total line emission by ray tracing the emission along the line of sight of the 3D disk model (e.g., Cataldi et al. 2014), assuming a Keplerian velocity field that is considered independent of z (the direction perpendicular to the midplane). However, for low gas densities, the gas might not be in Keplerian rotation, but rather be blown out as a wind (Kral et al. 2023). As discussed in Section 6.4, this caveat in unlikely to affect our results. Since we assume LTE, the level populations are computed from the Boltzmann equation. The atomic data is taken from the LAMDA database. ²⁰ The local intrinsic emission was assumed as a square line profile with a fixed width of 0.5 km s⁻¹. In Appendix B.2, we discuss how our fits change when the line width is changed.

We initially explored an even simpler model that ignored the velocity field of the disk, but found that such a model can easily overestimate the optical depth of the lines, resulting in highly uncertain column density estimates, particularly for highly inclined disks.

The parameter space is explored using the MCMC method implemented by the emcee package (Foreman-Mackey et al. 2013). The log-likelihood is given by

$$\begin{eqnarray}&&\mathrm{ln}{ \mathcal L }({ \mathcal P })=-\displaystyle \sum _{i}\displaystyle \frac{{\left({F}_{i}({ \mathcal P })-{D}_{i}\right)}^{2}}{2{\sigma }_{i}^{2}},\end{eqnarray} \tag{ 8 }$$

where the sum runs over all observed emission lines, ${ \mathcal P }$ represents the parameter values, F_i is the model flux, and D_i and σ_i are the observed flux and error as given in Table 2. The posterior probability is found by multiplying the likelihood ${ \mathcal L }$ with a flat prior for each parameter, with limits detailed in Table 3. The influence of the chosen temperature priors on our results is discussed in Appendix B.3.

Table 3. Parameters Used in Our Model Fitting

Parameter	Units	min	max
(1)	(2)	(3)	(4)
log M(CO)	$\mathrm{log}{M}_{\oplus }$	−8	1
log M(C0)	$\mathrm{log}{M}_{\oplus }$	−8	1
log M(C+) a	$\mathrm{log}{M}_{\oplus }$	−8	1
TCO	K	10	200
TC	K	10	200
log n(H2) b	log cm−3	1	7
log n(e−) b	log cm−3	−1	5

Notes. (3) Value below which the flat prior is zero. (4) Value above which the flat prior is zero.

^a Used only if C ii data were available. ^b Used for the non-LTE models only.

Download table as:  ASCIITypeset image

An additional prior is applied that ensures T_CO ≤ T_C. This represents our expectation that C is either well mixed with CO (T_CO = T_C) or preferentially occupies regions where CO is efficiently photodissociated such as the disk atmosphere (Marino et al. 2022) or the inner disk regions where the temperature is expected to be higher (e.g., Kral et al. 2019). Note though that our disk model implicitly assumes that the gas is well mixed by considering an identical spatial distribution for each species.

For computational efficiency, we precompute the model fluxes on a grid of parameter values and interpolate on that grid during the MCMC run. The grid has a resolution of 5 K (for T < 30 K), respectively, 25 K (for T > 30 K) for temperature, at least 0.5 dex for gas masses, and 0.5 dex for collider number densities (used for the non-LTE models, see below). We tested the grid by comparing interpolated and directly calculated fluxes and found excellent agreement. For the MCMC, we employ 1000 walkers, each taking 20,000 steps. We discard the first 15,000 steps as burn-in.

We explore additional models where the LTE assumption is dropped. We introduce an additional free parameter: either the H₂ number density n(H₂), or the electron number density n(e⁻). Note that in the latter case, we do not assume n(C⁺) = n(e⁻), that is, n(e⁻) is an independent parameter. The level populations are determined by using a statistical equilibrium calculation in the radiative transfer code pythonradex, ²¹ where the input for the pythonradex code is the column density along the line of sight N_los corresponding to a given gas mass. The radiative transfer is then calculated in the same way as for the LTE models.

In Table 4 we present summary statistics of the posterior probability distributions we derived for each parameter in the LTE case. We also transformed the derived masses into column densities. We derive the averaged column density in the z direction N_⊥ by dividing the mass by the face-on projection of the disk $\pi ({r}_{\mathrm{out}}^{2}-{r}_{\mathrm{in}}^{2})$ and the molecular mass. We also derive the averaged column density along the line of sight N_los by considering the area of the disk projected onto the sky plane. To determine the latter, we project our 3D disk model along the line of sight and then sum the area of all nonzero pixels.

Table 4. Median as Well as the 16th and 84th Percentiles of the Posterior Probability Distribution Derived with our MCMC Runs Assuming LTE

Star	$\mathrm{log}M(\mathrm{CO})$	$\mathrm{log}M({{\rm{C}}}^{0})$	$\mathrm{log}M({{\rm{C}}}^{+})$	TCO	TC	$\mathrm{log}{N}_{\perp }(\mathrm{CO})$	$\mathrm{log}{N}_{\perp }({{\rm{C}}}^{0})$	$\mathrm{log}{N}_{\perp }({{\rm{C}}}^{+})$	$\mathrm{log}{N}_{\mathrm{los}}(\mathrm{CO})$	$\mathrm{log}{N}_{\mathrm{los}}({{\rm{C}}}^{0})$	$\mathrm{log}{N}_{\mathrm{los}}({{\rm{C}}}^{+})$
	($\mathrm{log}{M}_{\oplus }$)	($\mathrm{log}{M}_{\oplus }$)	($\mathrm{log}{M}_{\oplus }$)	(K)	(K)	($\mathrm{log}\,{\mathrm{cm}}^{-2}$)	($\mathrm{log}\,{\mathrm{cm}}^{-2}$)	($\mathrm{log}\,{\mathrm{cm}}^{-2}$)	($\mathrm{log}\,{\mathrm{cm}}^{-2}$)	($\mathrm{log}\,{\mathrm{cm}}^{-2}$)	($\mathrm{log}\,{\mathrm{cm}}^{-2}$)
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)	(10)	(11)	(12)
49 Ceti	$-{2.15}_{-0.07}^{+0.07}$	$-{2.05}_{-0.17}^{+0.2}$	$-{1.1}_{-1.0}^{+0.8}$	${15.2}_{-0.6}^{+0.6}$	${24}_{-2}^{+3}$	${16.74}_{-0.07}^{+0.07}$	${17.21}_{-0.17}^{+0.2}$	${18.1}_{-1.0}^{+0.8}$	${17.41}_{-0.07}^{+0.07}$	${17.88}_{-0.17}^{+0.2}$	${18.8}_{-1.0}^{+0.8}$
			>−3.6					>15.7			>16.3
β Pictoris	$-{4.32}_{-0.15}^{+0.10}$	$-{3.75}_{-0.09}^{+0.08}$	$-{3.45}_{-0.16}^{+1.5}$	${11.8}_{-1.1}^{+2}$	${92}_{-60}^{+70}$	${14.58}_{-0.15}^{+0.10}$	${15.52}_{-0.09}^{+0.08}$	${15.82}_{-0.16}^{+1.5}$	${15.59}_{-0.15}^{+0.10}$	${16.52}_{-0.09}^{+0.08}$	${16.82}_{-0.16}^{+1.5}$
			>−3.7					>15.5			>16.5
HD 21997	$-{1.32}_{-0.09}^{+0.07}$	$-{2.90}_{-0.09}^{+1.0}$	<−0.33	${10.25}_{-0.18}^{+0.3}$	${67}_{-50}^{+90}$	${17.71}_{-0.09}^{+0.07}$	${16.50}_{-0.09}^{+1.0}$	<19.1	${17.77}_{-0.09}^{+0.07}$	${16.55}_{-0.09}^{+1.0}$	<19.1
		>−3.1					>16.3			>16.3
HD 32297	$-{1.28}_{-0.15}^{+0.12}$	$-{2.53}_{-0.05}^{+0.08}$	$-{2.4}_{-0.5}^{+1.6}$	${25.3}_{-1.8}^{+1.8}$	${56}_{-16}^{+90}$	${17.72}_{-0.15}^{+0.12}$	${16.84}_{-0.05}^{+0.08}$	${16.9}_{-0.5}^{+1.6}$	${18.26}_{-0.15}^{+0.12}$	${17.37}_{-0.05}^{+0.08}$	${17.5}_{-0.5}^{+1.6}$
			>−3.8					>15.6			>16.1
HD 48370	<−4.9	<−4.4	<−3.9 a	${69}_{-40}^{+60}$	${150}_{-60}^{+40}$	<14.2	<15.1	<15.5 a	<14.5	<15.4	<15.9 a
HD 61005	<−5.4	<−4.9	<−4.6 a	${69}_{-40}^{+60}$	${150}_{-60}^{+40}$	<14.2	<15.1	<15.4 a	<15.1	<15.9	<16.3 a
HD 95086	<−5.4	<−3.9	<−3.0 a	${70}_{-40}^{+60}$	${150}_{-60}^{+40}$	<12.9	<14.8	<15.6 a	<12.9	<14.8	<15.6 a
HD 110058	$-{2.03}_{-0.12}^{+0.14}$	<−3.5	<−4.1 a	${102}_{-8}^{+8}$	${150}_{-30}^{+30}$	${18.63}_{-0.12}^{+0.14}$	<17.6	<16.9 a	${19.62}_{-0.12}^{+0.14}$	<18.5	<17.9 a
HD 121191	$-{2.25}_{-0.3}^{+0.20}$	$-{3.68}_{-0.13}^{+0.11}$	$-{5.01}_{-0.04}^{+0.03}$ a	${91}_{-8}^{+9}$	${140}_{-40}^{+40}$	${18.16}_{-0.3}^{+0.20}$	${17.10}_{-0.13}^{+0.11}$	${15.78}_{-0.04}^{+0.03}$ a	${18.21}_{-0.3}^{+0.20}$	${17.15}_{-0.13}^{+0.11}$	${15.82}_{-0.04}^{+0.03}$ a
HD 121617	$-{1.59}_{-0.07}^{+0.07}$	$-{2.42}_{-0.06}^{+0.12}$	$-{3.49}_{-0.09}^{+0.07}$ a	${42}_{-4}^{+4}$	${100}_{-50}^{+70}$	${17.78}_{-0.07}^{+0.07}$	${17.32}_{-0.06}^{+0.12}$	${16.24}_{-0.09}^{+0.07}$ a	${17.85}_{-0.07}^{+0.07}$	${17.39}_{-0.06}^{+0.12}$	${16.31}_{-0.09}^{+0.07}$ a
HD 131835	$-{1.76}_{-0.08}^{+0.06}$	$-{2.48}_{-0.11}^{+0.8}$	<−0.45	${12.0}_{-0.3}^{+0.3}$	${74}_{-60}^{+90}$	${17.02}_{-0.08}^{+0.06}$	${16.67}_{-0.11}^{+0.8}$	<18.7	${17.41}_{-0.08}^{+0.06}$	${17.05}_{-0.11}^{+0.8}$	<19.1
HD 146897	$-{4.4}_{-0.2}^{+1.2}$	<−3.9	<−4.2 a	${43}_{-30}^{+70}$	${140}_{-60}^{+50}$	${15.3}_{-0.2}^{+1.2}$	<16.2	<15.8 a	${15.8}_{-0.2}^{+1.2}$	<16.7	<16.4 a
HD 181327	$-{5.72}_{-0.12}^{+0.3}$	<−4.5	<−1.4	${39}_{-20}^{+60}$	${130}_{-70}^{+50}$	${14.00}_{-0.12}^{+0.3}$	<15.6	<18.7	${13.99}_{-0.12}^{+0.3}$	<15.6	<18.7
HR 4796	<−6.0	<−4.0	<−0.75	${65}_{-40}^{+60}$	${140}_{-70}^{+40}$	<13.9	<16.3	<19.5	<14.2	<16.6	<19.8

Notes. Upper limits corresponding to the 99th percentile are given if the disk-integrated flux has an S/N < 3. Lower limits corresponding to the 1st percentile are given if Equation (9) is satisfied. (7), (8), (9): column density in the z direction (perpendicular to the midplane). (10), (11), (12): column density along the line of sight.

^a Estimated from an ionization calculation, because no C ii data were available. See the text for details.

Download table as:  ASCIITypeset image

In Table 4, we list upper limits (99th percentile) for the mass and column density whenever the corresponding disk-integrated fluxes have an S/N < 3. The table also lists lower limits whenever a parameter has a substantial probability to be close to the upper bound of the prior. This is formalized as

$$\begin{eqnarray}&&{\int }_{{p}_{\max }-0.2W}^{{p}_{\max }}P(p){dp}\gt 5 \% ,\end{eqnarray} \tag{ 9 }$$

where p_max is the maximum parameter value allowed by the prior, W is a proxy for the width of the posterior distribution calculated as the difference between the 99th and the 1st percentiles, and P is the posterior probability density.

If the C ii line was not observed, we estimate the amount of ionized carbon with a simple ionization equilibrium calculation. For a given C⁰ mass, we calculate the ionization balance by assuming uniform density and that the electron density equals the C⁺ density. We use ionization cross sections (Nahar & Pradhan 1991) and recombination coefficients (Nahar 1995, 1996) from the NORAD-Atomic-Data database. ²² For the stellar flux, we use ATLAS9 models (Castelli & Kurucz 2003) with the parameters listed in Table 1 and scaled by $\exp (-{\alpha }_{{\rm{C}}}n({{\rm{C}}}^{0})({r}_{\mathrm{mean}}-{r}_{\mathrm{in}}))$ to take into account extinction by C⁰ ionization. Here, α_C is the ionization cross section and n(C⁰) the uniform C⁰ number density. For the interstellar radiation, we use the Draine field (Draine 1978; Lee 1984) scaled by $\exp (-{\alpha }_{{\rm{C}}}{N}_{\perp }({{\rm{C}}}^{0})/2)$. This calculation implicitly assumes that C⁰ and C⁺ are well mixed. Marino et al. (2022) found that C⁺ tends to form an upper layer in the disk independently of the strength of vertical mixing. Therefore, our simple model may underestimate the total C⁺ mass because it only accounts for the C⁺ colocated with C⁰.

In Table 5, we present the optical depth we derived for each emission line.

Table 5. Median as Well as the 16th and 84th Percentiles of the Optical Depth Along the Line of Sight, Derived from MCMC Runs Assuming LTE

Star	Emission Line	τ(ν0)
49 Ceti	C i 3P1–3P0	${11}_{-5}^{+8}$
	13C i 3P1–3P0	${0.15}_{-0.06}^{+0.09}$
	C ii 2P3/2–2P1/2	${110}_{-100}^{+600}$ (>0.084)
	12CO 2–1	${121}_{-14}^{+16}$
	12CO 3–2	${120}_{-20}^{+20}$
	13CO 2–1	${1.53}_{-0.18}^{+0.20}$
	C18O 2–1	${0.20}_{-0.02}^{+0.03}$
β Pictoris	C i 3P1–3P0	${0.080}_{-0.04}^{+0.3}$
	C ii 2P3/2–2P1/2	${0.45}_{-0.3}^{+30}$ (>0.12)
	12CO 2–1	${2.5}_{-0.9}^{+1.3}$
	12CO 3–2	${1.8}_{-0.5}^{+0.3}$
HD 21997	C i 3P1–3P0	${0.10}_{-0.07}^{+10}$ (>0.025)
	C ii 2P3/2–2P1/2	<260
	12CO 2–1	${450}_{-90}^{+60}$
	12CO 3–2	${280}_{-70}^{+60}$
	13CO 2–1	${6.0}_{-1.4}^{+1.0}$
	13CO 3–2	${3.7}_{-1.0}^{+0.8}$
	C18O 2–1	${0.64}_{-0.09}^{+0.18}$
	C18O 3–2	${0.46}_{-0.07}^{+0.05}$
HD 32297	C i 3P1–3P0	${1.2}_{-0.9}^{+1.6}$
	C ii 2P3/2–2P1/2	${3.9}_{-3}^{+200}$ (>0.087)
	12CO 2–1	${580}_{-200}^{+200}$
	13CO 2–1	${7.1}_{-2}^{+2}$
	C18O 2–1	${0.86}_{-0.3}^{+0.5}$
HD 48370	C i 3P1–3P0	<6.1 × 10−3
	12CO 2–1	<0.042
	13CO 2–1	<7.0 × 10−4
	C18O 2–1	<3.3 × 10−4
HD 61005	C i 3P1–3P0	<0.022
	12CO 2–1	<0.18
HD 95086	C i 3P1–3P0	<1.2 × 10−3
	12CO 1–0	<4.0 × 10−4
	12CO 2–1	<9.2 × 10−4
	12CO 3–2	<9.6 × 10−4
HD 110058	C i 3P1–3P0	<5.9
	12CO 2–1	${940}_{-200}^{+500}$
	12CO 3–2	(${1.8}_{-0.4}^{+1.2}$) × 10−3
	13CO 2–1	${12}_{-4}^{+7}$
	13CO 3–2	${24}_{-8}^{+12}$
HD 121191	C i 3P1–3P0	${0.16}_{-0.06}^{+0.09}$
	12CO 2–1	${39}_{-18}^{+20}$
	13CO 2–1	${0.48}_{-0.2}^{+0.3}$
	C18O 2–1	${0.060}_{-0.03}^{+0.04}$
HD 121617	C i 3P1–3P0	${0.42}_{-0.2}^{+1.0}$
	12CO 2–1	${83}_{-15}^{+16}$
	13CO 2–1	${1.1}_{-0.2}^{+0.2}$
	C18O 2–1	${0.13}_{-0.03}^{+0.03}$
HD 131835	C i 3P1–3P0	${0.30}_{-0.18}^{+19}$
	C ii 2P3/2–2P1/2	<240
	12CO 2–1	${170}_{-40}^{+40}$
	12CO 3–2	${114}_{-13}^{+11}$
	13CO 3–2	${1.52}_{-0.17}^{+0.2}$
	C18O 3–2	${0.21}_{-0.05}^{+0.04}$
HD 146897	C i 3P1–3P0	<0.37
	12CO 2–1	${0.98}_{-0.7}^{+130}$
HD 181327	C i 3P1–3P0	<0.016
	C ii 2P3/2–2P1/2	<78
	12CO 2–1	${0.012}_{-0.008}^{+0.02}$
	12CO 3–2	${0.017}_{-0.010}^{+0.015}$
	13CO 2–1	$({2.7}_{-2}^{+6})\times {10}^{-4}$
HR 4796	C i 3P1–3P0	<0.20
	C ii 2P3/2–2P1/2	<1.8 × 103
	12CO 2–1	<0.11
	12CO 3–2	<0.094
	12CO 6–5	<0.025

Note. Upper limits correspond to the 99th percentile and are given if the S/N of the disk-integrated flux is less than 3. Lower limits corresponding to the 1st percentile are given in parenthesis if Equation (9) is satisfied.

Download table as:  ASCIITypeset images: 1 2

Figures 4 and 5 show corner plots (1D and 2D histograms of the posterior probability distribution) for two example disks: the CO-rich disk around HD 32297 and the disk around HD 48370 where no gas emission was detected. For HD 32297, we see that the CO and C⁰ masses are well constrained, while the C⁺ mass could be very large if the temperature T_C is low. Indeed, for low temperatures, the C ii emission becomes optically thick and the mass cannot be meaningfully constrained anymore. High optical depth at low temperature and consequently large uncertainties on the underlying CO, C⁰, or C⁺ masses are an issue for several other disks. The CO temperature in the HD 32997 disk is well constrained. The posterior of the carbon temperature has a distinct peak, but also broad wings.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

View figure set (14 panels)

Tarball of figure set images

Zoom In Zoom Out Reset image size

For the HD 48370 disk, the posterior distributions of the masses clearly allow us to derive upper limits, while the temperatures are essentially unconstrained. Note that the shape of the posterior distributions of T_CO and T_C reflects our prior condition that T_CO < T_C.

The complete set of corner plots is available in the online journal. We emphasize that the values presented in Table 4 are merely summary statistics of the posterior distributions. The reader interested in a particular parameter or disk is invited to look at the corner plots to fully understand the posterior distributions and their mutual correlations.

Figure 6 compares the observed line fluxes F_obs to the median of the model line fluxes F_model, normalized by the observational error. We see that our model, in general, successfully reproduces the observed line emission. There are some exceptions, such as HD 21997 where the model underpredicts the ¹²CO 2–1 flux by 3.5σ and overpredicts the 2–1, and 3–2 fluxes of ¹³CO by 2σ.

Zoom In Zoom Out Reset image size

The derived masses, temperatures, and column densities for the non-LTE cases are presented in Tables 6 and 7 for collisions with H₂ and e⁻, respectively. The two cases deliver consistent gas masses, although we note that the case with e⁻ collisions generally gives less constraining results for the C⁰ masses (e.g., higher upper limits).

Table 6. Same as Table 4, but for the Non-LTE Case Assuming Collisions with H₂

Star	$\mathrm{log}M(\mathrm{CO})$	$\mathrm{log}M({{\rm{C}}}^{0})$	$\mathrm{log}M({{\rm{C}}}^{+})$	TCO	TC	$\mathrm{log}{N}_{\perp }(\mathrm{CO})$	$\mathrm{log}{N}_{\perp }({{\rm{C}}}^{0})$	$\mathrm{log}{N}_{\perp }({{\rm{C}}}^{+})$	$\mathrm{log}{N}_{\mathrm{los}}(\mathrm{CO})$	$\mathrm{log}{N}_{\mathrm{los}}({{\rm{C}}}^{0})$	$\mathrm{log}{N}_{\mathrm{los}}({{\rm{C}}}^{+})$
	($\mathrm{log}{M}_{\oplus }$)	($\mathrm{log}{M}_{\oplus }$)	($\mathrm{log}{M}_{\oplus }$)	(K)	(K)	($\mathrm{log}\,{\mathrm{cm}}^{-2}$)	($\mathrm{log}\,{\mathrm{cm}}^{-2}$)	($\mathrm{log}\,{\mathrm{cm}}^{-2}$)	($\mathrm{log}\,{\mathrm{cm}}^{-2}$)	($\mathrm{log}\,{\mathrm{cm}}^{-2}$)	($\mathrm{log}\,{\mathrm{cm}}^{-2}$)
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)	(10)	(11)	(12)
49 Ceti	$-{2.04}_{-0.15}^{+0.5}$	$-{2.0}_{-0.2}^{+0.3}$	$-{1.1}_{-1.1}^{+0.7}$	${15.7}_{-0.8}^{+1.1}$	${24}_{-2}^{+6}$	${16.86}_{-0.15}^{+0.5}$	${17.2}_{-0.2}^{+0.3}$	${18.2}_{-1.1}^{+0.7}$	${17.52}_{-0.15}^{+0.5}$	${17.9}_{-0.2}^{+0.3}$	${18.9}_{-1.1}^{+0.7}$
			>−3.5					>15.7			>16.4
β Pictoris	$-{4.36}_{-0.14}^{+0.14}$	$-{3.77}_{-0.09}^{+0.09}$	$-{3.2}_{-0.3}^{+0.5}$	${38}_{-30}^{+70}$	${130}_{-80}^{+50}$	${14.53}_{-0.14}^{+0.14}$	${15.50}_{-0.09}^{+0.09}$	${16.1}_{-0.3}^{+0.5}$	${15.54}_{-0.14}^{+0.14}$	${16.50}_{-0.09}^{+0.09}$	${17.1}_{-0.3}^{+0.5}$
HD 21997	$-{1.17}_{-0.14}^{+0.3}$	$-{2.91}_{-0.09}^{+0.9}$	<−0.33	${10.7}_{-0.5}^{+0.6}$	${72}_{-60}^{+90}$	${17.85}_{-0.14}^{+0.3}$	${16.48}_{-0.09}^{+0.9}$	<19.1	${17.91}_{-0.14}^{+0.3}$	${16.54}_{-0.09}^{+0.9}$	<19.1
HD 32297	$-{0.5}_{-0.6}^{+0.4}$	$-{2.46}_{-0.11}^{+0.19}$	$-{1.8}_{-0.7}^{+1.0}$	${26.3}_{-1.8}^{+1.8}$	${79}_{-30}^{+80}$	${18.5}_{-0.6}^{+0.4}$	${16.91}_{-0.11}^{+0.19}$	${17.6}_{-0.7}^{+1.0}$	${19.0}_{-0.6}^{+0.4}$	${17.44}_{-0.11}^{+0.19}$	${18.1}_{-0.7}^{+1.0}$
	>−1.5		>−3.4			>17.5		>15.9	>18.0		>16.5
HD 48370	<−4.1	<−4.1	<−3.7 a	${66}_{-40}^{+60}$	${140}_{-60}^{+40}$	<14.9	<15.3	<15.7 a	<15.2	<15.7	<16.0 a
HD 61005	<−4.7	<−4.7	<−4.5 a	${66}_{-40}^{+60}$	${140}_{-60}^{+40}$	<15.0	<15.3	<15.5 a	<15.8	<16.1	<16.4 a
HD 95086	<−4.5	<−3.5	<−2.9 a	${65}_{-40}^{+60}$	${140}_{-60}^{+40}$	<13.7	<15.1	<15.7 a	<13.8	<15.1	<15.7 a
HD 110058	$-{1.8}_{-0.3}^{+1.2}$	<−3.0	<−4.0 a	${104}_{-9}^{+9}$	${150}_{-30}^{+30}$	${18.8}_{-0.3}^{+1.2}$	<18.0	<17.0 a	${19.8}_{-0.3}^{+1.2}$	<19.0	<18.0 a
	>−2.3					>18.4			>19.4
HD 121191	$-{2.1}_{-0.4}^{+1.3}$	$-{3.61}_{-0.17}^{+0.7}$	$-{5.02}_{-0.2}^{+0.05}$ a	${102}_{-14}^{+30}$	${160}_{-40}^{+30}$	${18.3}_{-0.4}^{+1.3}$	${17.17}_{-0.17}^{+0.7}$	${15.76}_{-0.2}^{+0.05}$ a	${18.4}_{-0.4}^{+1.3}$	${17.21}_{-0.17}^{+0.7}$	${15.80}_{-0.2}^{+0.05}$ a
	>−3.2					>17.2			>17.3
HD 121617	$-{1.62}_{-0.09}^{+0.09}$	$-{2.42}_{-0.06}^{+0.14}$	$-{3.49}_{-0.09}^{+0.07}$ a	${43}_{-4}^{+4}$	${100}_{-50}^{+60}$	${17.75}_{-0.09}^{+0.09}$	${17.32}_{-0.06}^{+0.14}$	${16.25}_{-0.09}^{+0.07}$ a	${17.82}_{-0.09}^{+0.09}$	${17.39}_{-0.06}^{+0.14}$	${16.32}_{-0.09}^{+0.07}$ a
HD 131835	$-{1.70}_{-0.13}^{+0.7}$	$-{2.48}_{-0.12}^{+0.8}$	<−0.53	${12.0}_{-0.3}^{+0.5}$	${82}_{-70}^{+80}$	${17.08}_{-0.13}^{+0.7}$	${16.68}_{-0.12}^{+0.8}$	<18.6	${17.46}_{-0.13}^{+0.7}$	${17.06}_{-0.12}^{+0.8}$	<19.0
	>−3.3					>15.5			>15.9
HD 146897	$-{4.2}_{-0.4}^{+1.5}$	<−3.6	<−4.1 a	${45}_{-30}^{+70}$	${140}_{-60}^{+40}$	${15.5}_{-0.4}^{+1.5}$	<16.5	<15.9 a	${16.0}_{-0.4}^{+1.5}$	<17.0	<16.5 a
HD 181327	$-{5.3}_{-0.5}^{+0.9}$	<−3.9	<−1.2	${54}_{-30}^{+60}$	${140}_{-60}^{+40}$	${14.4}_{-0.5}^{+0.9}$	<16.1	<18.9	${14.4}_{-0.5}^{+0.9}$	<16.1	<18.9
HR 4796	<−4.1	<−3.5	<−0.66	${59}_{-40}^{+60}$	${140}_{-60}^{+40}$	<15.8	<16.8	<19.6	<16.1	<17.1	<19.9

Note.

^a Estimated from an ionization calculation because no C ii data was available. See the text for details.

Download table as:  ASCIITypeset image

Table 7. Same as Table 4, but for the Non-LTE Case Assuming Collisions with e⁻

Star	$\mathrm{log}M(\mathrm{CO})$	$\mathrm{log}M({{\rm{C}}}^{0})$	$\mathrm{log}M({{\rm{C}}}^{+})$	TCO	TC	$\mathrm{log}{N}_{\perp }(\mathrm{CO})$	$\mathrm{log}{N}_{\perp }({{\rm{C}}}^{0})$	$\mathrm{log}{N}_{\perp }({{\rm{C}}}^{+})$	$\mathrm{log}{N}_{\mathrm{los}}(\mathrm{CO})$	$\mathrm{log}{N}_{\mathrm{los}}({{\rm{C}}}^{0})$	$\mathrm{log}{N}_{\mathrm{los}}({{\rm{C}}}^{+})$
	($\mathrm{log}{M}_{\oplus }$)	($\mathrm{log}{M}_{\oplus }$)	($\mathrm{log}{M}_{\oplus }$)	(K)	(K)	($\mathrm{log}\,{\mathrm{cm}}^{-2}$)	($\mathrm{log}\,{\mathrm{cm}}^{-2}$)	($\mathrm{log}\,{\mathrm{cm}}^{-2}$)	($\mathrm{log}\,{\mathrm{cm}}^{-2}$)	($\mathrm{log}\,{\mathrm{cm}}^{-2}$)	($\mathrm{log}\,{\mathrm{cm}}^{-2}$)
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)	(10)	(11)	(12)
49 Ceti	$-{1.7}_{-0.4}^{+0.2}$	$-{1.0}_{-0.9}^{+0.4}$	$-{1.2}_{-0.9}^{+0.8}$	${17.0}_{-1.6}^{+1.6}$	${24.4}_{-1.9}^{+2}$	${17.2}_{-0.4}^{+0.2}$	${18.2}_{-0.9}^{+0.4}$	${18.1}_{-0.9}^{+0.8}$	${17.8}_{-0.4}^{+0.2}$	${18.9}_{-0.9}^{+0.4}$	${18.8}_{-0.9}^{+0.8}$
		>−2.5	>−3.6				>16.8	>15.7		>17.4	>16.4
β Pictoris	$-{4.43}_{-0.12}^{+0.12}$	$-{3.77}_{-0.09}^{+0.09}$	$-{3.52}_{-0.11}^{+0.4}$	${44}_{-30}^{+70}$	${130}_{-70}^{+50}$	${14.47}_{-0.12}^{+0.12}$	${15.50}_{-0.09}^{+0.09}$	${15.75}_{-0.11}^{+0.4}$	${15.48}_{-0.12}^{+0.12}$	${16.50}_{-0.09}^{+0.09}$	${16.75}_{-0.11}^{+0.4}$
HD 21997	$-{1.22}_{-0.10}^{+0.14}$	$-{2.91}_{-0.10}^{+1.3}$	<−0.35	${10.7}_{-0.4}^{+0.6}$	${71}_{-60}^{+90}$	${17.81}_{-0.10}^{+0.14}$	${16.48}_{-0.10}^{+1.3}$	<19.0	${17.87}_{-0.10}^{+0.14}$	${16.54}_{-0.10}^{+1.3}$	<19.1
		>−3.1					>16.3			>16.3
HD 32297	$-{0.5}_{-0.5}^{+0.4}$	$-{1.9}_{-0.5}^{+0.5}$	$-{2.5}_{-0.4}^{+1.5}$	${26.7}_{-1.9}^{+1.9}$	${66}_{-20}^{+90}$	${18.5}_{-0.5}^{+0.4}$	${17.5}_{-0.5}^{+0.5}$	${16.9}_{-0.4}^{+1.5}$	${19.0}_{-0.5}^{+0.4}$	${18.0}_{-0.5}^{+0.5}$	${17.4}_{-0.4}^{+1.5}$
	>−1.4		>−3.6			>17.6		>15.8	>18.1		>16.3
HD 48370	<−3.6	<−2.9	<−3.2 a	${65}_{-40}^{+60}$	${140}_{-60}^{+40}$	<15.5	<16.5	<16.2 a	<15.8	<16.9	<16.5 a
HD 61005	<−4.1	<−3.7	<−4.0 a	${65}_{-40}^{+60}$	${140}_{-60}^{+40}$	<15.5	<16.3	<16.0 a	<16.4	<17.2	<16.8 a
HD 95086	<−4.4	<−2.2	<−2.3 a	${66}_{-40}^{+60}$	${140}_{-60}^{+40}$	<13.9	<16.4	<16.3 a	<13.9	<16.4	<16.3 a
HD 110058	$-{1.7}_{-0.4}^{+1.2}$	<−2.8	<−4.0 a	${105}_{-9}^{+9}$	${150}_{-30}^{+30}$	${18.9}_{-0.4}^{+1.2}$	<18.2	<17.0 a	${19.9}_{-0.4}^{+1.2}$	<19.2	<18.0 a
	>−2.3					>18.4			>19.4
HD 121191	$-{2.1}_{-0.5}^{+1.4}$	$-{3.56}_{-0.20}^{+1.1}$	$-{5.02}_{-0.2}^{+0.05}$ a	${111}_{-20}^{+30}$	${160}_{-40}^{+30}$	${18.4}_{-0.5}^{+1.4}$	${17.22}_{-0.20}^{+1.1}$	${15.76}_{-0.2}^{+0.05}$ a	${18.4}_{-0.5}^{+1.4}$	${17.26}_{-0.20}^{+1.1}$	${15.80}_{-0.2}^{+0.05}$ a
	>−3.1					>17.3			>17.3
HD 121617	$-{1.64}_{-0.09}^{+0.09}$	$-{2.41}_{-0.07}^{+0.19}$	$-{3.49}_{-0.09}^{+0.07}$ a	${43}_{-4}^{+5}$	${100}_{-50}^{+60}$	${17.72}_{-0.09}^{+0.09}$	${17.32}_{-0.07}^{+0.19}$	${16.25}_{-0.09}^{+0.07}$ a	${17.79}_{-0.09}^{+0.09}$	${17.39}_{-0.07}^{+0.19}$	${16.32}_{-0.09}^{+0.07}$ a
HD 131835	$-{1.69}_{-0.11}^{+0.5}$	$-{2.45}_{-0.14}^{+1.1}$	<−0.46	${12.0}_{-0.3}^{+0.4}$	${75}_{-60}^{+80}$	${17.10}_{-0.11}^{+0.5}$	${16.70}_{-0.14}^{+1.1}$	<18.7	${17.48}_{-0.11}^{+0.5}$	${17.09}_{-0.14}^{+1.1}$	<19.1
	>−3.4	>−2.7				>15.4	>16.4		>15.8	>16.8
HD 146897	$-{3.6}_{-1.0}^{+1.5}$	<−2.4	<−4.0 a	${47}_{-30}^{+70}$	${140}_{-60}^{+40}$	${16.1}_{-1.0}^{+1.5}$	<17.7	<16.0 a	${16.6}_{-1.0}^{+1.5}$	<18.2	<16.6 a
HD 181327	$-{4.4}_{-1.3}^{+0.7}$	<−2.5	<−1.8	${59}_{-40}^{+60}$	${140}_{-60}^{+40}$	${15.3}_{-1.3}^{+0.7}$	<17.6	<18.3	${15.3}_{-1.3}^{+0.7}$	<17.6	<18.3
HR 4796	<−3.4	<−2.2	<−0.66	${60}_{-40}^{+60}$	${140}_{-70}^{+40}$	<16.6	<18.1	<19.6	<16.9	<18.4	<19.9

Note.

^a Estimated from an ionization calculation, because no C ii data was available. See the text for details.

Download table as:  ASCIITypeset image

We find that in general, n(H₂) and n(e⁻) are not well constrained. By visual inspection of the posterior distributions, we identify the disks where lower limits can be inferred. Although not particularly constraining, they are presented in Table 8 for completeness. The strongest constraints can be derived for HD 21997, for which we show the corner plot in Figure 7 for the case of collisions with H₂ (the full set of corner plots for both H₂ and e⁻ collisions is available in the online journal). These collider densities are derived by assuming a single collider species, either H₂ or e⁻. In reality, several colliders might contribute to the excitation. Thus, our approach might overestimate the collider densities of H₂ and e⁻. In Table 8, we also list n(C⁺) as estimated from our LTE fits, for comparison with n(e⁻). The values are formally consistent, that is, it could indeed be the case that n(e⁻) ≈ n(C⁺), as is often assumed in simple ionization calculations.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

View figure set (28 panels)

Tarball of figure set images

Table 8. Lower Limits on the H₂ and e⁻ Density for Disks Where Meaningful Constraints Can Be Derived

Star	$\mathrm{log}n({{\rm{H}}}_{2})$	$\mathrm{log}n({{\rm{e}}}^{-})$	$\mathrm{log}n({{\rm{C}}}^{+})$
	($\mathrm{log}\,{\mathrm{cm}}^{-3}$)	($\mathrm{log}\,{\mathrm{cm}}^{-3}$)	($\mathrm{log}\,{\mathrm{cm}}^{-3}$)
49 Ceti	>2.1	>0.3	>1.4
β Pictoris		>−0.7	>1.2
HD 21997	>2.8	>1.7	<4.8
HD 32297	>1.8	>0.4	>1.2
HD 110058		>0.4	<3.8 a
HD 121191		>−0.1	${2.43}_{-0.04}^{+0.03}$ a
HD 121617	>1.9	>0.8	${2.14}_{-0.09}^{+0.07}$ a
HD 131835		>0.6	<4.3

Notes. We note that these values are derived from models that assume a single collider, either H₂ or e⁻. The third column shows the number density of C⁺ estimated from the LTE fits for comparison.

^a Estimated from an ionization calculation because no C ii data were available. See the text for details.

Download table as:  ASCIITypeset image

Figure 6 compares the residuals for the LTE and non-LTE fits. The non-LTE fits occasionally reproduce the data slightly better (for example for HD 32297), but overall the residuals are very similar. This suggests that dropping the LTE assumption is not enough to resolve the discrepancy between modeled and observed line fluxes seen for some disks (for example ¹²CO 2–1 toward HD 21997). One possible cause could be that our non-LTE model only considers excitation caused by collisions and cosmic microwave background radiation. However, stellar photons (or dust thermal emission) can sometimes be an important contribution to the excitation. For example, Matrà et al. (2018b) found that pumping to electronic and rovibrational levels by the stellar UV had an important effect on the population of rotational levels of CO in the disk around β Pic.

Figure 8 shows how the disks of our sample are distributed in the plane of C⁰ versus CO masses and vertical column densities. Looking at the masses (Figure 8, left), we note that for all disks where both masses are well constrained, M(CO) ≳ M(C⁰), with one notable exception: β Pic. Second, we can recognize two groups of disks: CO-rich disks with M(CO) > 10⁻³ M_⊕ and CO-poor disks with M(CO) < 10⁻⁴ M_⊕. Only one disk (HD 146897) might be located in the CO mass range between 10⁻⁴ and 10⁻³ M_⊕. In a secondary gas scenario, this might indicate that the transition from a CO-poor to a CO-rich disk occurs fast once C shielding becomes significant. We also note that there are no disks with high C⁰ mass and low CO mass (bottom right of the plot). In the secondary gas production scenario, one might expect such disks to exist. Once the CO production rate declines, the CO mass can decrease rapidly while C stays in the system (Figure 2 in Marino et al. 2020). However, our sample is biased and therefore we cannot infer the nonexistence of such systems. The analysis of a statistically complete sample will be presented in a forthcoming paper. In Section 6.1, we discuss in detail how our results compare to models of secondary gas production.

Zoom In Zoom Out Reset image size

Figure 8 (left) also shows masses for the three protoplanetary disks DM Tau, LkCa 15, and TW Hya with red data points. While there are more protoplanetary disks with C i data (Kama et al. 2016a, 2016b; Alarcón et al. 2022), to the best of our knowledge, these are the only disks with published C⁰ masses (Tsukagoshi et al. 2015). The CO masses are from Zhang et al. (2019, their Figure 4) for DM Tau and from Jin et al. (2019) for LkCa 15. For TW Hya, we consider two CO mass estimates: the smaller is derived by combining the minimum CO abundance (relative to H₂) from Favre et al. (2013) with the minimum disk mass from Bergin et al. (2013). The larger CO mass is from Zhang et al. (2019, their Figure 4). The error bar for TW Hya encompasses the range between these two mass estimates. The three protoplanetary disks considered here have a substantially higher CO mass compared to the CO-rich debris disks, although this might not be true in general (e.g., Smirnov-Pinchukov et al. 2022). On the other hand, the protoplanetary C⁰ masses are comparable to (or possibly lower than) the debris disk C⁰ masses. This probably reflects the typical chemical structure of protoplanetary disks where atomic C is only expected to be present in the disk atmosphere, while most of the mass is located in the molecular layer closer to the midplane (e.g., Kama et al. 2016a). We also note that all three protoplanetary disks considered here are found around low-mass stars, complicating the comparison to our sample, which consists mostly of A-type stars.

Looking at the vertical column densities (Figure 8, right), we see the same two groups of disks as for the masses. The vertical dashed line marks a C⁰ column density of 10¹⁷ cm⁻², corresponding to a reduction of the CO photodissociation rate by roughly a factor of 1/e (Heays et al. 2017), assuming that C⁰ forms a distinct layer above and below CO (if they are well mixed instead, the shielding efficiency is reduced, Marino et al. 2022). We see that the CO-rich disks are consistent with being significantly shielded by neutral carbon, as required by the secondary gas scenario. The CO-poor disks on the other hand are not C⁰ shielded. The horizontal dashed line in Figure 8 (right) marks a CO column density of 10¹⁷ cm⁻² where self-shielding reduces the ISRF-induced dissociation rate to 1.6% of the unshielded rate (Visser et al. 2009). It shows that CO-rich disks are self-shielded, while CO-poor disks are not (at a CO column density of 10¹⁵ cm⁻², the dissociation rate is still at 40% of the unshielded value, Visser et al. 2009). We emphasize that from a secondary gas perspective, CO self-shielding is not sufficient to explain CO-rich disks. Indeed, unless the CO production rate is very large, CO can never reach a self-shielding column density without the help of C⁰ shielding (e.g., Cataldi et al. 2020).

We test whether the carbon ionization fraction can be constrained by our observations. In Table 9, we compute two ionization fractions: χ_ana is purely based on our derived C⁰ masses and corresponds to the analytically calculated ionization fraction of a gas parcel that is subject to the full stellar and interstellar radiation (that is, we neglect any optical depth). Other than that, the calculation was performed as described in Section 5.1 (ionization equilibrium with n(e⁻) = n(C⁺)). χ_ana should give a reasonable estimate for the disks with low C⁰ column densities (HD 48370, HD 61005, HD 95086, HD 146897, and HD 181327), but will overestimate the ionization for other disks. We also calculated χ_obs for the targets with C ii data. It is simply χ_obs = M(C⁺)/(M(C⁰) + M(C⁺)). Unfortunately, we are unable to provide strong constraints for χ_obs. The strongest constraint is found for the disk around β Pic where χ_obs > 0.4.

In this section, we compare our results to the model by Marino et al. (2020) that considers the production and evolution of secondary gas in debris disks. The Marino et al. (2020) model is one dimensional and assumes that CO is produced in a collisional cascade of planetesimals that releases CO ice into the gas phase. The radial variation of the CO production rate is modeled as a Gaussian centered at r_belt. The CO production rate is proportional to the rate at which mass is input by catastrophic collisions of the largest bodies in the cascade. This rate decreases with time as the disk depletes its mass. The model considers CO photodissociation as well as CO self-shielding and shielding by C⁰ and follows the viscous radial evolution of CO and C⁰ assuming a standard α-disk model. C⁰ can only be removed by accretion onto the star, while CO can be removed by both accretion and photodissociation. CO shielding by C⁰ is assumed to be maximally efficient, that is, C⁰ is assumed to be in a layer above and below CO, rather than well mixed (see Marino et al. 2022, for a discussion of how vertical mixing affects the CO lifetime).

Marino et al. (2020) constructed two synthetic populations of A-type star debris disks with gas: one assuming a low viscosity (α = 10⁻³) and one with a high viscosity (α = 0.1). To construct each synthetic population, they evolve a sample of model systems with randomly chosen stellar masses, disk masses, and disk radii. The systems are evolved to a random age between 3 and 100 Myr. By comparing their synthetic populations to CO observations, they conclude that a high viscosity (α = 0.1) is required. Indeed, the low viscosity population produces too many disks with very large (≳1 M_⊕) CO masses, inconsistent with observations. This is because for low viscosity, the removal of gas by accretion is slow, allowing a large column density to accumulate, which leads to efficient shielding.

Marino et al. (2020) also compared their model population to observations of C⁰, but they had only four disks with C⁰ data available: 49 Ceti, β Pic, HD 32297, and HD 131835. Considering the low number statistics, they concluded that the α = 0.1 model is in reasonable agreement with the C⁰ data.

We are now in a position to provide an updated comparison of the Marino et al. (2020) simulations with observations: we now have 10 disks around A-type stars with C⁰ data. Furthermore, new CO isotopologue data allow more accurate CO mass measurements compared to the values used by Marino et al. (2020). In particular, our CO masses for 49 Ceti and HD 32297 are larger by factors of ∼50 and ∼40, respectively, compared to the values adopted by Marino et al. (2020). Figure 9 shows the comparison of our derived CO and C⁰ masses to three population-synthesis models: the original models presented by Marino et al. (2020) with α = 10⁻³ and 10⁻¹, and a new model with α = 10⁻². One complication is that the total C⁰ mass of the simulations is often dominated by low-density gas at large radii that is not easily observable, or that might be ionized in reality (the model ignores ionization). To avoid this issue, the simulated C⁰ masses shown in Figure 9 correspond to what Marino et al. (2020) called the observable carbon mass defined by ${{\rm{\Sigma }}}_{{{\rm{C}}}^{0}}({r}_{\mathrm{belt}}){r}_{\mathrm{belt}}^{2}\pi$ with ${{\rm{\Sigma }}}_{{{\rm{C}}}^{0}}$ the C⁰ surface density and r_belt the mid-radius of the gas-producing belt.

Zoom In Zoom Out Reset image size

In Figure 9, each small dot (colored by the dust fractional luminosity) represents a single simulation. To take into account observational biases, Marino et al. (2020) consider the subpopulation of simulated systems that corresponds to the sample proposed by Moór et al. (2017), which is defined by the following criteria: (1) A-type host star, (2) fractional dust luminosity between 5 × 10⁻⁴ and 10⁻², (3) a dust temperature below 140 K, (4) a detection with Spitzer or Herschel at λ ≥ 70 μm, and (5) an age between 10 and 50 Myr. The black contours enclose 68%, 95%, and 99.7% of that subpopulation, but considering model and observational uncertainties, we shall consider the outermost contour as a reference to compare models and observations. We have 10 A-type stars in our sample, all of which are part of the Moór et al. (2017) sample. They are also shown in Figure 9 and should fall within the black contours if the model and data agree. We confirm the finding by Marino et al. (2020) that the low viscosity (α = 10⁻³) model overpredicts the observed masses. However, in contrast to Marino et al. (2020), we find that the high viscosity model (α = 10⁻¹) also struggles to reproduce the observed masses. In particular, the model overpredicts the C⁰ masses. The case with α = 10⁻² turns out to be qualitatively similar to α = 10⁻³.

As mentioned before, comparing observed and simulated C⁰ masses is not trivial. Instead, comparing column densities is more promising because it is the column densities of CO and C⁰ that determine the CO photodissociation timescale (Marino et al. 2020). Thus, in Figure 10, we instead compare the synthetic population and the observations in terms of vertical column density. More precisely, we compare the model column density at r_belt (i.e., the radius where CO production peaks) to the observed column density. Again comparing the observations to the black contours, the conclusion remains essentially unchanged: the models tend to overpredict the C⁰ column density.

Zoom In Zoom Out Reset image size

One possible explanation for the mismatch between the data and the secondary gas model is that the gas is primordial. We will discuss a primordial scenario in Section 6.2. Here, we consider possible explanations for the mismatch in the framework of the secondary gas scenario.

6.1.1. Additional Shielding Agent

6.1.2. Increased CO Content of Gas-producing Comets

6.1.3. Radiation Pressure

6.1.4. High C Ionization Fraction

6.1.5. Underestimation of Observed C⁰ Masses

6.1.6. Transient Events

6.1.7. Removal of C by Planets or Grain Surface Chemistry

One possible way to distinguish primordial, H-rich gas from secondary, H-poor gas is to study the gas chemical composition. However, chemical modeling by Smirnov-Pinchukov et al. (2022) showed that molecular emission (besides CO) is expected to be undetectable regardless of the hydrogen content of the gas (with the possible exception of HCO⁺). Therefore, we instead consider how the CO fraction (that is, the fraction of the total carbon that is in CO) can inform us about the chemistry and hydrogen content of debris disk gas. To this end, we compare our results to the thermochemical model by Iwasaki et al. (2023). They consider only A-type stars, so we limit the discussion to the A-type stars in our sample. Iwasaki et al. (2023) compute the chemistry of debris disk gas using a photon-dominated region (PDR) code (the Meudon code, Le Petit et al. 2006). The model solves for the thermal and chemical equilibrium of a stationary, semi-infinite, plane-parallel slab of gas illuminated from one side. It considers various heating and cooling processes, gas phase, and grain surface chemistry as well as photoionization, photodissociation, and self-shielding, taking into account the radiative transfer. In this model, there is no gas production from comets or gas loss from accretion onto the star, and the dust mass (and dust surface) is constant. Thus, the setup is reminiscent of a primordial gas origin. Iwasaki et al. (2023) fit an analytical formula to the CO fraction computed with their numerical model. The formula is given by

$$\begin{eqnarray}&&\displaystyle \frac{n(\mathrm{CO})}{{n}_{{\rm{C}}}}={\left(1+{\left({10}^{-14}{\eta }^{1.8}+6\times {10}^{-11}\eta \right)}^{-1}\right)}^{-1}.\end{eqnarray} \tag{ 10 }$$

Here n(CO) is the CO number density, while n_C is the number density of carbon nuclei (in practice n_C ≈ n(CO) + n(C⁰) + n(C⁺)). The parameter η is given by

$$\begin{eqnarray}&&\eta ={n}_{{\rm{H}}}{Z}^{0.4}{\chi }^{-1.1}.\end{eqnarray} \tag{ 11 }$$

Here n_H is the number density of H nuclei and Z is the gas metallicity, where Z = 1 corresponds to solar metallicity. For a given metallicity Z, we can estimate n_H using our observational results as

$$\begin{eqnarray}&&{n}_{{\rm{H}}}\approx \displaystyle \frac{n(\mathrm{CO})+n({{\rm{C}}}^{0})}{Z{\zeta }_{{\rm{C}}}},\end{eqnarray} \tag{ 12 }$$

where ζ_C = 1.32 × 10⁻⁴ is the relative elemental abundance of C adopted by Iwasaki et al. (2023). We neglect n(C⁺) in our estimation because it is generally poorly constrained by our fits. The parameter χ describes the strength of the UV field. It is given by

$$\begin{eqnarray}&&\chi =\sqrt{{\chi }_{\mathrm{CO}}{\chi }_{\mathrm{OH}}},\end{eqnarray} \tag{ 13 }$$

where χ_CO and χ_OH are the normalized fluxes of UV photons in the wavelength ranges that are important for the chemistry of CO (91.2–110 nm) and OH (160–170 nm). The normalization constants are 1.2 × 10⁷ cm⁻² s⁻¹ (the Habing field) for χ_CO and 10¹² cm⁻² s⁻¹ for χ_OH. Splitting out the contribution from the ISRF and the star gives

$$\begin{eqnarray}&&{\chi }_{\mathrm{CO}}=\displaystyle \frac{{\chi }_{\mathrm{CO},\mathrm{ISRF}}}{2}+{\chi }_{\mathrm{CO},\mathrm{star}}({r}_{\mathrm{mean}}).\end{eqnarray} \tag{ 14 }$$

The factor 1/2 for χ_CO,ISRF accounts for the fact that the ISRF penetrates the modeled gas slab only from one side. Iwasaki et al. (2023) find χ_CO,ISRF = 1.3. To compute χ_CO,star, we use ATLAS9 stellar atmosphere models (Castelli & Kurucz 2003), with the exception of β Pic, which observationally shows additional emission in the UV above the predictions from a standard stellar atmosphere model. Therefore, we use a PHOENIX model as described in Fernández et al. (2006) complemented with UV data from the Hubble Space Telescope (Roberge et al. 2000) and the Far Ultraviolet Spectroscopic Explorer (Bouret et al. 2002; Roberge et al. 2006). χ_OH,star is computed in the same way as χ_CO,star, while χ_OH,ISRF is negligible. In Table 10 we show the computed values of χ_CO and χ_OH.

Table 10. The Normalized Flux of UV Photons Affecting the CO (χ_CO) and OH (χ_OH) Chemistry of Debris Disk Gas, Evaluated at the Midplane at r = r_mean

Star	χCO	${\chi }_{\mathrm{CO}}^{\mathrm{ext}}$	χOH
(1)	(2)	(3)	(4)
49 Ceti	3.0	5.1 × 10−3	1.7
β Pictoris	3.3	5.6 × 10−1	8.2 × 10−2
HD 21997	8.9 × 10−1	4.9 × 10−3	4.5 × 10−1
HD 32297	7.2 × 10−1	3.1 × 10−3	1.6 × 10−1
HD 95086	6.5 × 10−1	5.7 × 10−1	1.4 × 10−2
HD 110058	2.2	6.5 × 10−5	4.6
HD 121191	1.2	1.0 × 10−3	1.4
HD 121617	1.0 × 101	8.3 × 10−4	3.4
HD 131835	8.0 × 10−1	9.8 × 10−3	4.0 × 10−1
HR 4796	2.1 × 102	8.2 × 101	9.4

Note. (2) Normalized UV flux between 91.2 and 110 nm without extinction. (3) Same as (2), but with extinction by C⁰ and CO. (4) Normalized UV flux between 160 and 170 nm.

Download table as:  ASCIITypeset image

To account for UV extinction by C⁰ and CO, we follow Iwasaki et al. (2023) and multiply χ_CO,ISRF and χ_CO,star by a shielding factor. For χ_CO,ISRF, this shielding factor is given by

$$\begin{eqnarray}&&{f}_{{\rm{shield}}}^{{ISRF}}={e}^{-{\alpha }_{{\rm{C}}}{N}_{\perp }({{\rm{C}}}^{0})/2}{\left(1+{\left(\displaystyle \frac{{N}_{\perp }({\rm{CO}})/2}{{10}^{14}{{\rm{cm}}}^{-2}}\right)}^{0.6}\right)}^{-1},\end{eqnarray} \tag{ 15 }$$

where N_⊥ is the column density perpendicular to the disk midplane and α_C is the carbon ionization cross section. For the latter, we adopt the same value as Iwasaki et al. (2023): α_C = 1.777 × 10⁻¹⁷ cm² (Heays et al. 2017). The first term in Equation (15) represents extinction by C⁰, while the second term is a function fitted by Iwasaki et al. (2023) to the CO shielding factors tabulated in Visser et al. (2009). For the shielding factor for χ_CO,star, we replace N_⊥/2 by N_∥(r_mean) in Equation 15, that is, the column density parallel to the disk midplane from the star to r_mean. We note that χ_OH is not affected by extinction. The values for ${\chi }_{{\rm{CO}}}$ with extinction taken into account (denoted ${\chi }_{{\rm{CO}}}^{{\rm{ext}}}$) are shown in Table 10.

We then compute the predicted CO fraction for all A-type stars in our sample using Equation (10) for metallicities ranging from 1 (primordial gas) to 10³ (high metallicity as expected for secondary gas). We emphasize that the high metallicity case does not directly correspond to the secondary gas models that were discussed in the previous section, because there is no CO production from comets included in the Iwasaki et al. (2023) model discussed here. Rather, the high metallicity case would correspond to an equilibrium state that secondary gas would attain in the absence of gas production from comets, for example, if secondary gas was produced in a single burst rather than continuously.

The predicted CO fractions are shown with the blue lines in Figure 11. We also show the CO fraction estimated from our observations with black lines, where we again made the approximation n_C ≈ n(CO) + n(C⁰). Note that for HD 95086 and HR 4796, neither CO nor C i emission is detected and thus the CO fraction remains observationally unconstrained.

Zoom In Zoom Out Reset image size

We see that among the eight targets with an observationally constrained CO fraction, seven are consistent with the predictions from the thermochemical models. For most of these objects, both low and high metallicity models are consistent with the data. The only object clearly incompatible with the model predictions is β Pic where even for a solar metallicity gas, the gas density is too low to sustain the observed CO fraction purely by chemistry. This is consistent with the view that the gas in the disk around β Pic is currently produced from the destruction of cometary material (e.g., Dent et al. 2014; Matrà et al. 2017a; Iwasaki et al. 2023).

In summary, we find that the CO fraction of the CO-rich debris disks around 49 Ceti, HD 21997, HD 32297, HD 110058, HD 121191, HD 121617, and HD 131835 is consistent with the chemistry expected for a primordial gas. However, the analytical formula we applied assumes a simplified geometry (plane-parallel slab) and a simplified chemical network. Therefore, the uncertainties remain high and we do not interpret our results as a confirmation of the primordial gas origin. Rather, our results encourage further investigation of the primordial scenario. The logical next step would be to apply the full numerical model including a more realistic geometry to each of the disks presented here. Iwasaki et al. (2023) already present such detailed models for β Pic and 49 Ceti, reaching conclusions similar to ours: the disk around β Pic cannot be explained by steady-state chemistry, but the model can explain the disk around 49 Ceti if Z = 1–10.

If the gas is produced from the destruction of solid material in a steady-state collisional cascade, the current gas production rate can be computed using the equations presented by Matrà et al. (2017b). The steady-state condition implies that the rate at which mass is input to the cascade by catastrophic collisions of the largest bodies ${\dot{M}}_{{D}_{\max }}$ equals the sum of the rate at which CO and CO₂ are outgassed ${\dot{M}}_{\mathrm{CO}+{\mathrm{CO}}_{2}}$ and the rate at which mass is lost via radiation pressure on the smallest grains ${\dot{M}}_{{D}_{\min }}$:

$$\begin{eqnarray}&&{\dot{M}}_{{D}_{\max }}={\dot{M}}_{\mathrm{CO}+{\mathrm{CO}}_{2}}+{\dot{M}}_{{D}_{\min }}.\end{eqnarray} \tag{ 16 }$$

Assuming that CO and CO₂ are outgassed before reaching the smallest grain size in the system, we have

$$\begin{eqnarray}&&{\dot{M}}_{\mathrm{CO}+{\mathrm{CO}}_{2}}={f}_{\mathrm{CO}+{\mathrm{CO}}_{2}}{\dot{M}}_{{D}_{\max }},\end{eqnarray} \tag{ 17 }$$

where ${f}_{\mathrm{CO}+{\mathrm{CO}}_{2}}$ is the CO+CO₂ fraction (by mass) of the colliding bodies. Combining these two equations yields

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

However, Matrà et al. (2017b) show that ${\dot{M}}_{{D}_{\min }}$ can be estimated from observable quantities as follows:

$$\begin{eqnarray}&&{\dot{M}}_{{D}_{\min }}=1.2\times {10}^{3}{R}^{1.5}{\rm{\Delta }}{R}^{-1}{f}^{2}{L}_{* }{M}_{* }^{-0.5},\end{eqnarray} \tag{ 19 }$$

where R (in au) is the belt radius, ΔR (in au) the belt width, f the fractional luminosity, L_* (in L_⊙) the stellar luminosity, M_* (in M_⊙) the stellar mass, and ${\dot{M}}_{{D}_{\min }}$ is in M_⊕ Myr⁻¹.

Naively, we might expect that the current gas production rate ${\dot{M}}_{\mathrm{CO}+{\mathrm{CO}}_{2}}$ should correlate with the observed gas mass of a disk. However, this is not the case, as can be seen in Figure 12 showing the observed CO and CO+C⁰ masses as a function of ${\dot{M}}_{\mathrm{CO}+{\mathrm{CO}}_{2}}$. We follow Matrà et al. (2019b) and consider a range of possible CO+CO₂ mass fractions ${f}_{\mathrm{CO}+{\mathrm{CO}}_{2}}$ between 0.8% and 80%. Figure 12 demonstrates that ${\dot{M}}_{\mathrm{CO}+{\mathrm{CO}}_{2}}$ is not a good predictor for gas-rich debris disks. Figure 12 (left) also shows that stellar age is not the dividing factor between gas-rich and gas-poor systems. On the other hand, the host stars of gas-rich disks tend to have higher effective temperatures in our sample. This reflects the fact that gas-rich disks have so far only been found around A-type stars, but a more detailed interpretation is complicated because of the biases in our sample.

Zoom In Zoom Out Reset image size

As pointed out by Marino et al. (2020), it is theoretically expected that the gas viscous evolution is slower than the steady-state collisional evolution of the dust. Thus, the gas can retain a memory of previous states, and therefore the gas content depends on the time-integrated gas production rate of the disk rather than the current rate. Furthermore, a steady-state collisional cascade has the interesting property that the mass at late times (i.e., at times exceeding the collisional timescale of the largest planetesimal in the cascade) is independent of the initial disk mass because massive disks process their mass faster (Wyatt et al. 2007). In other words, two disks with the same current mass-loss rate (and therefore the same ${\dot{M}}_{\mathrm{CO}+{\mathrm{CO}}_{2}}$) can potentially have strong differences in their mass-loss rate history, and thus gas content.

Kral et al. (2023) discuss the conditions that determine whether debris disk gas is in Keplerian rotation or in an outflowing wind. They find that for L_* > 20 L_⊙, a wind driven by the stellar radiation is expected. For lower luminosities, the gas is in Keplerian rotation as long as the gas density is larger than a critical density n_crit. For n < n_crit, a wind driven by the stellar wind is expected. The critical density is given by

$$\begin{eqnarray}&&{n}_{\mathrm{crit}}\approx 7\,{\mathrm{cm}}^{-3}{\left(\displaystyle \frac{{\rm{\Delta }}r}{50\,\mathrm{au}}\right)}^{-1}{\left(\displaystyle \frac{{\alpha }_{{\rm{X}}}}{{\alpha }_{\mathrm{CO}}}\right)}^{-2/3},\end{eqnarray} \tag{ 20 }$$

where Δr = r_out − r_in and α_X is the polarizability of species X.

Among the stars we consider here, only HR 4796 has a luminosity larger than 20 L_⊙. We thus expect a radiative wind in this system. Gas removal with the wind might be the reason why we do neither detect CO nor C i emission from this system. For the other disks, we compared the critical density and the observationally estimated density for CO and C⁰. We classify a disk as being in the wind regime if both CO and C⁰ have densities smaller than n_crit. This is the case for the disks around HD 48370, HD 61005, and HD 95086, all of which remain undetected in both CO and C i. Furthermore, the data allow either regime for the disks around HD 146897 and HD 181327. Despite assuming that the gas is in Keplerian rotation (Section 4), our model still gives the correct total line flux for a given gas mass even in the wind regime because the emission for low gas density disks is optically thin. However, our flux measurements are affected by our assumption of a Keplerian disk. Indeed, line emission from a wind would come from a different spatial region, although the solid angle of that region should be roughly comparable to the Keplerian case. Furthermore, wind emission should extend over a significantly larger velocity range (Kral et al. 2023). We may roughly estimate the factor by which our upper limits on the line emission (and therefore gas mass) should be increased if gas emission arises from a wind. This factor is given by ${f}_{\mathrm{wind}}=\sqrt{{\rm{\Delta }}{v}_{\mathrm{wind}}/{\rm{\Delta }}v}$ where Δv is the size of the velocity range over which the integrated flux was measured and Δv_wind is the velocity range over which wind emission arises. Using Δv_wind = 20 km s⁻¹ (Kral et al. 2023), we find f_wind ≤ 2.3 for the three disks in the wind regime mentioned above. Such a small increase in the upper limits is inconsequential for our analysis.

To facilitate a uniform and simple analysis, we have integrated over the spatial and spectral dimensions of our data, that is, we compared our model to a single number extracted from the data: the disk-integrated flux. Our approach assumed that all the tracers we consider (CO, C⁰, and C⁺) are co-spatial and well mixed. As far as we can tell from the available images, co-spatiality is a reasonable assumption to derive the total gas masses. Indeed, the CO and C i emission morphologies look generally rather similar. For example, the emission of both CO and C i appear centrally peaked and with a similar radial extent for HD 21997 (Figure 1). Similarly, both CO and C i show a ring morphology for HD 121617 (see Figure 3 for C i and Figure 1 in Moór et al. 2017 for CO).

However, the fact that a significant fraction of our targets is spatially (and spectrally) resolved clearly suggests a path for future work that will model the data in its entirety (such as was performed by, for example, Kral et al. 2019; Cataldi et al. 2020) to compare the spatial distribution of the different gas tracers and the dust continuum. For example, to test the secondary gas production scenario, one could investigate whether the gas has viscously spread with respect to its production location traced by the dust continuum. Similarly, C⁰ could be more spread out than CO, since the latter can only exist in regions with a high C⁰ column density (e.g., Marino et al. 2020). Comparing the morphologies of different disks could also be instructive. For example, the ring morphology of the C i emission toward HD 121617 (Figure 3) seems inconsistent with a disk that has viscously spread all the way to the star (accretion disk), although detailed modeling is required to exclude the possibility of an accretion disk that is ionized and thus not observable in C i in the inner region. Cataldi et al. (2020) also concluded that there is no accretion disk morphology for the edge-on HD 32297 disk. In that latter case, the velocity information was used to constrain the inner disk regions. Several possible mechanisms could prevent the gas from accreting, for example, planets (Marino et al. 2020). On the other hand, the disk around HD 21997 is centrally peaked in both CO and C i, while the dust continuum shows a ring morphology (see Figure 1 and Kóspál et al. 2013; Moór et al. 2013b). This disk might be consistent with the picture where gas is produced in a planetesimal belt and then accretes onto the star. Finally, we note that disks produced from transient events (e.g., giant collisions) are expected to show asymmetries (e.g., Jackson et al. 2014; Kral et al. 2015; Jones et al. 2023), providing another motivation to study the spatial distribution of the gas and dust.

For the targets with C ii data, we ran additional LTE fits where the temperatures of C⁰ and C⁺ are separated. In other words, instead of two temperatures, we now consider three temperatures for CO, C⁰, and C⁺, where we impose ${T}_{\mathrm{CO}}\lt {T}_{{{\rm{C}}}^{0}}\lt {T}_{{{\rm{C}}}^{+}}$. We find that, with the exception of ${T}_{{{\rm{C}}}^{0}}$ for 49 Ceti, the C⁰ and C⁺ temperatures are not well constrained. Compared to the fits with two temperatures, we do not find any significant differences in the derived gas masses.

For our standard LTE fits, we assumed a square line profile with a fixed width of 0.5 km s⁻¹. We performed additional MCMC runs where we changed the line width. When decreasing the line width to 0.1 km s⁻¹, the derived CO masses change significantly (i.e., beyond their error bars) for HD 21997 (decrease by a factor of ∼2), HD 121617 (increase by a factor of ∼2) and HD 131835 (decrease by a factor of ∼2). These models reproduce the data similarly well as the standard fits. For the other disks, the masses remain within the error bars.

When instead increasing the line width to 2 km s⁻¹, we observe the following significant changes. For 49 Ceti, the C⁰ mass decreases by a factor of ∼3, but the model provides a significantly worse fit to the ¹³C i emission. For HD 21997, the CO mass decreases by two orders of magnitude, but the model provides a much worse fit to the ¹³CO and C¹⁸O data. For HD 121617, the CO mass decreases by a factor of ∼1.5 while the fit quality remains unchanged. For HD 131835, the CO mass decreases by about two orders of magnitude, but the model provides a much worse fit to the ¹³CO and C¹⁸O data. For the other disks, no significant changes are observed. In conclusion, these tests confirm that a line width of 0.5 km s⁻¹ is a reasonable choice.

For our standard LTE fits, we assumed flat priors for the CO and C temperatures between 10 and 200 K (Table 3). Here we perform additional MCMC runs where we changed the limits of the temperature priors. If we increase the lower limit to 20 K, the derived CO and C⁰ masses do not change significantly, except for β Pic, HD 21997, and HD 131835 for which the standard fit indicates a CO temperature below 20 K. However, for these three disks (as well as 49 Ceti which also has a CO temperature below 20 K in the standard fit), models with T_CO ≥ 20 K provide a significantly worse fit to the observed fluxes. This result strengthens our choice of 10 K as the lowest temperature.

When decreasing the upper limit of the temperature priors from 200 to 100 K, the derived CO and C⁰ masses do not change significantly.

The edge-on, gas-rich debris disk around 49 Ceti has been the subject of multiple detailed studies over the past years (some of the more recent work includes Higuchi et al. 2019a, 2019b, 2020; Moór et al. 2019; Pawellek et al. 2019; Klusmeyer et al. 2021). Hughes et al. (2017) presented two non-LTE models of the CO disk around 49 Ceti with CO/H₂ abundances of 10⁻⁴ (representing a primordial gas) and 1 (representing a secondary gas). They used observations of ¹²CO 2–1 and 3–2 from the SMA and ALMA, respectively. These models yield CO masses of 10^−3.8 and 10^−3.5 M_⊕, respectively, more than an order of magnitude smaller than our CO mass (10^−2.15 M_⊕). A similar discrepancy applies to the column density derived by Nhung et al. (2017) based on ¹²CO 3–2. We find that the ¹²CO 2–1 and 3–2 emission is strongly optically thick (τ ∼ 120), while Hughes et al. (2017) concluded that it is optically thin. A key difference is that we have CO isotopologue data available, while Hughes et al. (2017) were relying on ¹²CO data only. To test the influence of the CO isotopologue data, we ran additional non-LTE fits (with H₂ collisions) considering only our ¹²CO data. In that case, we are only able to derive a lower limit on the CO mass (>10^−3.4 M_⊕) because of the high ¹²CO optical depth.

On the other hand, when including the isotopologue data, our results are roughly consistent with the mass determined by Moór et al. (2019, 10^−2.0 M_⊕), but an order of magnitude lower than the lowest value given by Higuchi et al. (2020, 10^−1.2 M_⊕). Both of these works included CO isotopologues. This highlights the importance of CO isotopologue observations in determining accurate CO masses: the ¹³CO and C¹⁸O masses can be determined accurately thanks to their optically thin emission. The total CO mass then follows from our assumption that the ¹²CO/¹³CO and ¹²CO/C¹⁸O ratios are equal to the C isotope ratios in the ISM. We note, however, that our work and the analysis by Moór et al. (2019) and Higuchi et al. (2020) are based on disk-integrated fluxes, while Hughes et al. (2017) modeled spatially and spectrally resolved data with a 3D disk model. Thus, part of the difference in the mass estimates might also be explained by model differences.

Our CO temperature appears well constrained at 15–17 K for both the LTE and non-LTE models, indicating a low gas temperature. This is consistent with the results by Higuchi et al. (2020, 8–11 K) and Hughes et al. (2017, 14 K at 100 au for their model that assumes CO/H₂ = 1). As discussed by Higuchi et al. (2020), such a low gas temperature might not be unreasonable, although more theoretical modeling considering the heating and cooling processes of the gas (e.g., Zagorovsky et al. 2010; Kral et al. 2016, 2019) is required. We note that the dust temperature has been measured to be 59 K (Holland et al. 2017), that is, above the CO freeze-out temperature of 20 K (e.g., Yamamoto 2017). Therefore, we do not expect CO freeze-out to occur even though the gas temperature is low. The model by Kral et al. (2019) showed that T_dust > T_gas can indeed occur in debris disks. Similarly, the thermochemical debris disk model (PDR model) by Iwasaki et al. (2023) also gives a gas temperature below the dust temperature.

For C⁰, the temperature of ∼24 K as well as the column density we derive are consistent with the results by Higuchi et al. (2019a). The fact that T_CO < T_C might suggest that C and CO occupy different layers in the disk. The C⁺ mass is not well constrained because of high optical depth at low temperature, as was already found by Roberge et al. (2013) who derived a lower limit consistent with our result.

The edge-on disk around β Pic is among the best-studied debris disks. It hosts two giant planets, β Pic b (Lagrange et al. 2009) and β Pic c (Lagrange et al. 2019). The system shows time-variable absorption features (e.g., Beust et al. 1998; Kiefer et al. 2014; Welsh & Montgomery 2016) and transit events (Zieba et al. 2019; Lecavelier des Etangs et al 2022; Pavlenko et al. 2022) attributed to exocomet activity. The debris disk can be traced out to ∼2000 au in scattered light (Janson et al. 2021), while the millimeter-sized dust grains are located in a ring centered at ∼100 au (Dent et al. 2014).

Gas emission is detected from CO (e.g., Dent et al. 2014; Matrà et al. 2017a), C i (Cataldi et al. 2018), C ii (Cataldi et al. 2014), O i (Brandeker et al. 2016) as well as various heavier elements such as Fe i, Na i, and Ca ii (Brandeker et al. 2004; Nilsson et al. 2012). Various atomic lines are also seen in absorption (e.g., Lagrange et al. 1998; Roberge et al. 2006), including nitrogen (Wilson et al. 2019) and hydrogen (Wilson et al. 2017), where the latter is hypothesized to come from the dissociation of water originating from evaporating exocomets.

The masses we derive are consistent with the results by Matrà et al. (2017a, CO), Cataldi et al. (2018, C⁰), and Cataldi et al. (2014, 2018, C⁺). The disk around β Pic is the only disk in our sample where we can be sure that the C⁰ mass ($({1.4}_{-0.3}^{+0.4})\times {10}^{-4}$ M_⊕) is clearly larger than the CO mass ($({4.8}_{-1.4}^{+1.2})\times {10}^{-5}$ M_⊕).

The C temperature is not well constrained by our fits. On the other hand, the CO temperature we derive under the LTE assumption (${11.8}_{-1.1}^{+2}$ K) is in general significantly lower than previous estimates from the literature. Indeed, the theoretical model by Zagorovsky et al. (2010) predicts a gas temperature of ∼60 K at 100 au. Kral et al. (2016) derived a temperature of 50 K (at 100 au) from a PDR model fitted to CO, C i, and C ii data. Matrà et al. (2017a) empirically derive ∼160 K (again at 100 au) from the CO scale height. However, our temperature is consistent with the excitation temperature derived by Matrà et al. (2017a) from the CO 3–2/2–1 line ratio. It is also roughly consistent with the rotational excitation temperature of 15.8 ± 0.6 K measured by Roberge et al. (2000) from CO absorption. Therefore, our low temperature likely indicates that CO is not in LTE (Matrà et al. 2017a), and therefore that this temperature does not correspond to the kinetic temperature.

The debris disk around HD 21997 is among the disks with the highest CO masses in our sample. The disk served as a prototype for the proposed class of hybrid disks, that is, disks where secondary dust and primordial gas coexist (Kóspál et al. 2013). ALMA observations of the CO emission were presented by Kóspál et al. (2013), who found that the gas and dust are not colocated: there is a dust-free inner gas disk. If the gas was secondary, then this would suggest viscous spreading of the gas with respect to the dust (Kral et al. 2019; Marino et al. 2020).

Compared to the CO analysis by Kóspál et al. (2013) and Higuchi et al. (2020) that included ¹²CO 2–1 and 3–2, ¹³CO 2–1 and 3–2, and C¹⁸O 2–1, we added observations of C¹⁸O 3–2 (as well as additional observations of ¹³CO 3–2). The disk around HD 21997 thus has the most CO lines observed in our sample. We also publish the first measurements of C i and C ii, and thus the masses of C⁰ and C⁺ are constrained for the first time in this paper.

Our derived CO masses of ${4.8}_{-0.9}^{+0.8}\times {10}^{-2}$ M_⊕ (LTE) and ${6.8}_{-1.9}^{+7}\times {10}^{-2}$ M_⊕ (non-LTE, H₂ colliders) are consistent with the ranges reported by Kóspál et al. (2013, 4–8 × 10⁻² M_⊕) and Higuchi et al. (2020, 5.5–85 × 10⁻² M_⊕). There is also agreement in the low temperature of the gas: we derive ∼10 K (the lowest temperature allowed in our fits), while Kóspál et al. (2013) and Higuchi et al. (2020) derive 6–9 K and 8–12 K, respectively.

The debris disk around HD 32297 has been studied extensively in scattered light (Bhowmik et al. 2019; Duchêne et al. 2020, to cite some recent work), in the far-IR with Herschel (Donaldson et al. 2013) and in the submillimeter/millimeter regime (e.g., MacGregor et al. 2018). We derive a CO mass of ${5.2}_{-1.5}^{+1.7}\times {10}^{-2}$ M_⊕, consistent with the value derived by Moór et al. (2019). This is the largest CO mass in our sample and about two orders of magnitude larger than the CO mass derived by MacGregor et al. (2018). As for 49 Ceti, the reason for this discrepancy is that our work and Moór et al. (2019) include CO isotopologue data, in contrast to MacGregor et al. (2018).

This disk is also detected in C i (Cataldi et al. 2020) and C ii (Donaldson et al. 2013). Our estimate of the C⁰ mass is consistent with the value derived by Cataldi et al. (2020) and roughly one order of magnitude smaller than the CO mass. Unfortunately, the C⁺ mass is not well constrained due to high optical depth at low temperature.

The debris disks around HD 48370, HD 61005, HD 95086, and HR 4796 are the disks in our sample that remain undetected in both CO and carbon. Our LTE upper limit on the CO mass in the HD 61005 disk is a factor of ∼6 higher than the value derived by Olofsson et al. (2016) assuming optically thin emission and LTE. Since Olofsson et al. (2016) did not publish the line flux upper limit they used for their estimate, it is difficult to find the reason for this difference. For HR 4796, our LTE upper limit on the CO mass is roughly consistent with the upper limit derived by Kennedy et al. (2018). Comparing instead to the LTE and non-LTE results by Kral et al. (2020b), our LTE upper limit is about an order of magnitude more constraining, while the higher non-LTE upper limit is roughly consistent with their derivation.

The derived upper limits show that these disks all have a CO content smaller than β Pic. The same is true for C⁰, except maybe for HD 95086 where the upper limit on the C⁰ mass is relatively high. The two G-stars HD 48370 and HD 61005 are the two coolest stars in our sample. Matrà et al. (2019a) suggested that for a fixed fractional luminosity, the production rate of secondary CO gas is proportional to the stellar luminosity, which might be the reason why we do not detect any gas. Moreover, these two stars are also among the oldest (∼40 Myr) in our sample, which might be another reason. HD 95086 does not stand out in any particular way in our sample in terms of stellar effective temperature, age, or fractional luminosity. On the other hand, HR 4796 has the highest effective temperature (spectral type A0, Houk 1982) and is the youngest system (8 Myr) in our sample. Kennedy et al. (2018) suggest that the non-detection of CO emission cannot be solely explained by the intense UV radiation field of an A0 star that limits the lifetime of CO molecule to 40 yr or less. In addition, the planetesimals in the HR 4796 system need to be CO-poor, with a CO+CO₂ ice mass fraction of <1.8%, which is smaller than solar system comets. It is indeed possible that the planetesimals we observe today at a radius of ∼80 au formed interior to the CO snowline, due to the high luminosity of HR 4796 (Marino et al. 2020). From Figure 1 in Matrà et al. (2018a), the CO snowline in the HR 4796 protoplanetary disk could have been located as far out as ∼100 au. On the other hand, if the CO emission is not in LTE, the CO upper limit becomes less constraining and the planetesimals could have a CO+CO₂ content similar to or even larger than solar system comets (Kral et al. 2020b). Another reason for the absence of gas emission from the HR 4796 system could be that gas is removed by radiation pressure, because the luminosity of HR 4796 is larger than 20 L_⊙ and therefore a wind driven by the stellar radiation can be expected (Kral et al. 2017, 2023).

The edge-on debris disk around HD 110058 was discovered by Kasper et al. (2015) in scattered light. The disk shows absorption features most probably due to hot circumstellar gas (Hales et al. 2017; Iglesias et al. 2018; Rebollido et al. 2018). Emission from ¹²CO and ¹³CO has been detected (Lieman-Sifry et al. 2016; Hales et al. 2022). Hales et al. (2022) derive a CO mass of $\mathrm{log}({M}_{\mathrm{CO}}\,[{M}_{\oplus }])=-{1.16}_{-0.77}^{+1.72}$ (95% credible interval), consistent with our $-{2.03}_{-0.12}^{+0.14}$. However, Hales et al. (2022) derive a significantly lower CO temperature (∼18 K at the radius of the debris belt, compared to our ∼100 K). They consider a power law for the radial temperature dependence, which is likely more realistic than our assumption of a constant temperature.

We present the first observations of C i toward this target. No C i emission was detected. This is an interesting result: among the disks with a high CO mass (see Figure 8), HD 110058 is the only one without C i detected. Hales et al. (2022) applied the Marino et al. (2020) secondary gas model to the debris disk around HD 110058 and successfully reproduced the observed CO mass. However, they predict a C⁰ mass of ∼10⁻² M_⊕, which exceeds our upper limit of M(C⁰) < 3.2 × 10⁻⁴ M_⊕. Thus, our data allow us to exclude their secondary gas model.

The debris disk around HD 121191 was first identified by Melis et al. (2013). The disk shows strong CO emission (Moór et al. 2017). The CO mass we derive is consistent with the masses published by Moór et al. (2017) and the LTE mass by Kral et al. (2020b), but our lower limit on the non-LTE CO mass is ∼30% larger than the largest non-LTE mass derived by Kral et al. (2020b). Our work presents the first detection of C i toward this target. We find that the C⁰ mass is about an order of magnitude smaller than the CO mass.

Dust emission from the debris disk around HD 121617 was first reported by Mannings & Barlow (1998), while CO emission was first detected by Moór et al. (2017), who estimated a CO mass roughly consistent with our results. Our work presents the first detection of C i emission. Similarly to HD 121191, the C⁰ mass is smaller than the CO mass by roughly an order of magnitude.

The debris disk around HD 131835, first detected by Moór et al. (2006), is known to show ¹²CO (Moór et al. 2015b; Lieman-Sifry et al. 2016; Hales et al. 2019) and C i (Kral et al. 2019) emission. For our model, we use additional observations of ¹³CO and C¹⁸O that will be presented in more detail in A. Moor (2023, in preparation). Compared to the CO mass published by Moór et al. (2017; ∼3.8 × 10⁻² M_⊕ after correcting for the updated distance of the target), we find a CO mass roughly a factor of 2 smaller (∼1.7 × 10⁻² M_⊕). On the other hand, our CO mass is larger by a factor of ∼4 compared to the CO mass derived by Hales et al. (2019) with a 2D thermochemical model (4.6 × 10⁻³ M_⊕), potentially because Hales et al. (2019) did not consider the CO isotopologue data and/or because they assumed that the dust temperature equals the gas temperature.

Our C⁰ mass is consistent with the value derived by Kral et al. (2019), although the mass is not well constrained due to optical depth and could be much higher if the gas temperature is low (see the corresponding corner plot).

The edge-on debris disk around HD 146897 has been resolved at both millimeter wavelengths with ALMA (Lieman-Sifry et al. 2016) as well as scattered light (e.g., Goebel et al. 2018). The disk shows CO emission (Lieman-Sifry et al. 2016) as well as absorption by hot circumstellar gas (Rebollido et al. 2018). In this work, we present ACA observations of C i that did not detect any emission. Thus, HD 146897 is another example of a disk with CO where no C i was detected, although it might be less surprising than in the case of HD 110058 because (1) the CO mass is likely much smaller (although it could be large if the temperature is low as can be seen in the corner plot) and (2) the ACA C i observations are much less sensitive.

HD 181327 is a member of the β Pic moving group and hosts a thin debris ring that exhibits weak CO emission (Marino et al. 2016). We derive a small CO mass of ∼2 × 10⁻⁶ M_⊕, consistent with the analysis by Marino et al. (2016).

Among the disks with a CO detection in our sample, this disk has by far the smallest CO mass (about 1.5 orders of magnitude smaller than the second smallest CO mass found in the β Pic debris disk). There are other debris disks with similarly small CO masses such as Fomalhaut (Matrà et al. 2017b) or TWA 7 (Matrà et al. 2019a), but they were not included in our sample due to a lack of C i data. In this work, we present the first observations of C i toward HD 181327. The line remains undetected, as does C ii (Riviere-Marichalar et al. 2014).