H$_2$ ortho-para spin conversion on inhomogeneous grain surfaces

TL;DR

This study models how the nuclear spin conversion of H$_2$ on water ice-coated grain surfaces varies with physical conditions, impacting star and planet formation chemistry.

Contribution

It introduces a detailed rate equation model for H$_2$ spin conversion on inhomogeneous grain surfaces, considering adsorption, thermal hopping, and desorption processes.

Findings

Spin conversion efficiency depends on H$_2$ density and surface temperature.

Maximum efficiency occurs at a critical surface temperature.

Lower activation barriers enhance spin conversion efficiency.

Abstract

We investigate the evolution of the ortho-to-para ratio of overall (gas + ice) H$_2$ via the nuclear spin conversion on grain surfaces coated with water ice under physical conditions that are relevant to star- and planet-forming regions. We utilize the rate equation model that considers adsorption of gaseous H$_2$ on grain surfaces which have a variety of binding sites with a different potential energy depth, thermal hopping, desorption, and the nuclear spin conversion of adsorbed H$_2$. It is found that the spin conversion efficiency depends on the H$_2$ gas density and the surface temperature. As a general trend, enhanced H$_2$ gas density reduces the efficiency, while the temperature dependence is not monotonic; there is a critical surface temperature at which the efficiency is the maximum. At low temperatures, the exchange of gaseous and icy H$_2$ is inefficient (i.e., adsorbed H$_2$ does not desorb and hinders another gaseous H$_2$ to be adsorbed), while at warm temperatures, the residence time of H$_2$ on surfaces is too short for the spin conversion. Additionally, the spin conversion becomes more efficient with lowering the activation barriers for thermal hopping. We discuss whether the spin conversion on surfaces can dominate over that in the gas-phase in star- and planet-forming regions. Finally, we establish a simple but accurate way to implement the H$_2$ spin conversion on grain surfaces in existing astrochemical models.

Tables

Table 1. Table 1: η op subscript 𝜂 op \eta_{\rm op} evaluated by Equation A11

	$f=0.5$				w/o thermal hopping
	104 cm-3	106 cm-3	108 cm-3		104 cm-3	106 cm-3	108 cm-3
8 K	9.998(-1)	9.959(-1)	9.781(-1)		9.493(-1)	8.563(-1)	7.181(-1)
10 K	9.972(-1)	9.906(-1)	7.948(-1)		7.162(-1)	5.839(-1)	2.687(-1)
12 K	9.959(-1)	9.879(-1)	7.164(-1)		4.427(-1)	3.687(-1)	1.907(-1)
14 K	9.942(-1)	9.409(-1)	3.912(-1)		1.841(-1)	1.568(-1)	8.188(-2)
16 K	5.227(-1)	4.815(-1)	1.266(-1)		5.443(-2)	5.256(-2)	2.984(-2)
18 K	1.590(-2)	1.586(-2)	1.328(-2)		9.187(-3)	9.173(-3)	8.052(-3)
20 K	5.646(-4)	5.646(-4)	5.613(-4)		5.484(-4)	5.484(-4)	5.453(-4)

Equations

$$\frac{1}{4}(1-\theta (E_{\mathrm{b}}))S{v}_{\mathrm{th}}{n}_{\mathrm{site}}^{-1}n({\mathrm{H}}_2)=\nu \theta (E_{\mathrm{b}})\exp (-E_{\mathrm{b}}/T),$$

$$\theta (E_{\mathrm{b}},t)={\theta }_o(E_{\mathrm{b}},t)+{\theta }_p(E_{\mathrm{b}},t).$$

$$\int g(E_{\mathrm{b}})d{E}_{\mathrm{b}}=1.$$

$$\Theta (t)=\int \theta (E_{\mathrm{b}},t)g(E_{\mathrm{b}})d{E}_{\mathrm{b}}.$$

$${\Theta }_{\alpha }(t)=\int {\theta }_{\alpha }(E_{\mathrm{b}},t)g(E_{\mathrm{b}})d{E}_{\mathrm{b}},$$

$$\frac{dn(\alpha \text{mathchar}45\text{ce}H2)}{dt}$$

$$\frac{d{\theta }_{\alpha }(E_{\mathrm{b}},t)}{dt}$$

$$-\int k_{\mathrm{hop}}(E_{\mathrm{b}}\to E_{\mathrm{b}}^{\prime }){\theta }_{\alpha }(E_{\mathrm{b}},t)[1-\theta (E_{\mathrm{b}}^{\prime },t)]g(E_{\mathrm{b}}^{\prime })d{E}_{\mathrm{b}}^{\prime }$$

$$+[1-\theta (E_{\mathrm{b}},t)]\int k_{\mathrm{hop}}(E_{\mathrm{b}}^{\prime }\to E_{\mathrm{b}}){\theta }_{\alpha }(E_{\mathrm{b}}^{\prime },t)g(E_{\mathrm{b}}^{\prime })d{E}_{\mathrm{b}}^{\prime }$$

$$+k_{\beta \alpha }^{surf}{\theta }_{\beta }(E_{\mathrm{b}},t)-k_{\alpha \beta }^{surf}{\theta }_{\alpha }(E_{\mathrm{b}},t),$$

$$R_{\mathrm{col}}(\alpha \text{mathchar}45\text{ce}H2)$$

$$R_{\mathrm{thdes}}(\alpha \text{mathchar}45\text{ce}H2)$$

$$k_{\mathrm{thdes}}=\nu \exp (-E_{\mathrm{b}}/T).$$

$$E_{\mathrm{hop}}(E_{\mathrm{b}}\to E_{\mathrm{b}}^{\prime })=f\times \min (E_{\mathrm{b}},E_{\mathrm{b}}^{\prime })+\max (0,E_{\mathrm{b}}-E_{\mathrm{b}}^{\prime }),$$

$$k_{\mathrm{op}}^{\mathrm{surf}}$$

$$k_{\mathrm{po}}^{\mathrm{surf}}$$

$$\theta (E_{\mathrm{b}})$$

$$E_{\mathrm{b}}^{\mathrm{des}}$$

$${\tau }_{op}=n(\mathrm{o}\text{mathchar}45{\mathrm{H}}_2)/[R_{\mathrm{thdes}}(\mathrm{p}\text{mathchar}45{\mathrm{H}}_2)-R_{\mathrm{ads}}(\mathrm{p}\text{mathchar}45{\mathrm{H}}_2)],$$

$$k_{\mathrm{des}}^{\mathrm{av}}\Theta N_{\mathrm{site}}=(1-\Theta )Sn({\mathrm{H}}_2)v_{\mathrm{th}}\sigma ,$$

$$\frac{dn(\mathrm{o}\text{mathchar}45{\mathrm{H}}_2)}{dt}$$

$$\frac{dn(\mathrm{p}\text{mathchar}45{\mathrm{H}}_2)}{dt}$$

$${\eta }_{\mathrm{op}}^{FS13}=\frac{k_{\mathrm{op}}^{\mathrm{surf}}}{k_{\mathrm{op}}^{\mathrm{surf}}+k_{\mathrm{thdes}}},$$

$$\frac{(1-\Theta )S[(1-{\eta }_{\mathrm{op}})f_{\mathrm{o}}+{\eta }_{\mathrm{po}}{f}_{\mathrm{p}}]}{(1-\Theta )S[{\eta }_{\mathrm{op}}{f}_{\mathrm{o}}+(1-{\eta }_{\mathrm{po}})f_{\mathrm{p}}]},$$

$$p_{\mathrm{op}}=\frac{k_{\mathrm{op}}^{\mathrm{surf}}}{k_{\mathrm{op}}^{\mathrm{surf}}+k_{\mathrm{des}}^{\mathrm{av}}}\cdot \frac{k_{\mathrm{des}}^{\mathrm{av}}}{k_{\mathrm{po}}^{\mathrm{surf}}+k_{\mathrm{des}}^{\mathrm{av}}},$$

