A Constraint on the Amount of Hydrogen from the CO Chemistry in Debris Disks

One of the differences of the debris disk from the ISM is absence of small dust grains. The small dust grains with radii smaller than $\sim 1~{}\mu$m are blown out by the radiation pressure of stellar photons (Burns et al., 1979; Kobayashi et al., 2008, 2009). Absence of the small dust grains reduces the visual extinction, the H₂ formation rate, and the photo-electric heating rate for a given dust-to-gas mass ratio. For simplicity, we consider spherical dust grains whose radii are denoted by $a$. The power-law grain size distribution $n(a)\propto a^{-3.5}$ is assumed, where $n(a)da$ is the number density of the dust grains with the radius range from $a$ to $a+da$ (Mathis et al., 1977). This power-law size distribution is controlled by the collisional cascade (Dohnanyi, 1969; Tanaka et al., 1996), although the power-law index could be modulated due to the size dependence of collisional strength (Kobayashi & Tanaka, 2010). In this work, the minimum and maximum radii of the dust grains are fixed to $a_{\mathrm{min}}=1~{}\mu$m and $a_{\mathrm{max}}=10~{}\mu$m. We will discuss that our results do not depend directly on the size distribution in Section 3.1.

The dust grains are assumed to be a mixture of graphite and silicon. The absorption and scattering coefficients of dust grains are calculated by averaging those of the graphite and silicon in a ratio of $7:3$, taking into account the size distribution. The heat capacity of the dust grains, which is used to determine the dust temperature at each dust size and position, is computed assuming the mixed composition. The absorption/scattering coefficients and heat capacity at each dust size are computed by Laor & Draine (1993) and Draine & Li (2001).

In addition to the FUV flux, the charge of dust grains characterizes the photo-electric yield on the dust grain. The Meudon code takes into account the probability distribution function of the dust charge for each dust size $f(a,Q)$, which is determined by collisional charging and photo-electric ionization (Bakes & Tielens, 1994), where $Q$ is the dust charge.

As the grain size increases, the photo-electric ionization rate increases faster than the recombination rate. As a result, the charge of dust grains increases with $a$. Since the maximum grain size considered in this work is up to 1 cm, a great number of the charge grid is required.

To reduce the computational cost, we approximate the charge distribution function at each grain size $f(a,Q)$ by a delta function of $\delta(Q-\langle Q\rangle_{a})$, where $Q$ is the dust charge and $\langle Q\rangle_{a}$ is the mean charge for the grain size $a$. The validity of the approximation $f(a,Q)=\delta(Q-\langle Q\rangle_{a})$ is confirmed.

Since the dust-to-gas mass ratio $f_{\mathrm{dg}}$ of debris disks is unknown observationally and theoretically, a wide range of $f_{\mathrm{dg}}$ is considered in the present work. Using the size distribution $n(a)$, the dust-to-gas mass ratio is expressed as

$$f_{\mathrm{dg}}\equiv\frac{1}{\mu_{\mathrm{C}}m_{\mathrm{C}}n_{\mathrm{C}}}\frac{4\pi}{3}\rho_{\mathrm{gr}}\langle a^{3}\rangle n_{\mathrm{d}},$$

(3)

where

$$\langle g(a)\rangle\equiv\frac{1}{n_{\mathrm{d}}}\int_{a_{\mathrm{min}}}^{a_{\mathrm{max}}}g(a){n}(a)da,$$

(4)

$n_{\mathrm{d}}=\int_{a_{\mathrm{min}}}^{a_{\mathrm{max}}}{n}(a)da$ is the total number density of dust grains, $m_{\mathrm{x}}$ and $n_{\mathrm{x}}$ are the mass and number density of element “x”, respectively, $\mu_{\mathrm{x}}$ is the mean molecular weight per element “x” nucleus, and $\rho_{\mathrm{gr}}=2.62~{}\mathrm{g~{}\mathrm{cm}^{-3}}$ is the internal density of the dust grains. We should note that the gas mass density is expressed as $\mu_{\mathrm{C}}m_{\mathrm{C}}n_{\mathrm{C}}$ because the gas component is considered on the basis of carbon. The gas metallicity is denoted by $Z$, and the elemental abundances of the gas phase used in the present work are shown in Table 1 for the solar metallicity $Z=1$. When $Z$ changes, the relative abundances among metals are fixed. From Table 1, $\mu_{\mathrm{C}}$ is given by

$$\mu_{\mathrm{C}}(Z)=883\left(Z^{-1}+6.6\times 10^{-3}\right).$$

(5)

In the context of PPDs, $n_{\mathrm{d}}$ is often determined from Equation (3) at a given $f_{\mathrm{dg}}$. By contrast, we use the fact that in debris disks, the amount of dust grains is constrained by the optical depth at V band, which is often estimated from the fraction of the disk luminosity to stellar luminosity. The most practical criterion for debris disks is $\tau<8\times 10^{-3}$ (Hughes et al., 2018), where $\tau$ denotes a typical vertical optical depth of debris disks. A typical optical depth along the mid-plane $\tau_{\mathrm{mid}}$ is given by $\tau\theta^{-1}$, where $\theta$ is the aspect ratio of debris disks.

The mean free path $\ell$ of a photon for dust absorption depends on the grain size distribution as follows:

$$\ell=\left(n_{\mathrm{d}}\langle Q_{\mathrm{abs}}\pi a^{2}\rangle\right)^{-1}.$$

(6)

In the integration, the absorption coefficient $Q_{\mathrm{abs}}$ can be set to be unity because $2\pi a/\lambda>1$ is satisfied for $a\geq 1~{}\mu$m at the V band. Eliminating $n_{\mathrm{d}}$ using Equations (3) and (6), one obtains

