Physical conditions of gas components in debris disks of 49 Ceti and HD 21997

Characterization of gas component in debris disks is of fundamental importance for understanding its origin. Toward this goal, we have conducted non-LTE (local thermodynamic equilibrium) analyses of the rotational spectral lines of CO including those of rare isotopologues ($^{13}$CO and C$^{18}$O) observed toward the gaseous debris disks of 49 Ceti and HD 21997 with the Atacama Large Millimeter/submillimeter Array (ALMA) and Atacama Compact Array (ACA). The analyses have been carried out for a wide range of the H$_{2}$ density, and the observed line intensities are found to be reproduced, as far as the H$_{2}$ density is higher than 10$^{3}$ cm$^{-3}$. The CO column density and the gas temperature are evaluated to be (1.8-5.9)$\times$10$^{17}$ cm$^{-2}$ and 8-11 K for 49 Ceti and (2.6-15)$\times$10$^{17}$ cm$^{-2}$ and 8-12 K for HD 21997, respectively, where the H$_{2}$ collision is assumed for the rotational excitation of CO. The results do not change significantly even if electron collision is considered. Thus, CO molecules can be excited under environments of no H$_{2}$ or a small number of H$_{2}$ molecules, even where the collision with CO, C, O, and C$^{+}$ would make an important contribution for the CO excitation in addition to H$_{2}$. Meanwhile, our result does not rule out the case of abundant H$_{2}$ molecules. The low gas temperature observed in the debris disks is discussed in terms of inefficient heating by interstellar and stellar UV radiation.