$$r=\frac{k_{\mathrm{op}}^{\mathrm{surf}}}{k_{\mathrm{op}}^{\mathrm{surf}}+k_{\mathrm{des}}^{\mathrm{av}}}\cdot \frac{k_{\mathrm{po}}^{\mathrm{surf}}}{k_{\mathrm{po}}^{\mathrm{surf}}+k_{\mathrm{des}}^{\mathrm{av}}}.$$

$${\eta }_{\mathrm{op}}^{(f)}$$

$$=p_{\mathrm{op}}/(1-r),$$

$$=\frac{k_{\mathrm{conv}}^{\mathrm{surf}}}{k_{\mathrm{conv}}^{\mathrm{surf}}+k_{\mathrm{des}}^{\mathrm{av}}}\cdot \frac{1}{1+\gamma },$$

$${\eta }_{\mathrm{po}}^{(f)}=\frac{k_{\mathrm{conv}}^{\mathrm{surf}}}{k_{\mathrm{conv}}^{\mathrm{surf}}+k_{\mathrm{des}}^{\mathrm{av}}}\cdot \frac{\gamma }{1+\gamma }.$$

$${\eta }_{\mathrm{op}}^{(s)}(E_{\mathrm{b}})$$

$${\eta }_{\mathrm{po}}^{(s)}(E_{\mathrm{b}})$$

$$⟨{\eta }_{\alpha \beta }^{(s)}⟩$$

$$k_{\mathrm{conv}}^{\mathrm{surf}}=[1-\theta (E_{\mathrm{b}}^{\mathrm{op}})]k_{\mathrm{hop}}(E_{\mathrm{b}}^{\mathrm{thresh}}\to E_{\mathrm{b}}^{\mathrm{op}}),$$

$${\Theta }_{\max }^{(f)}={\int }_0^{E_{\mathrm{b}}^{\mathrm{thresh}}}g(E_{\mathrm{b}}^{\prime })d{E}_{\mathrm{b}}^{\prime },{\Theta }_{\max }^{(s)}={\int }_{E_{\mathrm{b}}^{\mathrm{thresh}}}^{\infty }g(E_{\mathrm{b}}^{\prime })d{E}_{\mathrm{b}}^{\prime },$$

$${\Theta }^{(f)}={\int }_0^{E_{\mathrm{b}}^{\mathrm{thresh}}}\theta (E_{\mathrm{b}}^{\prime })g(E_{\mathrm{b}}^{\prime })d{E}_{\mathrm{b}}^{\prime },{\Theta }^{(s)}={\int }_{E_{\mathrm{b}}^{\mathrm{thresh}}}^{\infty }\theta (E_{\mathrm{b}}^{\prime })g(E_{\mathrm{b}}^{\prime })d{E}_{\mathrm{b}}^{\prime },$$

Full text

editorials, notices — miscellaneous — catalogs — surveys

††journal: ApJ

1 Introduction

Hydrogen is the most abundant element in the universe. In star- and planet-forming regions, hydrogen is primarily present in \ceH2, which has two nuclear spin configurations, ortho and para. As the internal energy difference between ortho-\ceH2 and para-\ceH2 (170.5 K) is much higher than the typical temperature of star-forming regions ($\sim$10 K), the ortho-to-para ratio (OPR) of \ceH2 can significantly affect the molecular evolution, for example, deuterium fractionation (see, e.g., Pagani et al., 1992; Flower et al., 2006; Taquet et al., 2014; Furuya et al., 2016).

\ce

H2 molecules form on grain surfaces with the statistical ortho-to-para ratio of three (Watanabe et al., 2010). After the \ceH2 formation, the ortho-para spin conversion of \ceH2 proceeds through proton exchange reactions with \ceH+ and/or with \ceH3+ in the gas phase (Gerlich, 1990; Honvault et al., 2011). Laboratory experiments have found that the \ceH2 spin conversion can also occur on bare grain (\ceD2 on graphite surfaces; e.g., Yucel et al., 1990) and on amorphous water ice surfaces (e.g., Watanabe et al., 2010) in laboratory timescales (around a few hours), while the mechanism of the spin conversion on the surfaces is not fully understood (see, e.g., Fukutani & Sugimoto, 2013; Ilisca, 2018). Given this very short timescale, it is expected that the spin conversion on surfaces affects the \ceH2 OPR evolution in star- and planet-forming regions. However, its efficiency in the astronomical conditions remains unclear for the following two reasons. First, almost all \ceH2 is present in the gas phase rather than on grain surfaces. Then the spin conversion timescale of overall (gas + solid) \ceH2 via the conversion on surfaces depends on how efficiently gaseous and solid \ceH2 interact. Second, the probability for the nuclear spin state of an adsorbed \ceH2 molecule to be changed before it is desorbed depends on the residence time on surfaces (i.e., thermal desorption timescale) versus the spin conversion timescale. Interstellar dust grains are coated with ice mantles, the main component of which is water, in the cold ($\lesssim$100 K) gas of star-forming regions (see Boogert et al., 2015, for a recent review). The surface of the ice mantles would contain various binding sites with a different energy depth. This is relevant to both points, because in that case, the thermal desorption timescale depends on site.

In order to see the two points raised above more quantitatively, first, let us consider the balance between the adsorption rate of gaseous \ceH2 on water ice surfaces and the thermal desorption rate of adsorbed \ceH2;

[TABLE]

where $S$ is the sticking probability to the water ice surface, $n({\rm H_{2}})$ is the number density of \ceH2 in the gas phase, $v_{\rm th}$ is the thermal velocity of \ceH2, $n_{\rm site}$ is the density of binding sites on the surface ($1.5\times 10^{15}$ cm*-2*), $\nu$ is the vibrational frequency (typically 1012 s*-1*), $E_{\rm b}$ is the binding energy of \ceH2 on the water ice surface, and $T$ is the temperature of the surface. $\theta(E_{\rm b})$ is the fraction of binding sites occupied by \ceH2 with the potential energy depth of $E_{\rm b}$. We assume that only one \ceH2 is allowed per binding site, which leads to the factor $1-\theta$ in the left hand side of the equation. From this equation, we can define critical binding energy ($E_{\rm b}^{\rm\,des}$) such that all sites with $E_{\rm b}>E_{\rm b}^{\rm\,des}$ will be occupied by \ceH2, i.e., $\theta(E_{\rm b}>E_{\rm b}^{\rm\,des})=1$ (Dissly et al., 1994). At $T=10$ K and $n({\rm H_{2}})=10^{4}$ cm*-3*, $E_{\rm b}^{\rm\,des}$ is 440 K. Let us define another critical binding energy ($E_{\rm b}^{\rm\,op}$) such that the thermal desorption timescale of \ceH2 in binding sites with $E_{\rm b}>E_{\rm b}^{\rm\,op}$ is long enough for the conversion of ortho-\ceH2 to para-\ceH2. By considering the balance between the spin conversion rate of ortho-\ceH2 to para-\ceH2 ($k_{\rm op}^{\rm surf}$) and the thermal desorption rate ($\nu\exp(-E_{\rm b}/T)$), we obtain $E_{\rm b}^{\rm\,op}$ of $\sim$360 K at the surface temperature of 10 K and for $k_{\rm op}^{\rm surf}$ of $3\times 10^{-4}$ s*-1* (Ueta et al., 2016). Based on theses arguments, one may think that binding sites which satisfies $E_{\rm b}^{\rm\,op}\lesssim E_{\rm b}\lesssim E_{\rm b}^{\rm\,des}$ contribute to the evolution of the \ceH2 OPR most efficiently; for $E_{\rm b}<E_{\rm b}^{\rm\,op}$, the residence time is too short for the spin conversion, while for $E_{\rm b}^{\rm\,des}<E_{\rm b}$, adsorbed \ceH2 does not desorb efficiently and hinders another gaseous \ceH2 to be adsorbed. Then binding energy distribution does matter, and the question is what fraction of sites have binding energy in the range of $E_{\rm b}^{\rm\,op}\lesssim E_{\rm b}\lesssim E_{\rm b}^{\rm\,des}$. Note that $E_{\rm b}^{\rm\,op}$, $E_{\rm b}^{\rm\,des}$, and their inequality relation depend on physical conditions as shown in Figure 1. The above discussion, however, neglects thermal hopping of adsorbed \ceH2. As we will see later, thermal hopping changes the situation significantly, because it allows adsorbed \ceH2 molecules to visit various sites with a different potential energy depth. In summary, to understand the spin conversion efficiency of \ceH2 in the astronomical conditions, one has to consider \ceH2 adsorption on the surface, which contains a variety of sites, thermal desorption and hopping, and the nuclear spin conversion in a self-consistent way. In this work, we construct such a model for the first time.