$$f_{\mathrm{dg}}=1.74\times 10^{23}~{}\mathrm{cm}^{-3}~{}\frac{\langle a^{3}\rangle}{\langle a^{2}\rangle}\mu_{\mathrm{C}}^{-1}{\cal N}_{\mathrm{C,\tau_{\mathrm{mid}}=1}}^{-1},$$

(7)

where ${\cal N}_{\mathrm{C,\tau_{\mathrm{mid}}=1}}=n_{\mathrm{C}}\ell$ is the mid-plane carbon nucleus column density of the gas phase at $\tau_{\mathrm{mid}}=1$ of the dust opacity at V band. The values of $f_{\mathrm{dg}}$ adopted in our PDR calculations will be shown in Section 2.5.

As mentioned in Section 2.4.1, the size distribution $n(a)\propto a^{-3.5}$ is used, and $\langle a^{3}\rangle/\langle a^{2}\rangle=\sqrt{a_{\mathrm{min}}a_{\mathrm{max}}}$.

In the Meudon PDR code, the temperature of dust grains $T_{\mathrm{gr}}$ is computed from the energy balance between absorption of the radiation and thermal emission at each grain radius bin. For reference, we here show the grain temperature given by Hughes et al. (2008) as follows:

$$T_{\mathrm{gr}}=10^{2}~{}\mathrm{K}\left(\frac{L_{*}}{10L_{\odot}}\right)^{0.2}\left(\frac{a}{1~{}\mu\mathrm{m}}\right)^{-0.2}\left(\frac{r}{50~{}\mathrm{au}}\right)^{-0.4},$$

(8)

where $L_{*}$ is the stellar luminosity, $a$ is the grain radius, and $r$ is the distance from the central star. The dust temperatures obtained the PDR calculations are consistent with Equation (8).

The Meudon code takes into account H₂ and HD formation with the Langmuir-Hinshelwood (LH) and the Eley-Rideal (ER) mechanisms (Le Bourlot et al., 2012). The grain-surface chemistry for other species heavier than H₂ and HD is not included in this study because photo-desorption is so efficient that heavy species are difficult to freeze out on the grain surface (e.g., Grigorieva et al., 2007).

In debris disks, the most important H₂ formation mechanism on dust grains is the ER mechanism involved by chemisorption of H atoms. This is because the dust temperatures, which are as high as $\sim$100 K as shown in Equation (8), are so high that a hydrogen atom physisorbed on the grain surface evaporates before it combines with another hydrogen atom.

In the ER mechanism, there are several free parameters in the Meudon PDR code. One is a sticking coefficient $\alpha(T)$, which approaches zero for high temperatures because the hydrogen atom cannot stick on the dust grains owing to the excess kinetic energy. The Meudon PDR code adopts the empirical form $\alpha(T)=1/(1+(T/T_{\mathrm{stick}})^{\beta})$, where $T_{\mathrm{stick}}=464~{}$K and $\beta=3/2$ are adopted (Le Bourlot et al., 2012).

Another parameter is $T_{\mathrm{chem}}$ corresponding to the energy barrier that a hydrogen atom overcomes to reach a chemisorption site on the grain surface. The H₂ formation rate is highly sensitive to $T_{\mathrm{chem}}$ because the probability of overcoming the energy barrier is proportional to $\exp(-T_{\mathrm{chem}}/T)$. Nonetheless, $T_{\mathrm{chem}}$ is highly uncertain. Its value depends on the surface condition and composition of dust grains. The energy barrier of chemisorption onto a perfect graphite surface is as high as $\sim 2000$ K (e.g., Sha et al., 2002, by using quantum mechanics). This has been confirmed by the experimental study of a highly ordered pyrolytic graphite in the laboratory done by Zecho et al. (2002). By contrast, topological defects on the carbonaceous surface decrease $T_{\mathrm{chem}}$, and can make chemisorption almost barrierless ($T_{\mathrm{chem}}\sim 10-100~{}$K), depending on structures (Ivanovskaya et al., 2010). Experiments of chemisorption of H on porous, defective, aliphatic carbon surfaces have found that the activation energy is as low as $\sim 70$ K (Mennella, 2006). By contrast, in the cases with silicate, quantum chemical calculations revealed $T_{\mathrm{chem}}\sim 300~{}$K both for crystalline silicate (Navarro-Ruiz et al., 2014) and amorphous silicate (Navarro-Ruiz et al., 2015), but they have not been confirmed experimentally yet.

Uncertainty of the H₂ formation rate on warm dust grains also has been discussed in PDRs illuminated by strong radiation from massive stars. Except for the absence of small dust grains, the environments are similar to debris disks. In such PDRs, Habart et al. (2011) observationally found that rotationally-excited H₂ is more abundant than the predictions of PDR models, suggesting that the H₂ production rate on warm dust grains needs to be higher than a typical value (e.g, Jura, 1974). This result suggests that there is significant uncertainty in H₂ formation rates on warm dust grains.

As will be mentioned in Section 2.5, considering the uncertainty of $T_{\mathrm{chem}}$, $T_{\mathrm{chem}}$ is changed as one of the model parameters in our PDR calculations.

As mentioned in Section 2.3, in order to investigate how the CO chemistry depends on radiation fields, the weak-FUV and strong-FUV models are considered. They are characterized by $\chi_{\mathrm{co}}$ and $\chi_{\mathrm{oh}}$.

The amount of dust grains is characterized by ${\cal N}_{\mathrm{C,\tau_{\mathrm{mid}}=1}}^{\mathrm{fid}}$. We here estimate ${\cal N}_{\mathrm{C,\tau_{\mathrm{mid}}=1}}^{\mathrm{fid}}$ from observation results of $\beta$ Pictoris and 49 Ceti, which are examples of the stars having gas-poor and gas-rich debris disks, respectively. Since both the debris disks are almost edge-on, the observed column densities are comparable to those along the mid-plane. The total column density of species “A” integrated along the mid-plane is denoted by ${\cal N}(A)$.