INTRODUCTION
Gas component of debris disks has attracted many astronomers and planetary scientists in relation to evolutionary process from protoplanetary disks to planetary system. It has been detected in optical absorption lines of some atoms (e.g., Slettebak 1975;Roberge et al. 2000Roberge et al. , 2006, far-infrared emission lines of [O i] and [C ii] (e.g., Donaldson et al. 2013;Roberge et al. 2013Roberge et al. , 2014Riviere-Marichalar et al. 2012Cataldi et al. 2014;Brandeker et al. 2016), and millimeter/submillimeter emission lines of CO (Zuckerman et al. 1995;Dent et al. 2005;Moór et al. 2011Moór et al. , 2015. With the advent of the Atacama Large Millimeter/submillimeter Array (ALMA), the CO emissions have extensively been studied in many debris disks during recent years (e.g., Kóspál et al. 2013;Dent et al. 2014;Lieman-Sifry et al. 2016;Hughes et al. 2017;Moór et al. 2017;Matrà et al. 2017;Hughes et al. 2018;Moór et al. 2019). Not only 12 CO but also its isotopic species ( 13 CO, and C 18 O) have successfully been observed in a few sources, e.g., HD 21997, HD 121191, HD 121617, HD 131488, 49 Ceti, and HD 32297 (Kóspál et al. 2013;Moór et al. 2017Moór et al. , 2019. In addition, the [C i] emission has been detected around 49 Ceti (Higuchi et al. 2017(Higuchi et al. , 2019a, β Pictoris (Higuchi et al. 2017;Cataldi et al. 2018), HD 131835 , and HD 32297 (Cataldi et al. 2020). Very recently, detection of [ 13 C i] has been reported in 49 Ceti (Higuchi et al. 2019b). CO emissions in debris disks have often been assumed to be optically thin, because its intensities from several gaseous debris disks are faint even with ALMA observations. However, detection of 13 CO and C 18 O in a few sources Moór et al. 2017Moór et al. , 2019 clearly indicates that the 12 CO emission is no longer regarded as optically thin in these gas disks. This situation is also true for the [C i] emission (Higuchi et al. 2019a,b).
Based on the observed peak intensities of optically thick 12 CO and 13 CO lines, Kóspál et al. (2013) suggested excitation temperatures as low as 6-9 K in the disk of HD 21997. Assuming LTE (local thermodynamic equilibrium) they derived a total CO mass of about (4-8)×10 −2 M ⊕ . Since the rotational temperature is much lower than the dust temperature reported so far Holland et al. 2017), such a low excitation temperature raised an important issue on the physical condition and the amount of the gas component in debris disks. Later, similar trends of the low excitation temperature were reported for several debris disks (Flaherty et al. 2016;Moór et al. 2017;Di Folco et al. 2020). Since the temperature is derived under the LTE assumption, it was thought that the non-LTE effect may affect the excitation temperature (Matrà et al. 2015). For 49 Ceti, Hughes et al. (2017) analyzed the 12 CO(J=3-2) data observed with ALMA and the 12 CO(J=2-1) data observed with the Submillimeter Array (SMA) by using the non-LTE code along with the disk model, and derived the temperature structure of the form as: where R is the radial distance from the central star and T 100 is the gas temperature at R=100 au. T 100 is evaluated to be 40 K and 14 K for the H 2 /CO ratio of 10 4 and 1, respectively. They also reported that the CO lines in 49 Ceti are almost thermalized. Similar analyses have been reported for the other debris disks (Matrà et al. 2015;Hales et al. 2019). However, these studies make use of only the 12 CO data. For further constraint of the gas kinetic temperature, it is essential to involve the 13 CO and C 18 O data in the analysis.
Moreover, the CO abundance relative to H 2 is still uncertain, although it is deeply related to the origin of the gas (i.e., a primordial gas inherited from the protoplanetary disk and/or a secondary gas outgassing from dust grains and icy solids). Higuchi et al. (2017) pointed out, by using a simple chemical model, that the C/CO column density ratio (hereafter C/CO ratio) is sensitive to the H 2 number density. To make full use of this method, we need to derive the C/CO ratio accurately by eliminating the effect of the opacity of the CO and [C i] emission.
In order to constrain the gas kinetic temperature and the CO column density accurately, use of as many molecular lines as possible is essential. It is particularly important to include optically thin lines such as rare isotopologue lines. To demonstrate this, we here focus on 49 Ceti and HD 21997, which are nearby ( (Moór et al. 2006;Torres et al. 2008;Bell et al. 2015). There are several high-sensitivity ALMA archival data for these two sources listed in Table 1. In this study, we have reanalyzed these archival data sets of 49 Ceti and HD 21997 to derive the column density and the gas temperature for observational understandings of the gas properties in these sources apart from the disk modeling.

DATA REDUCTION 2.1. 49 Ceti
We used the ALMA archival data sets of 49 Ceti obtained by using the Band 6 and Band 7 receivers (Table 1). For a fair comparison among the observations with different spatial resolution, the synthesized beam size for the 12 m array data is adjusted to that for the ACA data (∼6 ′′ ) by using the CASA task imsmooth. The synthesized beam sizes are summarized in Table 1. The angular size of 6 ′′ corresponds to 340 au.
We used version 4.7.2 and 5.3.0 of the Common Astronomy Software Applications package (CASA) (McMullin et al. 2007) for calibration and imaging, respectively. The CASA task tclean was employed to Fourier-transform the visibility data and to deconvolve the dirty images at a velocity interval of 2 km s −1 . Briggs weighting of +0.5 was applied to all line images for the best compromise between resolution and sensitivity. Continuum subtraction was performed for all the data by using the CASA task uvcontsub.

HD 21997
For HD 21997, we used the ALMA archival data (reduced cube data) reported by Kóspál et al. (2013) (Table 1). The synthesized beam is adjusted to be 4. ′′ 0 × 4. ′′ 0 with P.A.= 0.0 • for all the CO lines to cover the most of the CO distribution. For this source, the angular size of 4 ′′ corresponds to 280 au. show the integrated intensity maps of the 12 CO(3-2), 12 CO(2-1), and 13 CO(2-1) emissions, respectively, at a 6 ′′ resolution, where the velocity range for integration is from −6 to 11.5 km s −1 . Figure 1 (d) shows the integrated intensity map of the C 18 O(2-1) emission at a 6 ′′ resolution, where the velocity range for integration is from −4 to 6 km s −1 . Figures 2 (a) -(d) show the averaged spectra of the 12 CO(3-2), 12 CO(2-1), 13 CO(2-1), and C 18 O(2-1) emissions, respectively. These spectra are prepared by averaging over the 6 ′′ area in diameter which corresponds to the spatial resolution. The velocity range for the integrated intensity maps mentioned above is justified by these spectra. Each spectrum reveals a double-peak profile. Except for the C 18 O(2-1) data, the intensity of the red-shifted component is brighter than that of the blue-shifted component. The C 18 O(2-1) line is so faint that the data quality is not as good as the other lines. Double-Gaussian fitting is performed on all the spectra except for C 18 O(2-1) to derive the line parameters. For C 18 O(2-1), only the blue shifted component is fitted by the single Gaussian function. The results are listed in Table  2 1 .

HD 21997
Figures 3 (a) -(e) show the integrated intensity maps of the 12 CO(3-2), 12 CO(2-1), 13 CO(3-2), 13 CO(2-1), and C 18 O(2-1) emissions, respectively, at a 4 ′′ resolution, where the velocity range for integration is from −4 to 6 km s −1 . Details of the dataset are described elsewhere . Figures 4 (a) -(e) show the averaged spectra of the 12 CO(3-2), 12 CO(2-1), 13 CO(3-2), 13 CO(2-1), and C 18 O(2-1) emissions, respectively. These spectra are prepared by averaging over the 4 ′′ area in diameter which corresponds to the spatial resolution. Since all the spectra show a double-peak emission, double-Gaussian fitting is performed on the spectra to derive the line parameters. The results are listed in Table 3. For the 13 CO(3-2), 13 CO(2-1) and C 18 O(2-1) emissions, there seem to be small residuals (e.g., Figures 4 (b) 0.03 K, (d) 0.02 K and (e) 0.03 K) around the systemic velocity, which cannot be reproduced in the double-Gaussian fitting. Since it is close to the noise level, its contribution to the total spectra is not large. Thus, we ignore this feature in the following analyses.

LTE analysis
According to the result of double-Gaussian fitting, the peak intensity ratio of 13 CO(2-1)/ 12 CO(2-1) is calculated to be 0.46±0.03 and 0.48±0.04 for 49 Ceti and HD 21997, respectively. This ratio clearly indicates that the 12 CO(2-1) emission is optically thick, if the 12 C/ 13 C ratio of 77 (Wilson & Rood 1994) is assumed. A large optical depth of the CO line was already reported by Kóspál et al. (2013) for HD 21997. It is now revealed in 49 Ceti for the first time.
For a more detailed treatment, the rotational temperature, the column density of 12 CO, and the beam filling factor are derived by using the following equation, where we take the effect of the optical depth into account: Here, T B is the observed brightness temperature, T rot is the rotational temperature which is equivalent with the excitation temperature in the LTE condition, B ν (T rot ) is the Planck function, T CMB is the temperature of the cosmic background radiation, and f d is the beam filling factor. The optical depth, τ , is written as below (e.g., Goldsmith & Langer 1999;Sakai et al. 2008): 1 The conversion equation from Jy beam −1 to K is given as : Jy beam −1 , where θ maj is the major axis, θ min is the minor axis of the synthesized beam, ν is the frequency, and S is the flux density. and where A ul is the Einstein A-coefficient of the transition between the states u and l, N is the column density along the line of sight, g u is the degeneracy of the upper state, E u is the upper state energy, ∆v is the velocity width, and U(T rot ) is the partition function. We employ the isotope ratios of 1/77 and 1/560 for 13 CO and C 18 O, respectively (Wilson & Rood 1994). The spatial distributions and the velocity structures of the 12 CO(3-2), 12 CO(2-1) and 13 CO(2-1) emissions are well correlated with one another (see Figures 1 -4), and hence, the rotational temperature, the column density, and the beam filling factor of CO for 49 Ceti and HD 21997 are evaluated by using the least-squares analysis on the observed intensities of CO and its isotopic species. We treat the data for the blue-shifted and red-shifted components separately. The fitting results (i.e., calculated intensities and residuals in the fit) are shown in Tables 2 and 3, for 49 Ceti and HD 21997, respectively, whereas the derived parameters (T rot , N, and f d ) are summarized in Tables 4 and 5, respectively. The observed peak intensities are well reproduced within the measurement uncertainties for both sources. The highest correlation among the three parameters in the fit is 0.99 (correlation coefficient between parameters of the least squares analysis) between T rot and f d . Nevertheless, these two parameters are separately determined in the fit, and the effect of the correlation is reflected in the quoted errors.
For 49 Ceti, the optical depths of the 12 CO(2-1) and 12 CO(3-2) emissions are evaluated to be 38-46 and 18-24, respectively (Table 4). These 12 CO lines are indeed optically thick. The excitation temperature is also derived to be ∼ 9 K from our LTE analysis. According to the high resolution image ( Figure 5) reported by Hughes et al. (2017), the FWHM (full width half maximum) size of the emitting region of the CO(3-2) lines is estimated to be 1. ′′ 9 × 0. ′′ 8 and 1. ′′ 8 × 0. ′′ 8 for the blueshifted and red-shifted components, respectively, by using 2D Gaussian fits. Hence, the beam filling factor is estimated to be 0.04 for both blue-shifted and red-shifted components. These values are almost consistent with those determined in the fit (0.06±0.01 and 0.08±0.04 for the blue-shifted and red-shifted components, respectively).
For HD 21997, the optical depths of the 12 CO(2-1) and 12 CO(3-2) emissions are 66 -76 and 26 -36, respectively (Table 5). The excitation temperature of gas component in HD 21997 is determined to be 8 K. Since the emitting region of the CO emission is derived to be 1. ′′ 1 × 0. ′′ 9 and 1. ′′ 0 × 0. ′′ 9 for the blue-shifted and red-shifted components, respectively, by using 2D Gaussian fits on the high resolution channel maps shown in Figure 3 of Kóspál et al. (2013), the beam filling factor is estimated to be 0.06. This value is consistent with those determined in the fit (0.08±0.06 and 0.07±0.05 for the blue-shifted and red-shifted components, respectively), although the observed beam filling factors suffer from a large error.
For both sources, the beam-averaged excitation temperature (i.e., rotational temperature) is confirmed to be lower than 10 K based on the multi-line LTE analysis. In general, the rotational temperature derived from the LTE analysis does not always correspond to the gas kinetic tempera-ture, if the rotational level population is not well thermalized. To derive the gas kinetic temperature, a non-LTE analysis is performed.

Non-LTE analysis
We employ the non-LTE code prepared by ourselves in the analysis, which is tested by using the results reported by Goldreich & Kwan (1974) and also those derived from the RADEX code. The collisional cross sections are taken from The Leiden Atomic and Molecular Database (LAMDA) (Schöier et al. 2005). In this case, we need to consider the major collision partner for the rotational excitation of CO. If the gas is mainly primordial (i.e., remnant of the protoplanetary disk), collision with H 2 , the way most abundant constituent of such gas material, would be dominant. On the other hand, if the gas has an almost secondary origin (i.e., releases from icy grains and planetesimals), collision with H and O produced by the photodissociation of H 2 O, the C atom, the C + ion, electron, and CO itself can be considered. Collisional rate coefficients for H are lower by an order of magnitude than that for H 2 (Yang et al. 2013;Walker et al. 2015), while those for an electron are higher by two orders of magnitude than that for H 2 (Dickinson et al. 1977). On the other hand, collisional rate coefficients for C, O, CO and C + are not available, but they are roughly assumed to be comparable to those for H 2 . Based on these considerations, we regard an electron as the major collisional partner for the secondary origin case. Thus, we conduct the non-LTE analysis for the two distinct cases, where the major collision partner is H 2 or an electron. Caveats for this simplification are discussed later.

The case of H 2 collision
First, we conduct the non-LTE analysis based on the LVG (large velocity gradient) model (Goldreich & Kwan 1974;van der Tak et al. 2007), assuming the collision with H 2 . In this analysis, the rotational level populations are determined by the balance between the collision and radiation processes, from which the line intensities are calculated. Hence, there are four parameters to be determined: the H 2 density, the column density of CO, the gas kinetic temperature, and the beam filling factor. We have initially tried to optimize these four parameters by using the least-squares fit to reproduce the 12 CO, 13 CO and C 18 O data. However, all the four parameters cannot be determined simultaneously in the fit. Specifically, the H 2 density is not well constrained. Hence, we fix the H 2 density and determine the remaining three parameters. The range of the H 2 density assumed in the fit is from 3×10 3 to 1×10 7 cm −3 . The data for the blue-shifted and red-shifted components are treated separately. The least-squares fit does not converge for any of the two sources for the H 2 density below 10 3 cm −3 . Thus, we set the lowest H 2 density to be 3×10 3 cm −3 , for which the fit is successful. The fitting results (i.e., calculated intensities and residuals in the fit) are shown in Tables 2 and 3, and derived parameters are summarized in Tables 6 and 7 for 49 Ceti and HD 21997, respectively. The optical depths of the 12 CO, 13 CO, and C 18 O lines and residuals of the fit are also listed. The fitting is successful, as shown by the residuals.
Here, we assume that the major collisional partner with CO is H 2 molecules. However, collision with CO, C, C + , and electron(e) would also contribute to the excitation and de-excitation in debris disks. This is particularly important for the low H 2 density case. Since the collisional rate coefficients for the collision with CO/C/C + /O are not available, we cannot distinguish the collision with H 2 and that of CO/C/C + /O in our calculation and just employ the collisional rate coefficients for the H 2 collision. The rate for the CO-H collision is much lower than that of the CO-H 2 collision (Yang et al. 2013;Walker et al. 2015), and thus we have ignored the CO-H collision in this study. For the contribution of CO, the above assumption that the collisional rate coefficients of H 2 are similar to those of CO could be justified by the collisional broadening experiment of the CO line in the laboratory: broadening of the CO line by H 2 and self-broadening are similar to each other, indicating that the collisional rate is not very much different between H 2 and CO (Dick et al. 2009). In short, we have to recognize the derived H 2 density as effective one involving the contribution of CO, C, C + and O: namely, it should approximately be regarded as n(H 2 )+n(CO)+n(C)+n(C + )+n(O). Hereafter, we denote it as the gas density. Note that if the abundance of electron is less than 1% of the above sum, its effect can practically be ignored (Section 3.3.2).
For 49 Ceti, the column density of CO ranges from 5.9×10 17 to 2.0×10 17 cm −2 for the blue-shifted component and from 5.4×10 17 to 1.8×10 17 cm −2 for the red-shifted component ( Table 6). The derived CO column density is higher by 4 orders of magnitude than the upper limit derived from the CO absorption ). The column density and optical depth of CO are higher for the lower gas density. The optical depth is as high as 110 -130 for the gas density of 3×10 3 cm −3 . On the other hand, the gas kinetic temperature is between 8 and 11 K and is not very sensitive to the assumed the gas density. It is similar to the rotational temperature obtained in the LTE analysis. This means that the level population is almost thermalized in the assumed range of the gas density. The beam filling factor is almost independent of the gas density and is consistent with that expected from the source size and the beam size described in Section 3.2. Although the maximum correlation in the fit is 0.99 between T kin and f d , these parameters are well constrained.
We derive the number density of CO from the column density by assuming the line-of-sight length of the emitting region. Since 49 Ceti has an almost edge-on configuration, the disk radius of 100 au is employed as the line-of-sight length. Then, the averaged CO number density is roughly estimated to be 400 cm −3 and 360 cm −3 for the blue-shifted and red-shifted components, respectively, for the gas density of 3×10 3 cm −3 (Table 6). Since the number densities of C and C + are expected to be higher than the CO density by about an order of magnitude (Higuchi et al. 2017), the gas density can be interpreted by contributions of CO, C, and C + even without H 2 molecules in this case. Note that the lower limit of the [C ii] mass is reported to be 2.15×10 −4 M ⊕ from the Herschel observation . This lower limit is lower than the CO mass derived above. Since the spatial and velocity resolutions of the Herschel observation are much different from our observation, it seems difficult to evaluate the contribution of C + . If the abundance of C + are low, CO and C would contribute to the excitation.
The number density of CO and H 2 /CO ratio is derived in this way for the other H 2 density cases, as summarized in Table 6. If the gas density is lower than 10 6 cm −3 , the H 2 /CO ratio is lower than the canonical interstellar value of 10 4 . Since the gas density effectively involves the density of CO, C, and C + , as described above, the actual H 2 /CO ratio would be even lower than those in Table 6. It is thus found that the H 2 /CO ratio can be lower than the canonical value (10 4 ) for interstellar clouds, and the observed line intensities can be reproduced even without H 2 .
A similar analysis is conducted for HD 21997. Since this source is nearly face-on, the line-of-sight length is approximated by the disk scale height H at 300 au, assuming (H/r ∼ 0.06): here we use 20 au in derivation of the averaged number density of CO. The results are summarized in Table 7. The optical depths of the CO(2-1) and CO(3-2) lines are much higher in HD 21997 than those found in 49 Ceti. For this reason, the fitting is not as good as for 49 Ceti, particularly for the blue-shifted component, where the gas kinetic temperature and the CO column density are not well constrained for the gas density below 10 4 cm −3 . Nevertheless, the gas kinetic temperature is below 12 K, which is similar to the 49 Ceti case. Again, the H 2 /CO ratio can be lower than the canonical value for interstellar clouds.
We here note two caveats for the above analysis. First, the optical depth of the CO lines are very high, and hence, the simple approximation by the LVG model (Goldreich & Kwan 1974) employed here does not perfectly describe the radiation transfer. In the LVG model, a photon emitted from a certain volume is assumed not to be absorbed in a different part. This situation may not always be fulfilled for very high optical-depth cases. This might be the reason for the poor fitting in the HD 21997 case. More rigorous treatments including dust opacity are left for future study. For this purpose, we need to resolve the disk structure, and hence, high resolution data of CO and its isotopic species are necessary. Second, the state-to-state collisional rate coefficients for the CO, C, and C + are not available, and therefore the contribution of their collisions are not explicitly considered. Theoretical and experimental studies on the collision rates are awaited.

The case of electron collision
Next, we conduct the non-LTE LVG analysis assuming the collision with an electron, considering that the gas has the secondary origin. Here, we employ the collisional rate coefficients calculated by using the method described by Dickinson et al. (1977). As in the H 2 collision case, we cannot determine the four parameters simultaneously in the analysis, and hence, we fix the electron density and determine the remaining three parameters. Because electrons mostly come from photoionization of C (and organic species), the maximum electron density would be around 10 4 cm −3 , considering the CO number density derived later and the assumption that the C and C + abundances are higher than the CO abundance by an order of magnitude as the case of 49 Ceti Higuchi et al. 2017). Since the collision with an electron is more efficient for rotational excitation/de-excitation by two orders of magnitude than that with H 2 (and possibly CO, C, C + , and O) (Dickinson et al. 1977), the electron abundance relative to CO, C, C + , and O needs to be higher than 0.01 in order for an electron to become a major collision partner. Hence, we set the minimum electron density of 30 cm −3 . The analysis is conducted for several electron densities from 30 to 10 4 cm −3 for 49 Ceti and HD 21997, as shown in Tables 8 and 9, respectively. For 49 Ceti, the column density of CO ranges from 2×10 17 to 1.2×10 18 cm −2 and from 1.8×10 17 to 9.8×10 17 cm −2 for the blue-shifted and red-shifted components, respectively. The derived gas kinetic temperature is around 10 K. These values resemble those found in the H 2 collision case, except for the electron density of 30 cm −3 where the optical depth of the CO lines is very high due to insufficient excitation. The beam filling factors also resemble those in the H 2 collision case. These similarities seem reasonable because the rotational population is almost in LTE regardless of the collision partner. For HD 21997, the ranges of the column density and the gas kinetic temperature are also similar to those of the H 2 collision case. These results indicate that the CO excitation is also possible without H 2 , if the electron density is higher than 30 cm −3 .

Nature of gas components in debris disks
Based on the LTE and non-LTE analyses, we confirm that a huge amount of CO gas is associated with the debris disks of 49 Ceti and HD 21997, and its gas temperature is quite low (∼10 K). This result is essentially consistent with the report based on the LTE analysis for HD21997 by Kóspál et al.
Physical conditions of gas components in debris disks of 49 Ceti and HD 21997 9 (2013), but is based on the more detailed analysis including non-LTE one. The CO mass 2 evaluated from the CO column density, which is corrected the optical depth effect is (6.1 -35)×10 −2 M ⊕ and (5.5 -85)×10 −2 M ⊕ for 49 Ceti and HD 21997, respectively. Our results at the H 2 density of 1.0×10 6 to 1.0×10 7 cm −3 and electron density of 1.0×10 3 to 1.0×10 4 cm −3 are roughly consistent with the mass estimate by Kóspál et al. (2013) and Moór et al. (2019).
The spatial resolution of the present analysis is 6 ′′ for 49 Ceti and 4 ′′ for HD 21997. Since high angular resolution data of the CO isotopologues are not available, we cannot directly derive the radial distribution of the temperature in the disk. We here discuss the effective temperatures derived above in terms of the temperature distribution of equation 1 reported so far.
First, we assume the heating by the central star. For 49 Ceti, T 100 of equation 1 is reported to be 40 K for the H 2 /CO ratio of 10 4 (Hughes et al. 2017). If the CO molecules distributed in the outer disk make a dominant contribution to the observed emission, the lower gas temperature would be expected. However, the gas temperature at R=300 au is estimated to be 23 K, according to Hughes et al. (2017), and is still higher than our result. On the other hand, T 100 is reported to be 14 K if the H 2 /CO ratio is set to be 1. If we employ this value, the gas temperature at R=300 au is 8 K. This estimate is close to the gas kinetic temperature derived in this study. Thus, the low H 2 /CO ratio and the low gas temperature inferred in our analyses (Section 3.3) might be related with each other. Nevertheless, the temperature structure may not be as simple as equation 1. We need more detailed analysis based on high spatial resolution observations. We then consider the effect of the heating by the stellar and interstellar UV radiation on the gas temperature. Both UV radiation fields play an important role in photodissociation of CO in the outer disk, as shown in 49 Ceti (Higuchi et al. 2019a). The [C i]/CO intensity ratio decreases with increasing radial distance, has a minimum at 100 au from the central star, and increases again. However, debris disks contain fewer small grains responsible for photoelectric heating than protoplanetary disks (Wilner et al. 2012;MacGregor et al. 2017). For this reason, the gas temperature would not be as high as those in the surface area of typical protoplanetary disks. Although Besla & Wu (2007) and Zagorovsky et al. (2010) argued that photoelectric heating plays an important role in heating the gas in a β Pictoris-like debris disk, Kral et al. (2016) suggest that it does not contribute to the gas heating. Indeed, the gas temperature in β Pictoris is estimated to be as low as 20 K (Roberge et al. 2000). This is lower than the typical temperature of a diffuse gas in the interstellar clouds illuminated by the interstellar UV radiation (∼ 100 K) (e.g., Spaans 1996;Snow & McCall 2006). Thus the low temperature of the gas components seems possible in debris disks even under the exposure to the interstellar UV radiation.
It is well known that CO molecules are depleted onto dust grains at the dust temperature below its desorption temperature (∼ 20 K) (Caselli et al. 1999;Yamamoto 2017). CO depletion means that CO in the gas phase is frozen out onto dust grains. If the dust temperature is lower than the desorption temperature of CO (∼ 20 K), the CO molecules adsorbed onto dust grain do not come out. If the dust temperature were low at the outer edge of the disk, CO could be depleted onto dust grains. The depletion time scale is approximately 0.1 Myr for the gas density of 10 4 cm −3 if the dust-to-gas mass ratio and the dust size distribution are the same as those in the interstellar clouds (e.g., Aikawa et al. 2012). However, even in an outer region, the small dust grains are less abundant in debris disks (Pawellek et al. 2019) than in the interstellar medium, and hence, the CO desorption timescale is extended (see equation 1 of Aikawa et al. (2012)). Hence, the effect of the CO depletion would not be significant.
For the origin of gas, there are mainly two possibilities. If the gas origin is a primordial case, the H 2 /CO ratio should have been close to the canonical be 10 4 for interstellar clouds or even higher (10 5−6 ) as revealed in TW Hya (Bergin et al. 2013;Bergin & Williams 2018). As for the origin of the low H 2 /CO ratio, we first consider the photoevaporation of the primordial gas from a protoplanetary disk. Photoevaporation (e.g., Nakatani et al. 2018) is in general a hydrodynamic escape, where both CO and H 2 flow out. According to the theoretical models, the dispersal timescale was thought to be 3-5 Myr (e.g., Clarke et al. 2001;Gorti et al. 2009). For forming low H 2 /CO ratio gaseous debris disks, it is widely considered that Jeans escape (e.g., Volkov et al. 2011) may play a role. Since the lighter molecules will readily exceed the escape speed, they will preferentially disappear. Thus, the H 2 molecules may selectively escape from the debris disk. This mechanism works better at the lower gas temperature, as found in this study. However, it is unclear whether it is really possible or not in the physical conditions of the debris disks within their ages.
Another possibility is that the gas has a secondary origin, for which several different mechanisms can be considered; e.g., the sublimation of dust grains or planetesimals, photo-sputtering of dust grains, collisional vaporization of dust grains, and collision of comets or icy planetesimals (e.g., Beust et al. 1990;Grigorieva et al. 2007;Czechowski & Mann 2007;Zuckerman & Song 2012;Matrà et al. 2017;Marino et al. 2020;Kral et al. 2020). Hughes et al. (2017) support this mechanism based on the small scale height of the disk, which is caused by the larger molecular weight due to absence of H 2 . A recent model explains the observed amount of gas by invoking UV shielding by the C atoms .
In order to examine these two possibilities and their relative importance, systematic observations are necessary. In this regard, we only have a limited amount of observational information. For example, available archival data for 49 Ceti are of 12 CO, 13 CO, and C 18 O in Band 6 with ACA observations, 12 CO in Band 7, and [C i] in Band 8 with ALMA observations. For HD 21997, there are 12 CO, 13 CO, and C 18 O in Band 6 and Band 7 with ALMA observations. In addition to them, it is necessary to perform the high-resolution observations of 13 CO(3-2) in Band 7, 12 CO(4-3) and 13 CO(4-3) in Band 8, and 12 CO(7-6) and 13 CO(7-6) in Band 10 to derive the distributions of the gas temperature and the CO column densities. The spatial variation of the CO abundance and its comparison with that of C will be useful to test various models (e.g., Kral et al. 2019) for the origin of the gas in debris disks.
Furthermore, there is a possibility that the gas would contain other molecules (CH 3 OH, SO, SiO, H 2 O, etc.). These molecules are known as important constituents of icy mantle of dust grains (e.g., Nomura et al. 2009;Walsh et al. 2016), and are readily destroyed by UV photodissociation. If they can be detected in debris disks, the contribution of the secondary gas will be evident. Although the outgassing mechanism is not well understood, it is thought to be similar to the process occurring in comets (e.g., 67P/Churyumov-Gerasimenko). Thus, volatiles could be ejected from dust grains to the gas component in debris disks. If the CH 3 OH abundance is derived, we can discuss it in relation to the observational results reported for various comets (Mumma & Charnley 2011). Since a debris disk corresponds to the last phase of planet formation, the result will open a new avenue to study how organic molecules in space are incorporated into baby planets.

SUMMARY
We have analyzed the ALMA archival data of 12 CO and its isotopologues observed toward the gaseous debris disks, 49 Ceti and HD 21997. The major results are summarized as follows: 1. We have conducted the LTE and non-LTE analyses to test the excitation of the CO gas by analyzing the 12 CO, 13 CO, and C 18 O data for 49 Ceti and HD 21997. For the non-LTE analyses, the two distinct cases of the H 2 collision and the electron collision are considered. We have examined a wide range of the gas density, and have determined the beam-averaged column density of CO and the gas temperature.
2. The gas kinetic temperature is derived from the non-LTE analyses to be 8-11 K for 49 Ceti for the first time, which is significantly lower than the dust temperature. A similar temperature (8-12 K) is also obtained for HD 21997. These results would indicate inefficient photoelectric heating due to less abundant small grains.
4. The CO number density averaged over the observation beam is estimated from the CO column density and the scale of the line of sight. It is (100-400) cm −3 and (860-5000) cm −3 for 49 Ceti and HD 21997, respectively.
5. CO molecules can be excited under the environments of no H 2 or a small number of H 2 molecules or even without H 2 . The wide range of H 2 or electron density can account for the CO excitation.
We thank the referee for the thoughtful and constructive comments. This paper makes use of the following ALMA data:    Table 2. Line parameters of 12 CO(3-2), 12 CO(2-1), 13 CO(2-1), C 18 O(2-1) from double Gaussian fits for 49 Ceti Note-The numbers in parentheses represent the 1σ error. T B is the brightness temperature, dv is the velocity dispersion. T CAL represents the calculated intensity by the analysis. For the LVG analysis, we adopted the results at the H 2 density of 1.0×10 7 cm −3 . Table 3. Line parameters of 12 CO(3-2), 12 CO(2-1), 13 CO(3-2), 13 CO(2-1) and C 18 O(2-1) from double Gaussian fits for HD 21997 Note-The numbers in parentheses represent the 1σ error. T B is the brightness temperature, dv is the velocity dispersion.
T CAL represents the calculated intensity by the analysis. For the LVG analysis, we adopted the results at the H 2 density of 1.0×10 7 cm −3 .  Note-The numbers in parentheses represent the 1σ error. N CO is the CO column density, Trot is the rotational temperature, f d is the filling factor. The maximum correlation coefficient is 0.99 between Trot and f d .  a This mimics the primordial gas case.
Note-The numbers in parentheses represent the 1σ error. N CO is the CO column density, T kin is the kinetic temperature, f d is the filling factor. The maximum correlation coefficient is 0.99 between T kin and f d . b This mimics the primordial gas case.
Note-The numbers in parentheses represent the 1σ error. N CO is the CO column density, T kin is the kinetic temperature, f d is the filling factor. The maximum correlation coefficient is 0.99 between T kin and f d . c This mimics the secondary gas case.
Note-The numbers in parentheses represent the 1σ error. N CO is the CO column density, T kin is the kinetic temperature, f d is the filling factor. The maximum correlation coefficient is 0.99 between T kin and f d . Note-The numbers in parentheses represent the 1σ error. N CO is the CO column density, T kin is the kinetic temperature, f d is the filling factor. The maximum correlation coefficient is 0.99 between T kin and f d .