The effect of the spin conversion on the surface on the \ceH2 OPR evolution was theoretically studied by Bron et al. (2016) in the context of photodissociation regions (PDRs) and by Bovino et al. (2017) in the context of dense molecular clouds. Both models did not consider the binding energy distribution of \ceH2, but used a single “representative” binding energy as commonly assumed in astrochemical models for simplicity. Bron et al. (2016) found that the fluctuation of dust temperature due to stochastic heating by UV photons is important for determining the spin conversion efficiency in PDRs. Bovino et al. (2017) discussed $E_{\rm b}^{\rm\,op}$ considering uncertainties of relevant parameters, but their discussion lacks another necessary condition, $E_{\rm b}^{\rm\,des}$.

This paper is organized as follows: our numerical model is described in Sect. 2 and the results are discussed in in Sect. 3. In Sect. 4, we propose a simple model that reproduces our numerical results, and discuss whether the spin conversion on surfaces can dominate over that in the gas phase in star- and planet-forming regions. Our findings are summarized in Sect. 5.

2 Methods

2.1 Basic equations

We adopt a rate equation approach to investigate the efficiency of the ortho-para conversion on grain surfaces in star and planet forming regions. We consider a typical interstellar grain with radius of 0.1 $\mu$m with the dust-to-gas mass ratio of $10^{-2}$. The grain is assumed to be covered by water ice mantles and the number of binding sites on the water ice surface per area ($n_{\rm site}$) is set to be $1.5\times 10^{15}$ cm*-2*. The total number of binding sites per grain is thus $N_{\rm site}\approx 2\times 10^{6}$. We consider the experimental fact that the water ice surface contains various sites with a different potential energy depth (e.g., Amiaud et al., 2006). For simplicity, we assume the following throughout this work: (1) only one molecule is allowed to be adsorbed per binding site, (2) ortho-\ceH2 (${\rm o\mathchar 45\relax H_{2}}$) and para-\ceH2 (${\rm p\mathchar 45\relax H_{2}}$) share common binding sites, following the same binding energy distribution (see Section 2.2), and (3) chemical properties of ${\rm o\mathchar 45\relax H_{2}}$ and ${\rm p\mathchar 45\relax H_{2}}$ are the same except that they convert to each other on the surface with different rates. We denote the fraction of biding sites, which are occupied by ${\rm o\mathchar 45\relax H_{2}}$ (${\rm p\mathchar 45\relax H_{2}}$) as $\theta_{o}$ ($\theta_{p}$). The following condition should be satisfied:

[TABLE]

We denote the binding energy distribution of \ceH2 on the surface as $g$, which satisfies

[TABLE]

The surface coverage of \ceH2 at a given time $t$, $\Theta(t)$, is defined as

[TABLE]

Similarly $\Theta_{\alpha}(t)$, where $\alpha$ is $o$ or $p$, is defined as

[TABLE]

and thus $\Theta(t)=\Theta_{o}(t)+\Theta_{p}(t)$.

We numerically solve the following rate equations, which describe adsorption of \ceH2, thermal desorption, thermal hopping, and spin conversion of adsorbed \ceH2, considering various binding sites with a different potential energy depth (cf. Li et al., 2010):

[TABLE]

where $\alpha$ and $\beta$ indicate either ortho ($o$) or para ($p$). The collision rates to dust grains and desorption rates from the whole surface of dust grains of ${\rm o\mathchar 45\relax H_{2}}$ and ${\rm p\mathchar 45\relax H_{2}}$ are given by

[TABLE]

where $v_{\rm th}$ is the thermal velocity, $\sigma$ is the cross section of a dust grain, $n_{\rm gr}$ is the number density of dust grains per unit gas volume, and $k_{\rm thdes}$ is the thermal desorption rate (s*-1*). We assume gas and surface temperatures are the same and do not distinguish them throughout this paper.

The first terms in Eqs. 6 and 7 represent adsorption with the sticking probability $S$ of \ceH2 to the water ice surface. We consider the factor $1-\Theta$ or the factor $1-\theta$, because only one molecule is allowed to be adsorbed per binding site in our models. Then the maximum value of $\Theta$ is unity and the formation of \ceH2 multilayers does not occur in our models. Indeed, laboratory experiments have found that no matter how large \ceH2 fluence deposited on a water ice substrate is, the \ceH2 coverage is in the submonolayer regime even at 10 K (e.g., Gavilan et al., 2012; Kuwahata et al., 2015). The second terms in Eqs. 6 and 7 represent thermal desorption. The third term in Eq. 7 represents thermal hopping from a binding site with $E_{\rm b}$ to a site with $E_{\rm b}^{\prime}$, while the forth term represents the reverse process. The hopping activation energy in our models is discussed later. The fifth and sixth terms are for ortho-para conversion on surfaces, the rates of which are discussed in Section 2.3. Initially, all \ceH2 are assumed to be present in the gas phase with the ${\rm OPR(\ce{H2})}$ of three (i.e., the statistical value).

2.2 Binding energy distribution, hopping activation energy, and sticking probability

We use binding energy distribution and sticking probability of \ceH2 that are appropriate for nonporous amorphous solid water (ASW) in this work. The degree of porosity of interstellar ice, which is mainly composed of water, remains unclear. There is no clear observational evidence that interstellar water ice has a porous structure; the OH dangling bonds of water ice have not been detected in the midinfrared spectrum in the interstellar matter (ISM; Keane et al., 2001), although the non-detection might be due to the sensitivity limitations of the Infrared Space Observatory (ISO). Oba et al. (2009) found that in their experiments, water ices formed from atomic hydrogen and molecular oxygen at low temperatures (10 K-40 K) present a nonporous structure compared to vapor deposited water ices at the low temperatures. Garrod (2013) found that in their off-lattice Monte-Carlo simulations, ices formed by surface chemistry under dark cloud conditions present a nonporous structure, being consistent with the experiments. In addition, laboratory experiments have found that the porosity of amorphous water ice decreases after UV photon irradiation and/or cosmic-ray impacts (e.g., Raut et al., 2008; Palumbo et al., 2010). Taken together, nonporous ASW could be more representative for interstellar water ice rather than porous ASW.

The thermal desorption rate depends on the binding energy of the species to the surface,

[TABLE]

Our binding energy distribution of \ceH2, which ranges from 290 K to 635 K, is divided into 100 equal intervals in our simulations (see black line in the top panel of Figure 2). The binding energy distribution of \ceD2 on nonporous ASW is available in the literature, which was obtained from direct inversion of temperature programmed desorption spectra (Amiaud et al., 2007; He & Vidali, 2014). We obtained the binding energy distribution of \ceH2 by considering zero-point energy difference between \ceD2 and \ceH2, 3.15 meV (Amiaud et al., 2015). Note that laboratory experiments have found that ${\rm o\mathchar 45\relax D_{2}}$ is bound to surfaces slightly more strongly than ${\rm p\mathchar 45\relax D_{2}}$ ($\sim$1 meV) (Amiaud et al., 2008; Tsuge et al., 2019), but we neglect the difference in this work for simplicity.

The hopping activation energy from a site with the binding energy of $E_{\rm b}$ to another site with the binding energy of $E_{\rm b}^{\prime}$ ($E_{\rm hop}$) is given as follows (Cazaux et al., 2017, see their Fig. 11):

[TABLE]