For $\beta$ Pictoris, several authors derive the C⁰ column densities from the [CI] ${}^{3}P_{1}$-${}^{3}P_{0}$ emission line (Higuchi et al., 2017; Cataldi et al., 2018) and ${}^{3}P$ absorption lines (Roberge et al., 2000). Their values range from ${\cal N}(\mathrm{C}^{0})\sim 2\times 10^{16}$ to $7\times 10^{16}~{}\mathrm{cm}^{-2}$. The C⁺ column densities are estimated to be $\sim 2\times 10^{16}-12\times 10^{16}~{}\mathrm{cm}^{-2}$ from the [CII] 158 $\mu$m emission line (Cataldi et al., 2014, 2018) and absorption lines (Roberge et al., 2006). The CO column densities are estimated to be $3\times 10^{14}-5\times 10^{14}~{}\mathrm{cm}^{-2}$ from emission lines (Higuchi et al., 2017) and $\sim 6\pm 0.3\times 10^{14}~{}\mathrm{cm}^{-2}$ from absorption lines (Roberge et al., 2000). These observations suggest that ${\cal N}(\mathrm{C}^{0})$ is comparable to ${\cal N}(\mathrm{C}^{+})$, and ${\cal N}(\mathrm{CO})$ is roughly two orders of magnitude smaller than ${\cal N}(\mathrm{C}^{0})$ and ${\cal N}(\mathrm{C}^{+})$. The carbon nucleus column density ${\cal N}_{\mathrm{C}}$ is estimated to be $\sim{\cal N}(\mathrm{C}^{0})+{\cal N}(\mathrm{C}^{+})\sim 10^{17}~{}\mathrm{cm}^{-2}$. The spatial distribution of the dust surface area obtained by HST observations of Heap et al. (2000) is fitted by Fernández et al. (2006). The mid-plane optical depth is about $\tau_{\mathrm{mid}}\sim 10^{-2}$, where the fitting formula is integrated over $50~{}\mathrm{au}\leq R\leq 120~{}\mathrm{au}$, which is the gas extent of the best fit model in Cataldi et al. (2018). Finally, ${\cal N}_{\mathrm{C,\tau_{\mathrm{mid}}=1}}^{\mathrm{\beta Pic}}\sim{\cal N}_{\mathrm{C}}/\tau_{\mathrm{mid}}\sim 10^{19}~{}\mathrm{cm}^{-2}$ is obtained for $\beta$ Pictoris.

For 49 Ceti, the C⁰ and CO column densities are estimated to be ${\cal N}(\mathrm{C}^{0})\sim 2\times 10^{18}~{}\mathrm{cm}^{-2}$ (Higuchi et al., 2019) and ${\cal N}(\mathrm{CO})\sim 1.8\times 10^{17}-5.9\times 10^{17}~{}\mathrm{cm}^{-2}$ (Higuchi et al., 2020), respectively. The C⁺ column density is not well constrained because only a lower limit of the C⁺ mass, which is smaller than the CO mass, was obtained (Roberge et al., 2013). We however expect that ${\cal N}_{\mathrm{C}}$ is comparable to ${\cal N}(\mathrm{C}^{0})$ because ${\cal N}(\mathrm{C}^{+})$ should not be much larger than ${\cal N}(\mathrm{C}^{0})$ if all CO are formed through chemical reactions. The vertical optical depth estimated from the fractional luminosity is $\tau\sim 10^{-3}$ (Jura et al., 1993, 1998). Assuming that the aspect ratio of the dust distribution is $\sim 0.1$, the mid-plane optical depth is $\tau_{\mathrm{mid}}\sim 10^{-2}$. Finally, we obtain ${\cal N}_{\mathrm{C,\tau_{\mathrm{mid}}=1}}^{\mathrm{49Ceti}}\sim 2\times 10^{20}~{}\mathrm{cm}^{-2}$ for 49 Ceti.

From ${\cal N}_{\mathrm{C,\tau_{\mathrm{mid}}=1}}^{\mathrm{\beta Pic}}=10^{19}~{}\mathrm{cm}^{-2}$ and ${\cal N}_{\mathrm{C,\tau_{\mathrm{mid}}=1}}^{\mathrm{49Ceti}}=2\times 10^{20}~{}\mathrm{cm}^{-2}$, a fiducial value of ${\cal N}_{\mathrm{C,\tau_{\mathrm{mid}}=1}}={\cal N}_{\mathrm{C}}/\tau_{\mathrm{mid}}$ is set to ${\cal N}_{\mathrm{C,\tau_{\mathrm{mid}}=1}}^{\mathrm{fid}}=2\times 10^{19}~{}\mathrm{cm}^{-2}$ (${\cal N}_{\mathrm{C,\tau_{\mathrm{mid}}=1}}^{\mathrm{\beta Pic}}=0.5{\cal N}_{\mathrm{C,\tau_{\mathrm{mid}}=1}}^{\mathrm{fid}}$, ${\cal N}_{\mathrm{C,\tau_{\mathrm{mid}}=1}}^{\mathrm{49Ceti}}=10{\cal N}_{\mathrm{C,\tau_{\mathrm{mid}}=1}}^{\mathrm{fid}}$) because the behavior of the CO formation changes around ${\cal N}_{\mathrm{C,\tau_{\mathrm{mid}}=1}}\sim{\cal N}_{\mathrm{C,\tau_{\mathrm{mid}}=1}}^{\mathrm{fid}}$ for weak-FUV and strong-FUV as will be shown in Section 3.

In this paper, instead of ${\cal N}_{\mathrm{C,\tau_{\mathrm{mid}}=1}}$, a normalized dust-to-gas mass ratio is used as the parameter indicating the amount of dust grains. Using Equation (7), one obtains

$$\frac{f_{\mathrm{dg}}}{f_{\mathrm{dg,fid}}}=\left(\frac{{\cal N}_{\mathrm{C,\tau_{\mathrm{mid}}=1}}}{{\cal N}_{\mathrm{C,\tau_{\mathrm{mid}}=1}}^{\mathrm{fid}}}\right)^{-1},$$