where $f$ is a free parameter. The parameter $f$, which is the hopping-to-binding energy ratio, is poorly constrained and values between 0.3 and 0.8 are normally assumed in the astrochemical community. We choose $f=0.5$ in our fiducial model. Given the expression of $E_{\rm hop}$, the thermal hopping rate, $k_{\rm hop}=\nu\exp(-E_{\rm hop}/T)$, obeys the microscopic reversibility, i.e., $k_{\rm hop}(E_{\rm b}\rightarrow E_{\rm b}^{\prime})/k_{\rm hop}(E_{\rm b}^{\prime}\rightarrow E_{\rm b})=\exp[-(E_{\rm b}-E_{\rm b}^{\prime})/T]$ (Cuppen et al., 2013).

He et al. (2016) experimentally investigated the sticking probability for stable molecules on nonporous ASW in low surface coverage regime (below 10 %). For sticking probability of \ceH2 onto water ice surfaces ($S$), we use the formula recommended by He et al. (2016) (see their Eq. 1). For example, $S$ is $\sim$0.7 at 10 K and $\sim$0.5 at 16 K. The experimental values may be considered as the surface averaged value, while the sticking probability for each site may depend on the energy depth of each site. Such (possible) complexity is not considered in our models, i.e, $S$ is set to be the same for all binding sites.

2.3 Ortho-para conversion rates on surfaces

The ortho-para conversion timescale of \ceH2 on amorphous water ice ($\tau^{\rm surf}_{\rm conv}$) in the temperature range between 9 K and 16 K was measured in laboratory by Ueta et al. (2016). From $\tau^{\rm surf}_{\rm conv}$, the rate of the conversion from ${\rm o\mathchar 45\relax H_{2}}$ to ${\rm p\mathchar 45\relax H_{2}}$ ($k_{\rm op}^{\rm surf}$) and that of the reverse process ($k_{\rm po}^{\rm surf}$) can be deduced to be

[TABLE]

where $\gamma$ is the thermalized value of ${\rm OPR(\ce{H2})}$, $9\exp(-170.5/T)$ (Bron et al., 2016), assuming the energy difference between ${\rm o\mathchar 45\relax H_{2}}$ and ${\rm p\mathchar 45\relax H_{2}}$ on water ice surfaces is the same as that in the gas phase. On water ice surfaces, \ceH2 molecules would not rotate freely and thus the energy difference between ${\rm o\mathchar 45\relax H_{2}}$ and ${\rm p\mathchar 45\relax H_{2}}$ would become smaller than that in the gas phase, but the exact value remains unclear (cf. see Hama et al., 2016, for the discussion on the energy difference between ortho-\ceH2O and para-\ceH2O on surfaces). Ueta et al. (2016) found that at the temperature lower than $\sim$12 K, $\tau_{\rm conv}$ is fitted by a power low $1/(AT^{n})$, where $A$ is $3.2\times 10^{-11}$ s*-1* and $n$ is 7.1. At the higher temperature, $\tau^{\rm surf}_{\rm conv}$ is almost constant with the value of around $1/(1.5\times 10^{-3})\approx 670$ s. We take $\tau^{\rm surf}_{\rm conv}$ from Ueta et al. (2016) with the lower limit of 670 s. At 10 K, for example, $k_{\rm op}^{\rm surf}$ and $k_{\rm po}^{\rm surf}$ are $3.1\times 10^{-4}$ s*-1* and $10^{-12}$ s*-1*, respectively.

3 Results

3.1 \ceH2 coverage

The \ceH2 coverage on the water ice surface is discussed in detail in our separate work (Furuya et al., in prep.), where the similar rate equations to Eqs. 6 and 7 are used, but without distinction of the \ceH2 nuclear spin states. We briefly summarize this here. The adsorption and desorption of \ceH2 reach the equilibrium in a very short timescale ($\lesssim$1 yr ($10^{4}$ cm*-3*/$n({\rm H_{2}})$)). Then only the equilibrium condition is relevant in the dense ISM. The occupation of sites with the potential energy depth of $E_{\rm b}$ is then determined by the balance between the adsorption rate of gaseous \ceH2 on each site and the thermal desorption rate of adsorbed \ceH2, and one can obtain $\theta(E_{\rm b})$ by solving Eq. 1 (see also Amiaud et al., 2006);

[TABLE]

$E_{\rm b}^{\rm\,des}$ is the critical binding energy such that a half of sites with $E_{\rm b}$ will be occupied by \ceH2, i.e., $\theta(E_{\rm b}^{\rm\,des})=0.5$, under the adsorption-desorption equilibrium. Thus all sites with $E_{\rm b}\gg E_{\rm b}^{\rm\,des}$ will be occupied by \ceH2, $\theta(E_{\rm b}\gg E_{\rm b}^{\rm\,des})=1$. As an example, the top panel of Figure 2 shows the equilibrium occupation distribution ($\theta\times g$) of \ceH2 at $n({\rm H_{2}})=10^{4}$ cm*-3* and $T=10$ K, where the \ceH2 surface coverage ($\Theta$) is $\sim$30 %. It shows that deeper sites are preferentially occupied by \ceH2. The equilibrium \ceH2 coverage as functions of the \ceH2 gas density and temperature are shown in Figure 3. The \ceH2 coverage increases with increasing the gas density and with decreasing the temperature. The occupation distribution and the \ceH2 surface coverage at the equilibrium do not depend on the hopping parameter $f$.

3.2 \ceH2 ortho-para spin conversion in the fiducial physical conditions

We first show model results at $n({\rm H_{2}})=10^{4}$ cm*-3* and $T=10$ K (our fiducial physical conditions) and discuss the dependence on the physical conditions later. The middle panel of Figure 2 shows the ${\rm OPR(\ce{H2})}$ on the surface as functions of $E_{\rm b}$ (i.e., $\theta_{o}/\theta_{p}$ for each $E_{\rm b}$), varying the parameter $f$. We chose the time when the \ceH2 fluence (the time integral of the \ceH2 flux) reaches $5\times 10^{16}$ cm*-2*, corresponding to the duration time of $\sim 10(10^{4}$ cm*-3*/$n({\rm H_{2}})$) yr. By that time, the \ceH2 coverage on the surface reaches the adsorption-desorption equilibrium at all physical conditions explored in this work, while the duration time is too short for the spin conversion of the overall (gas + solid) ${\rm OPR(\ce{H2})}$. Then the ${\rm OPR(\ce{H2})}$ in the gas phase remains unchanged from the initial value of three. The chosen duration time is shorter than $\tau^{\rm surf}_{\rm conv}$ for $n({\rm H_{2}})\geq 10^{9}$ cm*-3*; then we choose $t\sim$10*-3* yr ($\gg\tau^{\rm surf}_{\rm conv}$) at the higher densities, corresponding to the \ceH2 fluence of $5\times 10^{16}$ ($n({\rm H_{2}})/10^{8}$ cm*-3*) cm*-2*.

When the thermal hopping is turned off, $\theta_{o}/\theta_{p}$ is determined by the timescale of thermal desorption from the site ($k_{\rm thdes}^{-1}$) versus the spin conversion timescale; sites with higher $E_{\rm b}$ have lower $\theta_{o}/\theta_{p}$ due to the longer thermal desorption timescale (i.e., the longer resident timescale). When the thermal hopping is turned on, the situation changes; adsorbed \ceH2 can visit multiple sites via thermal hopping. In the fast hopping cases ($f\leq 0.5$), $\theta_{o}/\theta_{p}$ is almost constant across the surface. This indicates that the resident time of \ceH2 on the surface is independent of the energy depth of a site in which a \ceH2 molecule was initially adsorbed, due to the efficient thermal hopping after adsorption on the surface. Then $k_{\rm thdes}^{-1}$ is not a good measure of the resident time of adsorbed \ceH2, when thermal hopping is considered. While $\theta_{o}/\theta_{p}$ for given $E_{\rm b}$ is very different depending on the hopping rate, the OPR averaged on the whole surface ($\Theta_{o}/\Theta_{p}$) is similar regardless of the hopping rate; e.g., $\Theta_{o}/\Theta_{p}=1.7\times 10^{-4}$ in the model with $f=0.5$, while it is $1.2\times 10^{-4}$ in the model without thermal hopping.

In terms of the ${\rm OPR(\ce{H2})}$ evolution in the ISM, the ortho-para ratio of desorbing gas (i.e., $R_{\rm thdes}({\rm o\mathchar 45\relax H_{2}})/R_{\rm thdes}({\rm p\mathchar 45\relax H_{2}})$) is more relevant than that of the surface ($\Theta_{o}/\Theta_{p}$), as almost all \ceH2 is present in the gas phase rather than on the surface. We find that relation between $R_{\rm thdes}({\rm o\mathchar 45\relax H_{2}})/R_{\rm thdes}({\rm p\mathchar 45\relax H_{2}})$ and $\Theta_{o}/\Theta_{p}$ depends on the efficiency of thermal hopping; they are similar in the models with fast hopping ($1.9\times 10^{-4}$ versus $1.7\times 10^{-4}$ for $f=0.5$), while they are very different in the model without thermal hopping ($2.7\times 10^{-1}$ versus $1.2\times 10^{-4}$). Thermal desorption rates of ${\rm o\mathchar 45\relax H_{2}}$ and ${\rm p\mathchar 45\relax H_{2}}$ as functions of $E_{\rm b}$ are shown in the bottom panel of Figure 2. \ceH2 desorption predominantly occurs from sites with $E_{\rm b}\lesssim E_{\rm b}^{\rm\,des}$ (by definition). In the model without thermal hopping, only sites with $E_{\rm b}^{\rm\,op}\lesssim E_{\rm b}\lesssim E_{\rm b}^{\rm\,des}$ contribute to the decrease of $R_{\rm thdes}({\rm o\mathchar 45\relax H_{2}})/R_{\rm thdes}({\rm p\mathchar 45\relax H_{2}})$; for $E_{\rm b}<E_{\rm b}^{\rm\,op}$, the residence time is too short for the spin conversion, while for $E_{\rm b}^{\rm\,des}<E_{\rm b}$, adsorbed \ceH2 does not desorb efficiently. On the other hand, in the model with thermal hopping, sites with $E_{\rm b}\gtrsim E_{\rm b}^{\rm\,des}$ also contribute to the decrease of $R_{\rm thdes}({\rm o\mathchar 45\relax H_{2}})/R_{\rm thdes}({\rm p\mathchar 45\relax H_{2}})$; they trap \ceH2 molecules, the spin states of the \ceH2 molecules are converted, and after some time, the \ceH2 molecules hop to shallower sites and desorb to the gas phase. These results demonstrate that the binding energy distribution and the thermal hopping among sites are essentially important for the spin conversion on grain surfaces in the ISM.

Figure 4 shows the long term evolution of the ${\rm o\mathchar 45\relax H_{2}}$ abundance in the gas phase with respect to \ceH2 at $n({\rm H_{2}})=10^{4}$ cm*-3* and $T=10$ K. The gaseous ${\rm o\mathchar 45\relax H_{2}}$ abundance decreases with time due to the spin conversion on the surface. The spin conversion timescale of ${\rm o\mathchar 45\relax H_{2}}$ to ${\rm p\mathchar 45\relax H_{2}}$ is in the order of 105 yr, and the timescale is shorter in the model with $f=0.5$ than that in the model without hopping by a factor of $\lesssim$2. The spin conversion timescale of gaseous ${\rm o\mathchar 45\relax H_{2}}$ is given by

[TABLE]

where $R_{\rm ads}$ is the adsorption rate of \ceH2 on dust grain surfaces ($=(1-\Theta)SR_{\rm col}$). We confirmed that $x_{0}({\rm o\mathchar 45\relax H_{2}})\exp(-t/\tau_{op})$, where $x_{0}({\rm o\mathchar 45\relax H_{2}})$ is the initial abundance of gaseous ${\rm o\mathchar 45\relax H_{2}}$ with respect to \ceH2, reproduces the numerical results shown in Fig. 4. Note that if the spin conversion on the surface does not occur, $R_{\rm thdes}({\rm p\mathchar 45\relax H_{2}})=R_{\rm ads}({\rm p\mathchar 45\relax H_{2}})$ under the adsorption-desorption equilibrium. The steady-state abundance is given by $9\exp(-170.5/10)\approx 3\times 10^{-7}$.

3.3 Density dependence

Here we discuss the density dependence of the spin conversion timescale of gaseous ${\rm o\mathchar 45\relax H_{2}}$ via the spin conversion on the surface ($\tau_{op}$). Again, we focus on the results at $t\sim 10(10^{4}$ cm*-3*/$n({\rm H_{2}})$) yr for $n({\rm H_{2}})\leq 10^{8}$ cm*-3*, while at higher densities, we focus on the results at $t\sim 10^{-3}$ yr. As discussed in the Introduction, there are two main factors that control $\tau_{op}$: (i) the efficiency of the interaction between gaseous and solid \ceH2 and (ii) the probability of the spin-state conversion of an adsorbed \ceH2 molecule before it is desorbed. The top panel of Figure 5 shows the timescale for gaseous and solid \ceH2 interaction defined by $\tau_{\rm int}=n({\rm H_{2}})/R_{\rm thdes}(\ce{H2})$ (or equivalently $n({\rm H_{2}})/R_{\rm ads}(\ce{H2})$ at the adsorption-desorption equilibrium) as function of the \ceH2 gas density. For convenience, we normalize $\tau_{\rm int}$ by the collisional timescale, $\tau_{\rm col}=n({\rm H_{2}})/R_{\rm col}(\ce{H2})=(v_{\rm th}\sigma n_{\rm gr})^{-1}\approx 3\times 10^{9}/n({\rm\ce{H2}})$ yr. The normalized interaction timescale ($\tau_{\rm int}/\tau_{\rm col}=1/[S(1-\Theta)]$), which means the average number of collisions for an \ceH2 molecule required to be adsorbed on the water ice surface, becomes larger with increasing $n({\rm H_{2}})$, because $E_{\rm b}^{\rm\,des}$ becomes smaller (see Fig. 1) and thus $\Theta$ increases.

The middle panel of Figure 5 shows the OPR(\ceH2) in the desorbing gas, $R_{\rm thdes}({\rm o\mathchar 45\relax H_{2}})/R_{\rm thdes}({\rm p\mathchar 45\relax H_{2}})$; it is higher (i.e., the spin conversion probability upon adsorption becomes lower) with increasing the \ceH2 gas density. This trend does not depend on the parameter $f$, which is explained as follows. In the model without thermal hopping, only sites with $E_{\rm b}^{\rm\,op}\lesssim E_{\rm b}\lesssim E_{\rm b}^{\rm\,des}$ efficiently contribute to the spin conversion. At a given temperature, $E_{\rm b}^{\rm\,des}$ becomes smaller with increasing the \ceH2 gas density, but $E_{\rm b}^{\rm\,op}$ does not change (see Fig. 1); the number of sites that efficiently contribute to lowering the ${\rm OPR(\ce{H2})}$ in the desorbing gas becomes smaller with increasing the \ceH2 gas density. In the case with hopping, adsorbed \ceH2 can visit various potential sites before it is desorbed, and $\theta_{o}/\theta_{p}$ is similar across the surface regardless of $E_{\rm b}$ as discussed above. Let us define the averaged desorption rate of \ceH2 ($k_{\rm des}^{\rm av}$) on the surface, which should satisfy

[TABLE]

under the adsorption-desorption equilibrium (see also the Appendix). Then the averaged residence time of an adsorbed \ceH2 on the surface, 1/$k_{\rm des}^{\rm av}$, is proportional to $\Theta/(1-\Theta)$ and inversely proportional to $n({\rm H_{2}})$. As $\Theta$ depends only weakly on $n({\rm H_{2}})$ (Fig. 3), the average residence time on the surface is reduced with increasing $n({\rm H_{2}})$. Therefore, the probability of the ortho-para conversion upon adsorption is reduced with increasing $n({\rm H_{2}})$ in the case with hopping as well.