(9)

where $f_{\mathrm{dg,fid}}$ is the dust-to-gas mass ratio at ${\cal N}_{\mathrm{C,\tau_{\mathrm{mid}}=1}}={\cal N}_{\mathrm{C,\tau_{\mathrm{mid}}=1}}^{\mathrm{fid}}$ and given by

	$\displaystyle f_{\mathrm{dg,fid}}(a_{\mathrm{min}},a_{\mathrm{max}},Z)$	$\displaystyle=$	$\displaystyle 3\times 10^{-3}\left\{Z^{-1}+6.6\times 10^{-3}\right\}^{-1}$		(10)
			$\displaystyle\times\left(\frac{a_{\mathrm{min}}}{1~{}\mathrm{\mu m}}\right)^{1/2}\left(\frac{a_{\mathrm{max}}}{10~{}\mathrm{\mu m}}\right)^{1/2},$

where $n(a)\propto a^{-3.5}$ is used. We should note that although both $f_{\mathrm{dg}}$ and $f_{\mathrm{dg,fid}}$ depend on $a_{\mathrm{min}}$ and $a_{\mathrm{max}}$, the ratio $f_{\mathrm{dg}}/f_{\mathrm{dg,fid}}$ is independent of the size distribution for a given ${\cal N}_{\mathrm{C,\tau_{\mathrm{mid}}=1}}$.

To study how the amount of dust grains affects the thermal and chemical structures, the PDR calculations are performed with three different $f_{\mathrm{dg}}$ as listed in Table 2. For weak-FUV, since a typical $f_{\mathrm{dg}}$ is about $\sim 2f_{\mathrm{dg,fid}}$ for $\beta$ Pictoris, the cases with ($0.1f_{\mathrm{dg,fid}}$, $f_{\mathrm{dg,fid}}$, $10f_{\mathrm{dg,fid}}$) are considered. For strong-FUV, since a typical $f_{\mathrm{dg}}$ is about $\sim 0.1f_{\mathrm{dg,fid}}$, the cases with ($10^{-2}f_{\mathrm{dg,fid}}$, $0.1f_{\mathrm{dg,fid}}$, $f_{\mathrm{dg,fid}}$) are considered.

The parameters associated with the gas component are $n_{\mathrm{C}}$ and metallicity $Z$. Given $n_{\mathrm{C}}$ and $Z$, the corresponding hydrogen nucleus number density is $n_{\mathrm{H}}=n_{\mathrm{C}}({\cal A}_{\mathrm{C}}Z)^{-1}$, where ${\cal A}_{\mathrm{C}}=1.32\times 10^{-4}$ is the carbon elemental abundance of the gas for $Z=1$ (Table 1). The hydrogen gas density is difficult to constrain in most debris disks observationally although there are some successful exceptions which will be described in Section 4.4.

Fiducial carbon nucleus number densities are defined as $n_{\mathrm{C}}=1.32\times 10^{2}~{}\mathrm{cm}^{-3}$ for weak-FUV and $n_{\mathrm{C}}=1.32\times 10^{3}~{}\mathrm{cm}^{-3}$ for strong-FUV. The debris disk around $\beta$ Pictoris has a spatial extent of $\sim 70~{}\mathrm{au}$ in the best-fit model in Cataldi et al. (2018, also see Section 4.2.2). A mean density is $n_{\mathrm{C}}\sim 10^{17}~{}\mathrm{cm}^{-2}/70~{}\mathrm{au}\sim 100~{}\mathrm{cm}^{-3}$, which is comparable to the fiducial $n_{\mathrm{C}}$ for weak-FUV. For 49 Ceti, because the radial extent is $\sim 100~{}$au, which corresponds to a cut-off radius of the surface density distribution (Hughes et al., 2018, also see Section 4.2.3), the mean density is $n_{\mathrm{C}}\sim 2\times 10^{18}~{}\mathrm{cm}^{-2}/100~{}\mathrm{au}\sim 10^{3}~{}\mathrm{cm}^{-3}$, which is comparable to the fiducial $n_{\mathrm{C}}$ for strong-FUV. For comparison, the densities 10 times larger than the typical $n_{\mathrm{C}}$ are also considered for both the weak-FUV and strong-FUV models.

Next, we consider a range of the gas metallicity. The gas metallicity in debris disks depends on its origin. If the gas is a remnant of PPDs, the gas metallicity is expected not to be significantly different from the solar metallicity $Z=1$. It is maximized when all the gas comes from the secondary processes; hydrogen is provided by photo-dissociation of H₂O outgassed from dust grains. In this case, the number of hydrogen nuclei becomes comparable to that of oxygen nuclei, and the gas metallicity reaches $1.6\times 10^{3}$ since ${\cal A}_{\mathrm{O}}=3.2\times 10^{-4}$ at $Z=1$ (Table 1). In this paper, a range of $1\leq Z\leq 10^{3}$ is considered.

In order to investigate how the uncertainty of the H₂ formation affects CO formation (Section 2.4.6), the cases with $T_{\mathrm{chem}}=300$ K and $10$ K are considered. The fiducial case $T_{\mathrm{chem}}=300~{}$K leads to a moderate H₂ formation rate, corresponding to situations where chemisorption efficiently occurs on grains with plenty of surface defects (Le Bourlot et al., 2012). The cases with $T_{\mathrm{chem}}=10~{}$K give an upper limit of the H₂ formation rate. Other parameters $T_{\mathrm{stick}}$ and $\beta$ are fixed to $464$ K and 1.5, respectively, because the gas temperatures do not exceed $\sim 500$ K in most cases.