Finally, the spin conversion timescale of gaseous ${\rm o\mathchar 45\relax H_{2}}$ ($\tau_{op}$) normalized by $\tau_{\rm col}$ is shown in the bottom panel of Figure 5. The normalized conversion timescale ($\tau_{op}/\tau_{\rm col}$) corresponds to the average number of ${\rm o\mathchar 45\relax H_{2}}$ collision to the surface to be required to produce one ${\rm p\mathchar 45\relax H_{2}}$. $\tau_{op}/\tau_{\rm col}$ increases with increasing $n({\rm H_{2}})$, because the timescale for gaseous and solid \ceH2 interaction becomes longer (the top panel) and the spin conversion upon adsorption becomes less efficient (the middle panel) with increasing $n({\rm H_{2}})$. Note that $\tau_{\rm col}$ is inversely proportional to $n({\rm H_{2}})$ and thus $\tau_{\rm op}$ in fact drops with increasing $n({\rm H_{2}})$. The absolute value of $\tau_{op}/\tau_{\rm col}$ depends on the efficiency of thermal hopping; $\tau_{op}/\tau_{\rm col}$ in the models with ($f\leq 0.5$) is smaller than that in the model without thermal hopping by a factor of a few.

3.4 Temperature dependence

Figure 6 is similar to Figure 5, but shows dependencies on temperature. At given $n({\rm H_{2}})$, the normalized interaction timescale, $\tau_{\rm int}$/$\tau_{\rm col}$, becomes smaller with increasing temperature, because $\Theta$ decreases with increasing temperature. On the other hand, the OPR(\ceH2) in the desorbing gas becomes higher (i.e., the spin conversion probability upon adsorption becomes lower) with increasing temperature. As the temperature affects the interaction timescale and the conversion probability in the opposite direction, there is a critical temperature at which $\tau_{op}/\tau_{\rm col}$ is the smallest for given $n({\rm H_{2}})$: $\sim$12-14 K for $n({\rm H_{2}})=10^{4}$ cm*-3* and $\sim$10-14 K for $n({\rm H_{2}})=10^{8}$ cm*-3*. At the lower temperatures, the exchange of gaseous and icy \ceH2 is inefficient (i.e., adsorbed \ceH2 does not desorb and hinders another gaseous \ceH2 to be adsorbed), while at the higher temperatures, the residence time of \ceH2 on surfaces is too short for the spin conversion.

Figure 7 shows the long term evolution of the ${\rm o\mathchar 45\relax H_{2}}$ abundance in the gas phase with respect to \ceH2 in the models with (left) and without thermal hopping (right) at $n({\rm H_{2}})=10^{4}$ cm*-3*, varying temperature from 8 K to 20 K. Again, the temperature dependence of the spin conversion time scale is non-monotonic. Note that the steady-state abundances depends on the temperature and are given by $9\exp(-170.5/T)$.

4 Discussion

4.1 A simple model for the \ceH2 spin conversion rate

For the spin conversion of \ceH2 on grain surfaces in the ISM, the binding energy distribution and the thermal hopping among sites are essentially important. In astrochemical simulations of star- and planet-forming regions, the rate-equation approach is usually employed to describe the gas-phase and grain-surface chemistry. In rate equation models, binding energy distribution is normally neglected and the binding energy of each species is represented as a single “representative” value. The goal of this subsection (and Appendix A) is to derive simple equations that reproduce our full numerical simulations of the overall (gas + solid) ${\rm OPR(\ce{H_{2}})}$ evolution and can be implemented in existing astrochemical codes easily.

The spin conversion rate of overall \ceH2 via the conversion on surfaces can be described as $\eta_{\rm op}(1-\Theta)SR_{\rm col}({\rm o\mathchar 45\relax H_{2}})$ and $\eta_{\rm po}(1-\Theta)SR_{\rm col}({\rm p\mathchar 45\relax H_{2}})$, where $\eta_{\rm op}$ ($\eta_{\rm po}$) is the yield of gaseous ${\rm p\mathchar 45\relax H_{2}}$ (${\rm o\mathchar 45\relax H_{2}}$) per ${\rm o\mathchar 45\relax H_{2}}$ (${\rm p\mathchar 45\relax H_{2}}$) adsorption. $\eta_{\rm op}(1-\Theta)S$ expresses the yield of ${\rm p\mathchar 45\relax H_{2}}$ per ${\rm o\mathchar 45\relax H_{2}}$ collision to dust grain surfaces, which is equivalent to $\tau_{\rm col}/\tau_{op}$. If such $\eta_{\rm op}$ and $\eta_{\rm po}$ are given, the time evolution of the \ceH2 OPR in the gas phase via the spin conversion on grain surfaces can be obtained by solving simple rate equations, assuming the adsorption-desorption equilibrium of \ceH2:

[TABLE]

Once the binding energy distribution of \ceH2 is given, it is straightforward to calculate $\Theta$ using Eqs. 4, 14, and 15.

The expression of $\eta_{\rm op}$ has been proposed by Fukutani & Sugimoto (2013) as

[TABLE]

which describes the competition between the spin conversion and thermal desorption of adsorbed \ceH2. It has been used in astrochemical simulations (Bovino et al., 2017). While Eq. 20 is valid when surface property can be described by single binding energy, it does not take into account the binding energy distribution and the thermal hopping among various sites. We develop more rigorous expression of $\eta_{\rm op}$ and $\eta_{\rm po}$, which reproduces our numerical results. Our strategy is as follows: we first construct $\eta_{\rm op}$ and $\eta_{\rm po}$ that are adequate in two extreme cases, the fast hopping case and the slow (no) hopping case, and then combine the two extremes to obtain a general expression. The derivation and formulations of $\eta_{\rm op}$ and $\eta_{\rm po}$ are described in the Appendix.

Using $\eta_{\rm op}$ and $\eta_{\rm po}$, the ${\rm OPR(\ce{H_{2}})}$ of the desorbing gas from the surface can be expressed as

[TABLE]

where $f_{\rm o}$ and $f_{\rm p}$ are the fraction of ${\rm o\mathchar 45\relax H_{2}}$ and ${\rm p\mathchar 45\relax H_{2}}$, respectively, in adsorbing \ceH2 (or equivalently \ceH2 in the gas phase). The factor $1-\eta_{\rm op}$ indicates the probability that adsorbed ${\rm o\mathchar 45\relax H_{2}}$ desorbs as ${\rm o\mathchar 45\relax H_{2}}$. We set $f_{\rm o}$ and $f_{\rm p}$ to be 0.75 and 0.25, respectively, and compare Eq. 21 with the numerical results (i.e., $R_{\rm thdes}({\rm o\mathchar 45\relax H_{2}})/R_{\rm thdes}({\rm p\mathchar 45\relax H_{2}})$) at $t\sim 10(10^{4}$ cm*-3*/$n({\rm H_{2}})$) yr or $t\sim 10^{-3}$ yr, whichever is longer) in Figure 8. In the case where the thermal hopping is efficient ($f\leq 0.5$) or where thermal hopping is turned off, the ${\rm OPR(\ce{H_{2}})}$ of the desorbing gas obtained by Eq. 21 with our $\eta_{\rm op}$ and $\eta_{\rm po}$ almost perfectly agrees with the results of the full numerical simulation. In the case of $f=0.8$, the two results are deviate, but only by a factor of two at most. In the bottom panel of Figure 8, Eq. 21 evaluated with Eq. 20 and $k_{\rm thdes}=\nu\exp(-440/T)$ is also shown (gray dashed line). We chose 440 K as “representative” binding energy of \ceH2 (e.g., Cuppen& Herbst, 2007). In this case, the ${\rm OPR(\ce{H_{2}})}$ starts to sharply drop at $\sim$13 K, where $E_{\rm b}^{\rm\,op}\sim 440$ K (see Fig. 1), because of the exponential dependence of $k_{\rm thdes}$. The comparison clearly demonstrates that Eq. 20 does not reproduce our numerical results.

4.2 Conversion on the surface versus in the gas phase

So far, we investigated the efficiency of the \ceH2 spin conversion on grain surfaces by solving the rate equations of gas-phase and grain-surface \ceH2. Our detailed modeling has revealed that the efficiency of the \ceH2 spin conversion on grain surfaces depends on the temperature, the \ceH2 gas density, and the thermal hopping rates. In the ISM, the ${\rm OPR(\ce{H_{2}})}$ is mostly determined by the competition between the \ceH2 formation on surfaces and the spin conversion (i.e., thermalization) in the gas phase and that on surfaces. In the dense ISM, hydrogen is predominantly present in molecular form, and only small fraction of hydrogen is in atomic form, which is produced via a sequence of gas-phase reactions initiated by the cosmic-ray ionization of \ceH2 (e.g., Tielens, 2005). Atomic hydrogen can recombine on grain surfaces to reform \ceH2. The OPR of \ceH2 upon formation on surfaces is three, while the thermalized value of the OPR is on the order of 10*-7* at 10 K. This significant deference makes the \ceH2 formation important for the \ceH2 OPR evolution in the dense ISM, even if the rate of \ceH2 formation is much lower than the spin conversion rate (see e.g., Furuya, 2018). Here we investigate the evolution of the ${\rm OPR(\ce{H_{2}})}$ in full gas-grain chemical reaction network model, in which the three relevant processes, the \ceH2 formation and the spin conversion in the gas phase and on grain surfaces, are considered. The main question we would like to explore here is at which conditions the spin conversion on grain surfaces dominates over the conversion in the gas phase.

We run a grid of full gas-grain chemical reaction network model, which includes a variety of gaseous and icy species in addition to \ceH2. The model is based on Furuya et al. (2015), but additionally considers the spin conversion on grain surfaces using our $\eta_{op}$ and $\eta_{po}$. In the model of Furuya et al. (2015), the gas-ice chemistry is described by a three-phase model, in which three distinct phases, gas-phase, icy grain surface, and the bulk of ice mantle are considered (Hasegawa & Herbst, 1993). Gas-phase reactions, gas-surface interactions, and surface reactions are considered. The chemical reaction network includes nuclear spin states of \ceH2 and \ceH+3 and deuterated species. The \ceH2 spin conversion in the gas phase through proton exchange reactions with \ceH+ and with \ceH3+ is included (Gerlich, 1990; Honvault et al., 2011). For this work, we exclude deuterated species for simplicity.

We run a grid of pseudo-time dependent models (i.e., the gas density and the temperature are fixed in each model), varying $n({\rm H_{2}})$ from 104 cm*-3* to 108 cm*-3* and temperature from 8 K to 20 K. For each physical condition, we run three models, varying the treatment of the \ceH2 spin conversion on grain surfaces: the model without the conversion on the surfaces, the model in which $\eta_{op}$ and $\eta_{po}$ are calculated assuming $f=0.5$, and the model in which $\eta_{op}$ and $\eta_{po}$ are calculated neglecting thermal hopping. $\Theta$ is calculated using Eqs. 4 and 14, and the sticking probability of \ceH2 is taken from He et al. (2016). We assume uniform grain radius of 0.1 $\mu$m with the dust-to-gas mass ratio of $10^{-2}$. Elemental abundance ratios for H:He:C:N:O:Na:Mg:Si:S:Fe are 1.00:9.75(-2):7.86(-5):2.47(-5):1.80(-4):2.25(-9):1.09(-8):9.74(-9):9.14(-8):2.74(-9), where $a(-b)$ means $a\times 10^{-b}$ (Aikawa & Herbst, 1999). Initially, all elements except for hydrogen is in atomic form, while hydrogen is present as \ceH2 with the OPR of three. The cosmic-ray ionization rate of \ceH2 is set to be $\xi=1.3\times 10^{-17}$ s*-1*.

Figure 9 shows the temporal variation of the ${\rm OPR(\ce{H2})}$ in the gas phase. The steady-state value of the ${\rm OPR(\ce{H2})}$ in our models is higher than the thermalized value at $\lesssim$16 K due to the \ceH2 formation on grain surfaces (e.g., Furuya, 2018). We find that the conversion on the surfaces dominates over that in the gas phase at the temperatures below 20 K, regardless of the \ceH2 gas density and the thermal hopping rates of \ceH2. The rate of the spin conversion on the surface drops at temperatures higher than the critical temperature, while that in the gas phase is not sensitive to the temperature in the range of 8 K to 20 K. The impact of the spin conversion on the surfaces becomes more significant with increasing the \ceH2 gas density; the timescale of the spin conversion via the gas phase proton exchange reactions roughly scales with $(n({\rm H_{2}}))^{-0.5}$, while the collisional timescale of \ceH2 to the surface scales with $(n({\rm H_{2}}))^{-1}$ and $\eta$ depends only weakly on $n({\rm H_{2}})$ (see the bottom panel of Fig. 5). The rate of spin conversion in the gas phase and that on the grain surfaces depend differently on physical and chemical conditions ($T,n({\rm H_{2}}),\xi$, $R_{\rm col}(\ce{H2})$, etc.). Therefore, one has to consider the spin conversion both in the gas phase and on grain surfaces for accurate modeling of the ${\rm OPR(\ce{H2})}$ evolution in star- and planet forming regions, which cover wide ranges of the physical and chemical conditions.

Our model was constructed using the experimentally derived binding energy distribution of \ceH2 and the nuclear spin conversion rate on water ice surfaces. In star- and planet-forming regions, gaseous \ceH2 would interact with not only water ice surfaces, but also various types of surfaces, including silicates, graphites, and CO ices. \ceH2 molecules should be formed well before dust grains are coated by water ice mantles. Infrared ice observations have found that the catastrophic CO freeze out happens in dense cores, and ice layers, which mainly consist of CO and \ceCH3OH are formed on top of the water ice layers (e.g., Pontoppidan, 2006). To the best of our knowledge, similar experimental measurements adequate for bare dust grains and ices other than water are not available in the literature. Once such measurements become available, it is straightforward to apply our models to the other types of surfaces, and to simulate the evolution of the OPR(\ceH2) from the formation stage of molecular clouds to the dense core stage (e.g., Furuya et al., 2015), considering the \ceH2 spin conversion both on surfaces and in the gas phase. This is the necessary step for better understanding of the OPR(\ceH2) especially in the early stages of star formation.

5 Conclusion

The ortho-to-para ratio of \ceH2 can significantly affect the molecular evolution, for example deuterium fractionation, in the ISM. The main mechanism of the \ceH2 ortho-para conversion, i.e., whether in the gas phase or on grain surfaces, remains unclear, because the efficiency of the latter in the ISM is not well understood. In this work, we have studied the impact of the nuclear spin conversion of \ceH2 on the water ice surface on the evolution of the overall (gas+ice) \ceH2 under the physical conditions that are relevant to star- and planet-forming regions. We have constructed the rate equation model that considers adsortption of gaseous \ceH2, thermal hopping, desorption, and the nuclear spin conversion of adsorbed \ceH2. We have used the experimentally derived binding energy distribution of \ceH2 and the nuclear spin conversion rate on amorphous water ice surfaces. It was found that the spin conversion efficiency depends on \ceH2 gas density and surface temperature. There are two main factors that control the efficiency of the spin conversion: (i) the efficiency of gaseous and solid \ceH2 interaction and (ii) the probability of the spin-state conversion of an adsorbed \ceH2 molecule before it is desorbed. Enhanced \ceH2 gas density reduces the spin conversion efficiency, because the \ceH2 coverage on the surface increases with increasing the \ceH2 gas density, which hinders gaseous \ceH2 molecules to be adsorbed. The temperature dependence is not monotonic; there is a critical surface temperature at which the efficiency is the maximum. At low temperatures, the exchange of gaseous and icy \ceH2 is inefficient, while at warm temperatures, the residence time of \ceH2 on surfaces is too short for the spin conversion.

By constructing the full gas-ice chemistry model with the \ceH2 spin conversion on grain surfaces, we have found that the spin conversion on the surface dominates over that in the gas-phase at the temperatures below 20 K, regardless of the \ceH2 gas density and the thermal hopping rate of \ceH2 on the surface. We have developed a simple, but accurate formulation to implement the nuclear spin conversion on grain surfaces in existing gas-ice astrochemical models (see Appendix). Our formulation can be applied to any other types of surfaces (e.g., bare dust grain surfaces and CO ice surfaces), once the nuclear spin conversion rate, the sticking probability, and the binding energy distribution of \ceH2 on the other surfaces become available.