Figures 2a and 2b show that the gas temperatures increase with $f_{\mathrm{dg}}$ for both the weak-FUV and strong-FUV models because the photo-electric heating rate from dust grains, which is one of the dominant heating processes, increases in proportion to $f_{\mathrm{dg}}/f_{\mathrm{dg,fid}}$. A weak positive dependence of $T$ on $T_{\mathrm{chem}}$ is attributed to enhancement of the heating rate owing to the H₂ formation on the grain surface.

Comparison between Figures 2a and 2b shows that the gas temperatures increase with $\chi_{\mathrm{co}}$ because the photo-electric heating rate increases with $\chi_{\mathrm{co}}$.

Before presenting the results, we introduce the analytic formation rate of H₂ on grain surfaces,

$$F_{\mathrm{H_{2},dust}}=3\times 10^{-3}R_{\mathrm{H_{2},ISM}}~{}Z\left(\frac{f_{\mathrm{dg}}}{f_{\mathrm{dg,fid}}}\right)\sqrt{T}\kappa(T)n_{\mathrm{H}}n(\mathrm{H})$$

(12)

(Le Bourlot et al., 2012), where Equation (11) is used, $R_{\mathrm{H_{2},ISM}}=3\times 10^{-17}~{}$cm³ s^-1 is a typical rate coefficient in the ISM (Jura, 1974), and $\kappa(T)$ is the chemisorption efficiency⁵⁵5 We should note that the Meudon code does not solve Equation (12) directly, but solves a set of rate equations considering chemisorbed hydrogen atoms and the dust size distribution (see Le Bourlot et al., 2012, for details). In addition to the chemisorbed H atoms, the Meudon code treats H₂ formation with physisorbed H atoms. . Le Bourlot et al. (2012) derived

$$\kappa(T)=\frac{\alpha(T)\exp\left(-T_{\mathrm{chem}}/T\right)}{1+\alpha(T)\exp\left(-T_{\mathrm{chem}}/T\right)},$$

(13)

where $\alpha(T)$ and $T_{\mathrm{chem}}$ are defined in Section 2.4.6. For gas temperatures lower than $T_{\mathrm{stick}}=464~{}$K, which is satisfied in all the cases shown in Figure 2, the H₂ formation rate is proportional to $\exp(-T_{\mathrm{chem}}/T)$, and it depends sensitively on $T$.

First, the cases with $T_{\mathrm{chem}}=300~{}$K are considered. Although H₂ mainly forms on grain surfaces for the ISM, this may not be the case for debris disks because the H₂ formation rate is extremely small. In order to examine the contribution of the grain-surface chemistry to H₂ formation, we performed the additional PDR calculations where the grain-surface chemistry is switched off and plotted the results by the thin lines in the middle column of Figure 2. For $f_{\mathrm{dg}}\leq f_{\mathrm{dg,fid}}$ (weak-FUV) and $f_{\mathrm{dg}}\leq 0.1f_{\mathrm{dg,fid}}$ (strong-FUV), there are almost no changes in the H₂ fractions with and without the grain-surface chemistry (Figures 2c and 2d), indicating that the grain-surface chemistry does not contribute to H₂ formation. In this case, H₂ is mainly formed in the gas phase, and the dominant chemical reaction is $\mathrm{CH^{+}+H\rightarrow H_{2}+C^{+}}$, where CH⁺ is formed by the radiative association between C⁺ and H⁶⁶6 Especially for metal-poor environments, the so-called H^- channel is an important H₂ formation path (e.g., Sternberg et al., 2021). Additional PDR calculations were performed with chemical reactions associated with H^-, which had not been activated in the default setting of the Meudon code. The photo-detachment of H^- is computed directly from $\int\sigma_{\lambda}I_{\lambda}/h\nu d\lambda$ (Section 2.2), where the cross section obtained in Wishart (1979) is used while the Meudon code calculates it by using a fitting formula. We confirmed that the H^- channel does not affect the CO chemistry. .

The H₂ formation rate of the grain-surface chemistry is enhanced in proportion to $f_{\mathrm{dg}}$ (Equation (12)) while the rate in the gas phase does not change. As shown in Figures 2c and 2d, $f_{\mathrm{dg}}$ is high enough for H₂ to be formed mainly through the grain-surface chemistry when $f_{\mathrm{dg}}=10f_{\mathrm{dg,fid}}$ for weak-FUV and $f_{\mathrm{dg}}=f_{\mathrm{dg,fid}}$ for strong-FUV.

For $T_{\mathrm{chem}}=300~{}$K, the H₂ fractions stay low level even after they increase rapidly around $\sim 10^{14}~{}\mathrm{cm}^{-2}$ owing to self-shielding (Draine & Bertoldi, 1996). This is because the H₂ formation rate is so low that the destruction reactions in the gas phase and H₂ photo-dissociation keep the H₂ fractions low.

As mentioned in Section 2.4.6, there is large uncertainty in H₂ formation on the warm dust grains with $T_{\mathrm{gr}}\sim 100~{}$K. A decrease in $T_{\mathrm{chem}}$ from 300 K to 10 K significantly enhances the H₂ fractions in Figures 2c and 2d because $\kappa(T)$ depends sensitively on $T_{\mathrm{chem}}$. Considering the gas temperatures shown in Figures 2a and 2b, the H₂ formation rate increases by a factor of $\sim 1.6\times 10^{4}$ at $T=30$ K and $\sim 1.3\times 10^{2}$ at $T=60$ K.

First, the models with $T_{\mathrm{chem}}=300~{}$K are considered. Figures 2e and 2f show that $n(\mathrm{CO})/n_{\mathrm{C}}$ does not depend on $f_{\mathrm{dg}}$ significantly as long as $f_{\mathrm{dg}}\leq f_{\mathrm{dg,fid}}$ for both weak-FUV and strong-FUV. The CO fractions decrease with increasing $\chi_{\mathrm{co}}$ because stronger FUV flux destroys CO more efficiently.

When $f_{\mathrm{dg}}$ is increased from $f_{\mathrm{dg,fid}}$ to $10f_{\mathrm{dg,fid}}$, the weak-FUV models show that the CO fractions increase by a factor of 2 or 3.