We thank the anonymous referee for the useful comments. This work is partly supported by JSPS KAKENHI Grant numbers, 17K14245 and 17H06087.

Appendix A Construction of parameter $\eta$

Here we derive equations for the yield of gaseous ${\rm p\mathchar 45\relax H_{2}}$ per ${\rm o\mathchar 45\relax H_{2}}$ adsorption ($\eta_{\rm op}$) and the yield of gaseous ${\rm o\mathchar 45\relax H_{2}}$ per ${\rm p\mathchar 45\relax H_{2}}$ adsorption ($\eta_{\rm po}$). Our strategy is as follows. We first construct $\eta_{\rm op}$ and $\eta_{\rm po}$ that are adequate in two extreme cases: the fast hopping case, where thermal hopping of adsorbed \ceH2 is fast and the \ceH2 OPR is the same across the surface (Eqs. A3 and A4) and the slow hopping case, where thermal hopping is negligible (Eq. A7). Then we combine the two extremes to obtain a general expression (Eq. A11).

A.1 In the limit of fast hopping

As discussed in Sect. 3, when thermal hopping is considered, $k_{\rm thdes}^{-1}$ is not a good measure for the residence time, but $(k_{\rm des}^{\rm av})^{-1}$ is. We denote the conditional probability that the spin conversion on a surface from ${\rm o\mathchar 45\relax H_{2}}$ to ${\rm p\mathchar 45\relax H_{2}}$ occurs before ${\rm o\mathchar 45\relax H_{2}}$ desorption, and then ${\rm p\mathchar 45\relax H_{2}}$ desorption occurs before reconversion to ${\rm o\mathchar 45\relax H_{2}}$ as $p_{\rm op}$:

[TABLE]

We denote the conditional probability that the spin conversion on a surface from ${\rm o\mathchar 45\relax H_{2}}$ to ${\rm p\mathchar 45\relax H_{2}}$ occurs before ${\rm o\mathchar 45\relax H_{2}}$ desorption, and then reconversion from ${\rm p\mathchar 45\relax H_{2}}$ to ${\rm o\mathchar 45\relax H_{2}}$ occurs before ${\rm p\mathchar 45\relax H_{2}}$ desorption as $r$:

[TABLE]

Then, for example, $rp_{\rm op}$ is the probability that a sequence of spin conversion, ${\rm o\mathchar 45\relax H_{2}}$ $\rightarrow$ ${\rm p\mathchar 45\relax H_{2}}$ $\rightarrow$ ${\rm o\mathchar 45\relax H_{2}}$ $\rightarrow$ ${\rm p\mathchar 45\relax H_{2}}$, occurs on a surface and then ${\rm p\mathchar 45\relax H_{2}}$ desorbs. Using $p_{\rm op}$ and $r$, the yield of gaseous ${\rm p\mathchar 45\relax H_{2}}$ per ${\rm o\mathchar 45\relax H_{2}}$ adsorption in the limit of fast hopping ($\eta^{(f)}$) may be written as follows:

[TABLE]

where $k^{\rm surf}_{\rm conv}=1/\tau^{\rm surf}_{\rm conv}$. Note that $\eta_{\rm op}$ given above considers the possibility of multiple spin conversion on a surface. Similarly, $\eta^{(f)}_{\rm po}$ is given by

[TABLE]

A.2 In the limit of slow hopping

In the limit of slow hopping, we can treat sites with a different energy depth separately. Then we define $\eta^{(s)}$ for each $E_{\rm b}$ using $k_{\rm thdes}$ as follows:

[TABLE]

Again $\eta^{(s)}(E_{\rm b})$ considers multiple spin conversion on a surface, in contrast to Eq. 20, where only single spin conversion is considered.

We define $\langle\eta^{(s)}_{\alpha\beta}\rangle$, where $\alpha$ and $\beta$ are $o$ or $p$, as the average of $\eta^{(s)}_{\alpha\beta}(E_{\rm b})$ weighted by thermal desorption rates;

[TABLE]

where the integration range is from $E_{\rm b}^{\rm thresh}$ (defined below) to $\infty$.

A.3 General case

We denote the threshold binding energy as $E_{\rm b}^{\rm thresh}$; sites with the binding energy lower (higher) than $E_{\rm b}^{\rm thresh}$ is considered in the fast (slow) hopping regime. $E_{\rm b}^{\rm thresh}$ is defined as the binding energy that satisfies

[TABLE]

and $E_{\rm b}^{\rm thresh}>E_{\rm b}^{\rm\,op}$. Note that $k_{\rm hop}(E_{\rm b}^{\rm thresh}\rightarrow E_{\rm b}^{\rm\,op})$ is the lower limit of $k_{\rm hop}(E_{\rm b}\rightarrow E_{\rm b}^{\prime})$ where $E_{\rm b}^{\rm\,op}<E_{\rm b}<E_{\rm b}^{\rm thresh}$ and $E_{\rm b}^{\rm\,op}<E_{\rm b}^{\prime}<E_{\rm b}^{\rm thresh}$ (Eq. 11). For sites with $E_{\rm b}^{\rm\,op}<E_{\rm b}<E_{\rm b}^{\rm thresh}$, the thermal hopping rate is greater than the spin conversion rate, and the thermal desorption rate is smaller than the two rates (i.e., $k_{\rm thdes}<k^{\rm surf}_{\rm conv}<k_{\rm hop}$). Thus the ${\rm OPR(\ce{H_{2}})}$ in such sites are expected to be similar (see the middle panel of Fig. 2). At $n({\rm H_{2}})=10^{4}$ cm*-3*, for example, $E_{\rm b}^{\rm thresh}$ is 596 K for $f=0.5$ and 393 K for $f=0.8$. In sites with $E_{\rm b}<E_{\rm b}^{\rm\,op}$, thermal desorption is more efficient than the spin conversion. Thus the ${\rm OPR(\ce{H_{2}})}$ in sites with $E_{\rm b}<E_{\rm b}^{\rm\,op}$ is mostly determined by the competition between adsorption and thermal hopping from sites with $E_{\rm b}>E_{\rm b}^{\rm\,op}$. For simplicity, we assume that the ${\rm OPR(\ce{H_{2}})}$ in sites with $E_{\rm b}<E_{\rm b}^{\rm\,op}$ is the same as that in sites with $E_{\rm b}^{\rm\,op}<E_{\rm b}<E_{\rm b}^{\rm thresh}$.

Using $E_{\rm b}^{\rm thresh}$, we define $\Theta^{(f)}_{\rm max}$, $\Theta^{(s)}_{\rm max}$, $\Theta^{(f)}$, and $\Theta^{(s)}$ as

[TABLE]

where $\Theta^{(f)}_{\rm max}$ ($\Theta^{(s)}_{\rm max}$) is the fraction of sites that is considered in the fast (slow) hopping regime. $\Theta^{(f)}$ ($\Theta^{(s)}$) is a subset of the \ceH2 coverage that is considered in the fast (slow) hopping regime. Note that $\Theta^{(f)}_{\rm max}+\Theta^{(s)}_{\rm max}=1$ and $\Theta^{(f)}+\Theta^{(s)}=\Theta$.

Using $\eta^{(f)}$ and $\langle\eta^{(s)}\rangle$ defined above, we define the general expression of the spin conversion yield upon adsorption, $\eta$, as

[TABLE]

Note that $\eta_{\rm po}/\eta_{\rm op}=\gamma$, because $\eta^{(f)}_{\rm po}/\eta^{(f)}_{\rm op}=\gamma$ and $\langle\eta^{(s)}_{\rm po}\rangle/\langle\eta^{(s)}_{\rm op}\rangle=\gamma$. We realized that in the evaluation of $k_{\rm des}^{\rm av}$, using $\Theta^{(f)}$ is more reasonable rather than using $\Theta$. Then we redefine $k_{\rm des}^{\rm av}$ as

[TABLE]

and use this in the evaluation of $\eta^{(f)}$. $\eta_{\rm op}$ calculated by Eq. A11 is listed in Table 1.