Acceleration of the CO formation owing to H₂ is more significant for $T_{\mathrm{chem}}=10~{}$K than for $T_{\mathrm{chem}}=300$ K since a decrease in $T_{\mathrm{chem}}$ increases the H₂ fraction (Section 3.1.2). The CO fractions monotonically increase with $f_{\mathrm{dg}}$.

The $f_{\mathrm{dg}}$-dependence of the CO fractions disappears when the grain-surface chemistry is omitted from the PDR calculations. This clearly indicates that the enhancement of the CO fraction for larger $f_{\mathrm{dg}}$ is related to efficient H₂ formation on the grain surface.

How much H₂ is needed to accelerate CO formation? Comparing between the middle and right columns of Figure 2, one can see that H₂ fraction must be well above 0.1 for the CO fractions to increase by more than on order of magnitude. When the H₂ fraction is comparable to 0.1 as in the strong-FUV models with ($f_{\mathrm{dg}}=0.01f_{\mathrm{dg,fid}}$, $T_{\mathrm{chem}}=10~{}$K) and ($f_{\mathrm{dg}}=f_{\mathrm{dg,fid}}$, $T_{\mathrm{chem}}=300~{}$K), CO formation is not accelerated significantly.

In the deep interior $N_{\mathrm{C}}>5\times 10^{16}~{}\mathrm{cm}^{-2}$, the CO fractions increase owing to shielding effects, which will be discussed in Section 3.2.3.

Figure 3 shows the spatial distributions of the CO fractions for various $\chi$, $n_{\mathrm{C}}$, and $Z$ listed in Table 2. The dust-to-gas mass ratios are fixed to $f_{\mathrm{dg}}=0.1f_{\mathrm{dg,fid}}$ for weak-FUV and $f_{\mathrm{dg}}=0.01f_{\mathrm{dg,fid}}$ for strong-FUV, where the CO formation proceeds, regardless of H₂ (Figure 2).

We first focus on the CO fractions for $N_{\mathrm{C}}\lesssim 10^{17}~{}\mathrm{cm}^{-2}$. For given $n_{\mathrm{C}}$ and the FUV fluxes ($\chi_{\mathrm{co}},\chi_{\mathrm{oh}}$), the CO fractions decrease with increasing $Z$. Since $Z$ is changed keeping $n_{\mathrm{C}}$ constant, $n_{\mathrm{H}}$ decreases with increasing $Z$. The CO fractions thus decrease with decreasing $n_{\mathrm{H}}$. This behavior has been pointed out by Higuchi et al. (2017); hydrogen is involved in CO formation.

Comparing between the panels of Figure 3, one can see that the CO fractions increase as $n_{\mathrm{C}}$ increases and/or $\chi_{\mathrm{co}}$ decreases. This behavior was qualitatively expected because the CO formation rate will increase with increasing $n_{\mathrm{C}}$ and the CO destruction rate will decrease with decreasing $\chi_{\mathrm{co}}$.

For all the models shown in Figure 3, the CO fractions begin to increase in deep interiors owing to shielding effects, which will be investigated in Section 3.2.3.

In this section, we develop an analytic formula for the CO fractions from the most important chemical reactions associated with CO formation of the H₂-free gas in Figure 4. There are three main paths to form CO. At the high density limit (Equation (23)) , the CO formation is intermediated mainly by OH (Path_OH). For lower densities, the dominant CO formation path depends on $Z$. At $Z\leq 10$, the CO formation is intermediated by CH⁺ (Path${}_{\mathrm{CH^{+}}}$). For higher $Z$, the CO formation proceeds without hydrogen as shown by the green arrows in Figure 4 (Path_woH).

The analytic formula taken into account all the chemical reactions shown in Figure 4 is complicated and is not useful. In Appendix A, we construct a fitting formula of the CO fraction based on Path_OH, which is dominated for higher densities, as follows:

$$\frac{n(\mathrm{CO})_{\mathrm{ana}}}{n_{\mathrm{C}}}=\frac{{\cal A}_{\mathrm{O}}}{{\cal A}_{\mathrm{O,ism}}}\left(10^{-14}\eta^{1.8}+6.0\times 10^{-11}\eta\right),$$

(14)

where

$$\eta=n_{\mathrm{H}}Z^{0.4}\chi^{-1.1},$$

(15)

$\chi\equiv\sqrt{\chi_{\mathrm{co}}\chi_{\mathrm{oh}}}$, ${\cal A}_{\mathrm{O}}$ is the relative oxygen abundance with respect to hydrogen in the gas phase, and ${\cal A}_{\mathrm{O,ism}}=3.2\times 10^{-4}Z$ corresponds to ${\cal A}_{\mathrm{O}}$ shown in Table 1.

We should note that Equation (14) is based on the chemical reactions in regions where shielding effects are not important. However, by considering the shielding effect in $\chi$, we will show Equation (14) is applicable also in regions where shielding effects are effective in Section 3.2.3.

The analytic formula $n(\mathrm{CO})_{\mathrm{ana}}$ at a given position depends on $n_{\mathrm{H}}$, $Z$, and $\chi=\sqrt{\chi_{\mathrm{co}}\chi_{\mathrm{oh}}}$. Although $n_{\mathrm{H}}$ and $Z$ are given locally, $\chi$ is determined by the column densities integrated from sources of light. The FUV radiation that destroys CO can be attenuated by dust grains, C⁰, CO, and H₂. In debris disks, the dust extinction is negligible.

Which of C⁰, CO, and H₂ is more important for the shielding effects? In each shielding effect, one can define the critical column density at which $\chi_{\mathrm{co}}$ is halved. Their critical densities are

	$\displaystyle N(\mathrm{H_{2}})_{\mathrm{shld}}\sim 4\times 10^{20}~{}\mathrm{cm}^{-2},$
	$\displaystyle N(\mathrm{C^{0}})_{\mathrm{shld}}\sim 4\times 10^{16}~{}\mathrm{cm}^{-2}$		(16)
	$\displaystyle N(\mathrm{CO})_{\mathrm{shld}}\sim 10^{14}~{}\mathrm{cm}^{-2}.$

$N(\mathrm{H_{2}})_{\mathrm{shld}}$ is given by the results in van Dishoeck & Black (1988); we derive the H₂ column density at which the shielding factor becomes 0.5 by the linear interpolation using Table 5 of their paper at $\log_{10}N(\mathrm{CO})=0$ (also see Visser et al., 2009). $N(\mathrm{C}^{0})_{\mathrm{shld}}$ is given by $(\ln 2)\alpha_{\mathrm{C}}^{-1}$, where $\alpha_{\mathrm{C}}=1.777\times 10^{-17}~{}\mathrm{cm}^{2}$ is the cross section of photo-ionization of C⁰ (Heays et al., 2017), because the shielding factor owing to C⁰ is proportional to $\exp(-\alpha_{\mathrm{C}}N(\mathrm{C}^{0}))$. $N(\mathrm{CO})_{\mathrm{shld}}$ is taken from the CO column density at which the shielding factor is 0.5 at $\log_{10}N(\mathrm{H}_{2})=0$ in Table 5 of Visser et al. (2009, also see Equation (21)).

Let us estimate the importance of the shielding effect owing to H₂. Considering the abundance ratio between hydrogen and carbon, one finds that the H₂ column density is expressed as $N(\mathrm{H}_{2})=3.8\times 10^{3}Z^{-1}N(\mathrm{C}^{0})$, where it is assumed that all hydrogen is in H₂ ($n(\mathrm{H}_{2})=n_{\mathrm{H}}/2$) and all carbon is in C⁰ ($n(\mathrm{C}^{0})=n_{\mathrm{C}}$). By contrast, in order for the shielding effect owing to H₂ to be more important than the C⁰ attenuation, $N(\mathrm{H}_{2})$ should be larger than $10^{4}N(\mathrm{C}^{0})$, where the numerical factor $10^{4}$ corresponds to $N(\mathrm{H}_{2})_{\mathrm{shld}}/N(\mathrm{C}^{0})_{\mathrm{shld}}$. This condition is not satisfied even for $Z=1$, and the shielding effect owing to H₂ is not important compared with the C⁰ attenuation. For simplicity, the shielding by H₂ is neglected in the analytic formula. We should note that the photo-dissociation rate may be overestimated in significantly shielded regions with $N(\mathrm{H_{2}})>N(\mathrm{H_{2}})_{\mathrm{shld}}$.

In debris disks, C⁰ and CO contribute mainly to a decrease in $\chi_{\mathrm{co}}$. The problem is that the C⁰ and CO fractions both depend on local $\chi_{\mathrm{co}}$, which depends on their column densities. In this section, we construct a procedure to determine the the spatial distributions of $n(\mathrm{C^{0}})$, $n(\mathrm{CO})$, and $\chi_{\mathrm{co}}$ consistently.

Before presenting a way to evaluate the shielding effects, we should mention about the analytic formulae for $n(\mathrm{CO})$ and $n(\mathrm{C^{0}})$. The analytic formula for $n(\mathrm{CO})$ shown in Equation (14) is not applicable directly to a situation where carbon nuclei are mostly in CO because $n(\mathrm{CO})_{\mathrm{ana}}$ increases without limit as $\eta$ increases. To address this issue, Equation (14) is simply modified so that $n(\mathrm{CO})_{\mathrm{ana}}$ approaches $n_{\mathrm{C}}$ smoothly in $\eta\rightarrow\infty$ as follows:

	$\displaystyle\frac{n(\mathrm{CO})_{\mathrm{ana}}}{n_{\mathrm{C}}}=$
	$\displaystyle\left[1+\left\{\frac{{\cal A}_{\mathrm{O}}}{{\cal A}_{\mathrm{O,ism}}}\left(10^{-14}\eta^{1.8}+6.0\times 10^{-11}\eta\right)\right\}^{-1}\right]^{-1}.$		(17)

This modification may be ad hoc, but it predicts the CO fractions consistent with those obtained from the PDR calculations as will be shown in Figure 3.

An analytic formula for the C⁰ fraction is presented. The C⁰ fraction is determined by the balance between photo-ionization of C⁰ and radiative recombination of C⁺ (Kamp & Bertoldi, 2000) as follows:

$$\frac{n(\mathrm{C}^{0})_{\mathrm{ana}}}{n_{\mathrm{C}}-n(\mathrm{CO})_{\mathrm{ana}}}=\frac{2\xi+1-\sqrt{1+4\xi}}{2\xi},$$

(18)

where the denominator $n_{\mathrm{C}}-n(\mathrm{CO})_{\mathrm{ana}}$ guarantees that the sum of the C⁰, C⁺, and CO number densities is equal to $n_{\mathrm{C}}$, and $\xi$ is the ratio of the recombination to the ionization coefficients that is given by

	$\displaystyle\xi$	$\displaystyle\equiv$	$\displaystyle\frac{k_{\mathrm{rec}}(n_{\mathrm{C}}-n(\mathrm{CO})_{\mathrm{ana}})}{\alpha_{\mathrm{C}}\chi_{\mathrm{co}}F_{\mathrm{H}}}$		(19)
		$\displaystyle=$	$\displaystyle\chi_{\mathrm{co}}^{-1}\left(\frac{n_{\mathrm{C}}-n(\mathrm{CO})_{\mathrm{ana}}}{11~{}\mathrm{cm}^{-3}}\right)\left(\frac{T}{100~{}\mathrm{K}}\right)^{-0.82}.$

We confirmed that Equation (18) reproduced the C⁰ fractions obtained from the PDR calculations.

The normalized FUV flux $\chi_{\mathrm{co}}$ is expressed in terms of the unattenuated value $\chi_{\mathrm{co,thin}}$ as follows:

$$\chi_{\mathrm{co}}=\chi_{\mathrm{co,thin}}f_{\mathrm{shield}}(N(\mathrm{C}^{0})_{\mathrm{ana}},N(\mathrm{CO})_{\mathrm{ana}}),$$

(20)

where $f_{\mathrm{shield}}$ corresponds to the shielding factor, which depends on the C⁰ and CO column densities ($N(\mathrm{C^{0}})_{\mathrm{ana}}$ and $N(\mathrm{CO})_{\mathrm{ana}}$) computed from $n(\mathrm{C}^{0})_{\mathrm{ana}}$ and $n(\mathrm{CO})_{\mathrm{ana}}$, respectively. As the shielding factor, we adopt the following expression,

$$f_{\mathrm{shld}}=e^{-\alpha_{\mathrm{C}}N(\mathrm{C}^{0})_{\mathrm{ana}}}\left[1+\left(\frac{N(\mathrm{CO})_{\mathrm{ana}}}{10^{14}~{}\mathrm{cm}^{-2}}\right)^{0.6}\right]^{-1},$$

(21)

where the first factor on the right-hand side corresponds to the C⁰ attenuation, and the second factor is a fitting function of the self-shielding factor at $N(\mathrm{H}_{2})=0$ tabulated in Visser et al. (2009).

The spatial distributions of $\chi_{\mathrm{co}}$, $n(\mathrm{C}^{0})_{\mathrm{ana}}$, and $n(\mathrm{CO})_{\mathrm{ana}}$ are determined consistently in a iterative manner by using Equations (17), (18), and (20).

It is useful to derive $n_{\mathrm{H}}$ required to produce a specific value of the CO fraction. In the high density limit where Path_OH is dominated, since the second term in the right-hand side of Equation (14) is negligible, Equation (14) is solved for $n_{\mathrm{H}}$ as follows:

	$\displaystyle n_{\mathrm{H,req}}$	$\displaystyle=$	$\displaystyle 6\times 10^{7}~{}\mathrm{cm}^{-3}~{}Z^{-0.4}\chi^{1.1}$		(22)
			$\displaystyle\;\;\;\;\times\left(\frac{{\cal A}_{\mathrm{O}}}{{\cal A}_{\mathrm{O,ism}}}\right)^{-0.56}\left(\frac{n(\mathrm{CO)}}{n_{\mathrm{C}}}\right)^{0.56},$

which is valid when

$$n_{\mathrm{H,req}}>5.3\times 10^{4}~{}\mathrm{cm}^{-3}~{}Z^{-0.4}\chi^{1.1}.$$

(23)

Equation (22) is rewritten as

	$\displaystyle n_{\mathrm{C,req}}$	$\displaystyle=$	$\displaystyle 8\times 10^{3}~{}\mathrm{cm}^{-3}Z^{0.6}\chi^{1.1}$		(24)
			$\displaystyle\;\;\;\;\times\left(\frac{{\cal A}_{\mathrm{O}}}{{\cal A}_{\mathrm{O,ism}}}\right)^{-0.56}\left(\frac{n(\mathrm{CO})}{n_{\mathrm{C}}}\right)^{0.56},$

which corresponds to the carbon nucleius density required to reproduce a specific value of $n(\mathrm{CO})/n_{\mathrm{C}}$. From Equation (24), an important conclusion can be drawn that $n(\mathrm{CO})/n_{\mathrm{C}}$ decreases with increasing $Z$ for fixed $n_{\mathrm{C}}$ and $\chi$. In other words, the upper limit of the CO fraction is obtained at $Z=1$.

Figure 3 compares the results of the PDR calculations with the predictions from the analytic formula. $n(\mathrm{CO})_{\mathrm{ana}}$ reproduces the spatial profiles of $n(\mathrm{CO})$ reasonably well from the unattenuated regions ($N_{\mathrm{C}}\ll 10^{17}~{}\mathrm{cm}^{-2}$) to the shielded regions ($N_{\mathrm{C}}\gtrsim 10^{17}~{}\mathrm{cm}^{-2}$) for all the parameter sets.

For lower $n_{\mathrm{C}}$ shown in Figures 3a and 3c, the CO fractions begin to increase around $N_{\mathrm{C}}\sim 10^{16-17}~{}\mathrm{cm}^{-2}$, regardless of $Z$ and $\chi$. By contrast, Figures 3b and 3d show that the CO fractions with $Z=1$ begin to increase at lower $N_{\mathrm{C}}$. This comes from the fact that CO self-shielding only works if the CO fraction is high enough; otherwise, C⁰ attenuation becomes more important.

A critical CO fraction above which CO self-shielding effect works is given by

$$\left(\frac{n(\mathrm{CO})}{n_{\mathrm{C}}}\right)_{\mathrm{cri}}\sim\frac{N(\mathrm{CO})_{\mathrm{shld}}}{N(\mathrm{C}^{0})_{\mathrm{shld}}}\sim 3\times 10^{-3}.$$

(25)

The horizontal dashed lines in Figure 3 correspond to $(n(\mathrm{CO})/n_{\mathrm{C}})_{\mathrm{cri}}$. In order to illustrate the effect of the CO self-shielding in Figure 3, we overplot the profiles of $n(\mathrm{CO})_{\mathrm{ana}}$ without considering the CO self-shielding effect in Equation (21) by the thick dashed lines. It is clearly seen that the CO self-shielding effect increases the CO fractions at $N_{\mathrm{C}}<N(\mathrm{C}^{0})_{\mathrm{shld}}$ only when $n(\mathrm{CO})$ exceeds $n(\mathrm{CO})_{\mathrm{cri}}$.

The importance of C⁰ attenuation in CO formation has been pointed out by Kral et al. (2019).