Quantifying the Chemical Desorption of H2S and PH3 from Amorphous Water-ice Surfaces

Kenji Furuya; Yasuhiro Oba; Takashi Shimonishi

doi:10.3847/1538-4357/ac4260

1. Introduction

Various molecules have been detected in the cold (∼10 K) gas of star-forming regions (see e.g., McGuire 2018). Some molecules, such as CO and N₂H⁺, have gas-phase production pathways that are sufficiently efficient to explain their observed abundances, whereas others, such as CH₃OH and H₂S, do not form by gas-phase reactions efficiently enough to explain the observations (e.g., Geppert et al. 2006; Garrod et al. 2007). Molecules of the latter type are thought to form on dust grains via surface reactions, and to be released into the gas phase via nonthermal desorption processes because thermal desorption is negligible at 10 K. Thus far, several mechanisms of nonthermal desorption have been proposed, including UV photodesorption (e.g., Hama et al. 2009; Bertin et al. 2016; Fuente et al. 2017), stochastic heating and/or sputtering by cosmic rays (e.g., Dartois et al. 2015; Ivlev et al. 2015), and chemical desorption (e.g., Dulieu et al. 2013; Minissale et al. 2016; He et al. 2017; Chuang et al. 2018; Oba et al. 2018; Nguyen et al. 2020). A quantitative understanding of the nonthermal desorption processes is critical for understanding the interplay between gas-phase and ice-surface chemistry.

Among nonthermal desorption processes, chemical desorption is the topic of the present paper. Chemical desorption is caused by the energy released by chemical reactions; after an exothermic surface reaction, (the part of) excess energy goes into the products' translational energy in the direction perpendicular to the surface, leading to desorption (e.g., Fredon et al. 2017). Laboratory experiments have demonstrated that chemical desorption indeed occurs in several reaction systems on astrophysically relevant surfaces (e.g., Dulieu et al. 2013; Minissale et al. 2016; He et al. 2017; Chuang et al. 2018; Oba et al. 2018). In most of these studies, chemical desorption was quantified solely by gas-phase measurements using quadrupole mass spectroscopy (QMS) even though chemical desorption is a surface process. One important result obtained from these experiments is that the fraction of reaction products desorbed from the surface depends on both the reaction and the surface composition. The fraction on amorphous solid water (ASW) is lower than that on graphite or silicate surfaces (Minissale et al. 2016). For some of the reactions studied thus far (e.g., O + O → O₂), the fraction of reaction products released from the surface is reported to be several tens of percent on graphite/silicate surfaces; by contrast, on ASW, it falls below the upper limit of ∼10% measurable in the experiments (see Table 1 in Minissale et al. 2016, for summary). Quantitative understanding of the chemical desorption process, in particular, on ASW is limited because of its low efficiency, although ASW would be more representative of the surface of interstellar grains in star-forming regions than graphite/silicate. Because of our limited understanding of chemical desorption, current astrochemical models often assume that the desorption probability per reactive event is ∼1% (e.g., Garrod 2013; Taquet et al. 2014; Furuya et al. 2015).

For a better understanding of chemical desorption on ASW, researchers have conducted experimental studies in which time-resolved infrared measurements were used to monitor surface species (Oba et al. 2018, 2019; Nguyen et al. 2020, 2021; see also Chuang et al. 2018 for chemical desorption due to reactions involving CO and hydrogen). Oba et al. (2018, 2019) studied chemical desorption upon the reaction of H₂S with H atoms on ASW. They deposited atomic H onto ASW partly covered by H₂S. The following reactions can occur in this system:

$\begin{eqnarray}&&{{\rm{H}}}_{2}{\rm{S}}+{\rm{H}}\to \mathrm{HS}+{{\rm{H}}}_{2},\end{eqnarray} \tag{ 1 }$

$\begin{eqnarray}&&\mathrm{HS}+{\rm{H}}\to {{\rm{H}}}_{2}{\rm{S}}.\end{eqnarray} \tag{ 2 }$

Reactions (1) and (2) are exothermic, with exothermicities of 58 kJ mol⁻¹ (corresponding to ∼7000 K) and 374 kJ mol⁻¹ (∼45,000 K), respectively (Oba et al. 2018), which are much larger than the binding energies of H₂S and HS on ASW (2700 K; Collings et al. 2004; Wakelam et al. 2017). HS is formed from H₂S via a hydrogen abstraction reaction, and hydrogenation of HS leads to the (re)formation of H₂S. H₂S and HS can desorb by chemical desorption in each cycle of this H₂S–HS loop. Because the H₂S–HS loop can continue as long as H atoms are available, a substantial fraction of H₂S and HS can desorb from the ASW surface during an H-atom deposition, even if the probability of chemical desorption per reactive event is small. In Oba et al. (2019), during atomic-H deposition experiments, the authors monitored the abundance of H₂S on ASW by Fourier transform infrared (FTIR) spectroscopy. During H-atom deposition, a decrease in the intensity of the band at 2570 cm⁻¹, which was assigned to the S–H stretching band of H₂S, was observed. Reaction (1) has an activation energy barrier of 1560 K (Lamberts & Kästner 2017), whereas Reaction (2) is barrierless. Most S on ASW might exist as H₂S during the H₂S–HS loop (this is discussed in Section 3). If so, the decrease in the intensity of the band at 2570 cm⁻¹ should be due to chemical desorption of H₂S and/or HS. They found that the amount of H₂S on porous ASW gradually decreased with time (i.e., with increasing amount of H-atom deposition), that ∼60% of the initial amount of H₂S was eventually lost from the ASW surface in experiments at 10 and 20 K, and that ∼30% was lost in an experiment at 30 K. Although the dominant chemical form of S released from the surface (either H₂S or HS) could not be identified, H₂S could be the dominant form, because the exothermicity of Reaction (2) is greater than that of Reaction (1) (Oba et al. 2018).

Nguyen et al. (2020, 2021) experimentally studied the chemical desorption process upon reaction of PH₃ with H atoms on ASW in experiments similar to those of Oba et al. (2018, 2019). They deposited atomic H onto ASW, which was partly covered by PH₃. The following reactions can occur in this system:

$\begin{eqnarray}&&{\mathrm{PH}}_{3}+{\rm{H}}\to {\mathrm{PH}}_{2}+{\rm{H}}2,\end{eqnarray} \tag{ 3 }$

$\begin{eqnarray}&&{\mathrm{PH}}_{2}+{\rm{H}}\to {\mathrm{PH}}_{3}.\end{eqnarray} \tag{ 4 }$

Reactions (3) and (4) are exothermic, with exothermicities of 81.5 kJ mol⁻¹ and 337.5 kJ mol⁻¹, respectively (Molpeceres & Kästner 2021). Similar to the H₂S–HS loop, interconversion between PH₃ and PH₂ occurs through Reactions (3) and (4). In their experiments at 10 K, ∼80% of the initial amount of PH₃ deposited onto porous ASW was released from the surface.

Oba et al. (2019) and Nguyen et al. (2021) estimated a chemical desorption probability per reactive species (i.e., per incident H atoms) of ∼1% on porous ASW for the H₂S + H and PH₃ + H systems. As noted by Oba et al. (2018), the desorption probability per incident H atom corresponds to the lower limit of the desorption probability per reactive event, which astrochemical models require as inputs. This is because H atoms adsorbed onto ASW surfaces can be thermally desorbed or consumed by the H₂ formation via the recombination of two H atoms on the surface. The primary goal of this paper is to constrain the desorption probability of H₂S and PH₃ per reaction event on porous ASW on the basis of the experimental data reported by Oba et al. (2019) and Nguyen et al. (2021), respectively. To this end, we conduct kinetic Monte Carlo (kMC) simulations of their experiments.

The rest of this paper is organized as follows: our numerical model is described in Section 2, and the results are described in Section 3. We compare our derived desorption probability with the theoretical predictions reported in the literature in Section 4. The astrochemical implications of the chemical desorption of H₂S and PH₃ are briefly discussed. We summarize our findings in Section 5.

2. Methods

We intend to constrain the chemical desorption probability per reaction event (P_cd) for H₂S + H and PH₃ + H systems on porous ASW by numerically simulating the experiments of Oba et al. (2019) and Nguyen et al. (2021). For this purpose, we adopt an on-lattice kMC method (Gillespie 1976; Cuppen et al. 2013). In our kMC model, a sequence of processes—adsorption, thermal hopping, and desorption—are modeled as a Markov chain adopting the next reaction method (e.g., Gibson & Bruck 2000; Chang & Herbst 2012). In the case of thermal hopping, a hop in a different direction is treated as a distinct event. This sequence is chosen on the basis of random numbers in combination with the rates for these processes. In contrast to the rate-equation method, which is widely used in the astrochemical community, the position and movement of each chemical species on surfaces are tracked over time in our kMC simulations. Consideration of the binding energy distribution on surfaces in the kMC simulations is then straightforward. Water-ice surfaces are known to contain various adsorption sites with different energy depths (e.g., Amiaud et al. 2006; Hama et al. 2012; Karssemeijer et al. 2014). We will show that the binding energy distribution of atomic H is a key parameter for reproducing the experiments (see Section 3). Details of the numerical aspects of the kMC method can be found in Chang & Herbst (2012), Cuppen et al. (2013), and references therein.

Oba et al. (2019) experimentally investigated the surface reactions of H₂S with H atoms on porous ASW, nonporous ASW, and polycrystalline water ice. They found no strong correlation between the ice structure and the desorption probability per incident H atom in their experiments at 10 K. In the present work, we focus on the experiments on porous ASW because the experiments were conducted at three different surface temperatures: 10, 20, and 30 K. For nonporous ASW and polycrystalline water ice, the experiments were conducted only at 10 K. As shown later, we need experimental data at different surface temperatures to constrain the P_cd and resolve the degeneracy between parameters—in particular, the binding energy distribution for atomic H and the P_cd. Notably, the P_cd is assumed to be independent of the surface temperature throughout this work. This assumption would be reasonable because the exothermicity of Reactions (1)–(4) is much greater than the surface temperature.

For the surface reactions of PH₃ with H atoms, the experiments were conducted on porous ASW only at 10 K (Nguyen et al. 2021). Nonetheless, the P_cd for the PH₃ + H system can be constrained, because we can constrain the binding energy distribution of atomic H by modeling the H₂S + H system. Similar to the H₂S + H system, the desorption probability for PH₃ per incident H atom does not depend on the ice structure (Nguyen et al. 2021).

2.1. Numerical Setup

A square grid is considered with an n × n square lattice providing, at each point, an adsorption site for chemical species on porous ASW. Each site thus has four neighboring sites. Periodic boundary conditions are used. Six different chemical species—atomic H, H₂, H₂S, HS, PH₃, and PH₂—are considered in our models. These different chemical species are assumed to share adsorption sites, whereas the binding energy of each species is assumed to differ even at the same site. Initially, some of the adsorption sites are occupied by H₂S or PH₃, as in the experiments by Oba et al. (2019) and Nguyen et al. (2021). The fluxes of H atoms onto the ASW surface are 5.7 × 10¹³ cm⁻² s⁻¹ and 2.2 × 10¹⁴ cm⁻² s⁻¹ for simulations of the H₂S + H experiments and the PH₃ + H experiments, respectively. These fluxes are consistent with those used in Oba et al. (2019) and Nguyen et al. (2021). We set the flux of H₂ to be two-thirds of that of H atoms, because the H atoms are produced from H₂ with a dissociation probability of >60% (Oba et al. 2018). Refer to Section 2.2.1 for the sticking coefficients adopted in our models.

We consider thermal hopping and thermal desorption of H atoms and H₂ on the ASW surface, whereas those of H₂S, HS, PH₃, and PH₂ are ignored because their binding energies are relatively high (≳2000 K; see Section 2.2.2 for details). If an H atom is adsorbed onto or hops to a site already occupied by another adsorbate, the two can react. For surface reactions, H₂ formation via recombination of two H atoms and Reactions (1)–(4) are considered. Table 1 summarizes the exothermicity, the activation energy barrier (if one exists), and the rate coefficient (s⁻¹) for Reactions (1)–(4). If a H₂ molecule is adsorbed onto or hops to a site already occupied by S-bearing or P-bearing species, the two do not react, and we assume that H₂ is on the topmost layer at the site, with the adsorbate underneath. In this case, we assume that the binding energy of the H₂ is the same as that on the porous ASW. We do not allow the adsorbate underneath a H₂ molecule to experience any chemical processes until the H₂ molecule desorbs or hops to a neighboring site. The presence of H₂ on the surface thus reduces the efficiency of surface reactions in our models. The slowing of surface reactions due to the presence of H₂ in laboratory experiments has been reported by Hama et al. (2015), who studied addition reactions of H atoms with solid benzene. Multiple layers of solid H₂ cannot be created even at 10 K (e.g., Hama et al. 2012). Then we assume that neither H₂ nor H atoms can be adsorbed onto or hop to a site already occupied by either H₂ or atomic H in our simulations. The exception is the adsorption of atomic H onto a site occupied by another atomic H, which instantly leads to the formation of H₂.

Table 1. Summary of Reaction Properties

	((1)) H₂S + H → HS + H₂	((2)) HS + H → H₂S	((3)) PH₃ + H → PH₂ + H₂	((4)) PH₂ + H → PH₃
Reaction energy (kJ mol⁻¹)	−58^a	−374^a	−81.5^c	−337.5^c
Activation energy (K)	1560^b	barrierless	1280^c	barrierless
Rate coefficient (s⁻¹)	8.2 × 10⁶ ^d	⋯	8.7 × 10⁷ ^e	⋯

Notes.

^aGas-phase values (Oba et al. 2018).^bLamberts & Kästner (2017).^cCalculated values using a 20 water cluster model (Molpeceres & Kästner 2021).^dValue at 55 K, which is the lowest temperature studied in quantum chemistry calculations by Lamberts & Kästner (2017). Note that quantum tunneling becomes important at ≲300 K and that the rate constants at 70 and 55 K are almost identical (Lamberts & Kästner 2017).^eValue at 50 K, which is the lowest temperature studied in quantum chemistry calculations by Molpeceres & Kästner (2021).

Download table as: ASCII Typeset image

The thermal desorption rate (s⁻¹) of species i depends on its binding energy on the surface (E_b(i)),

$\begin{eqnarray}&&{k}_{\mathrm{des}}(i)={\nu }_{\mathrm{des}}\exp (-{E}_{{\rm{b}}}(i)/{{kT}}_{s}),\end{eqnarray} \tag{ 5 }$

where ν_des is the characteristic attempt frequency for thermal desorption (1.8 × 10¹² s⁻¹; Amiaud et al. 2015), T_s is the surface temperature, and k is the Boltzmann constant. The binding energies adopted in our models are described in Section 2.2.2. Similarly, the surface diffusion rate for species i via thermal hopping depends on the hopping activation energy (E_hop):

$\begin{eqnarray}&&{k}_{\mathrm{diff},\mathrm{hop}}(i)={\nu }_{\mathrm{hop}}\exp (-{E}_{\mathrm{hop}}/{{kT}}_{s}),\end{eqnarray} \tag{ 6 }$

where ν_hop is the attempt frequency for hopping in one direction, and we arbitrarily assume that ν_des = ν_hop. The activation energy for hopping from a site with a binding energy of E_b to another site with a binding energy of ${E}_{{\rm{b}}}^{\prime}$ is given as (Cazaux et al. 2017; Figure 1):

$\begin{eqnarray}&&{E}_{\mathrm{hop}}({E}_{{\rm{b}}}\to {E}_{{\rm{b}}}^{\prime} )=\alpha \times \min ({E}_{{\rm{b}}},\,{E}_{{\rm{b}}}^{\prime} )+\max (0,\,{E}_{{\rm{b}}}-{E}_{{\rm{b}}}^{\prime} ),\end{eqnarray} \tag{ 7 }$

where α, which is the hopping-to-binding energy ratio, is set to 0.65 for both atomic H and H₂, referring to the recommendation based on calculations of the diffusivity of atomic H on ASW (Asgeirsson et al. 2017). Given the expression for E_hop, the thermal hopping rate exhibits microscopic reversibility, i.e., ${k}_{\mathrm{diff},\mathrm{hop}}({E}_{{\rm{b}}}\to {E}_{{\rm{b}}}^{\prime} )$ / ${k}_{\mathrm{diff},\mathrm{hop}}({E}_{{\rm{b}}}^{\prime} \to {E}_{{\rm{b}}})$ $=\ \exp [-({E}_{{\rm{b}}}-{E}_{{\rm{b}}}^{\prime} )/{{kT}}_{s}]$ (Cuppen et al. 2013). Although Equation (6) is often used in kMC simulations (e.g., Cuppen & Herbst 2007; Garrod 2013), a different expression is proposed for the rate of hopping between adsorption sites with different potential-energy depths (Maté et al. 2020). We tested some models adopting Equation (14) in Maté et al. (2020) instead of Equation (6) and confirmed that the simulation results are not sensitive to the expression of the hopping rate. In our fiducial models, we do not consider the surface diffusion of atomic H by quantum tunneling. This effect will be discussed in Section 4.2.

**Figure 1.** Schematic of the hopping activation barrier from one site to another when the two sites have different binding energies. Based on Figure 11 in Cazaux et al. (2017).
Download figure:
Standard image High-resolution image

For surface reactions, both the Langmuir–Hinshelwood mechanism and the Eley–Rideal mechanism are considered. In our model, when two species arrive at the same adsorption site and can undergo a chemical reaction, we do not follow the details of the competition but determine immediately which process—reaction, hopping, or desorption—occurs by comparing their rates (e.g., Cuppen & Herbst 2007; Chang & Herbst 2012). For barrierless reactions, we assume that the reaction occurs before desorption and hopping, with a probability of unity. For reactions with an activation barrier, the probability of the reaction occurring before one of the species (here, atomic H) leaves to the site, where the two species meet, is given by

$\begin{eqnarray}&&{p}_{{\rm{A}}+{\rm{H}}}=\displaystyle \frac{{k}_{{\rm{A}}+{\rm{H}}}}{{k}_{{\rm{A}}+{\rm{H}}}+{k}_{\mathrm{des}}({\rm{H}})+{{\rm{\Sigma }}}_{i=1}^{4}{k}_{\mathrm{diff},\mathrm{hop}}^{i}({\rm{H}})},\end{eqnarray} \tag{ 8 }$

where A is either H₂S or PH₃, and k_A+H is the reaction rate coefficient. The term ${{\rm{\Sigma }}}_{i=1}^{4}{k}_{\mathrm{diff},\mathrm{hop}}^{i}$ originates from our assumption of a square lattice arrangement of surface sites; each site has four neighboring sites to which atomic H in the site can hop. Equation (8) is implemented in the Monte Carlo algorithm as follows (Chang et al. 2007): A random number between 0 and 1 is generated. If this number is smaller than p_A+H, the reaction occurs. Otherwise, thermal desorption or hopping of the atomic H occurs on the basis of their rates and the random number.

Figure 2 shows ${p}_{{{\rm{H}}}_{2}{\rm{S}}+{\rm{H}}}$ as functions of E_hop for atomic H at 10, 20, and 30 K under the assumption that E_hop is the same for hopping to the four neighboring sites. The value for parameter ${k}_{{{\rm{H}}}_{2}{\rm{S}}+{\rm{H}}}$ was taken from Lamberts & Kästner (2017), who studied Reaction (1) using quantum chemical calculations (Table 1). At 10 K, ${p}_{{{\rm{H}}}_{2}{\rm{S}}+{\rm{H}}}$ is unity (i.e., the reaction is effectively barrierless) when E_hop(H) ≳ 150 K. This threshold E_hop is higher for higher T_s because k_diff,hop for a given E_hop increases with increasing T_s: it is 300 K and 450 K at T_s = 20 K and 30 K, respectively. Notably, as a result of the loop of H₂S–HS interconversion via Reactions (1) and (2), the HS/H₂S abundance ratio approaches ${p}_{{{\rm{H}}}_{2}{\rm{S}}+{\rm{H}}}$ .

After the surface reactions, the products are released from the surface via chemical desorption with a probability of P_cd (see Section 2.2.3). A random number between 0 and 1 is generated; if this number is smaller than the P_cd, the products are instantly released from the surface. Throughout this paper, we use the term chemical desorption rate as the rate of a surface reaction multiplied by the P_cd.

2.2. Chemical Parameters

2.2.1. Sticking Coefficient

The sticking coefficient depends on various parameters, such as the incident kinetic energy of colliding atoms/molecules, surface temperature, and the surface composition. The temperature of atomic H (and H₂) deposited on the ASW was 100 K in the experiments (Oba et al. 2019; Nguyen et al. 2021). In our models, the sticking coefficient for atomic H is taken from Dupuy et al. (2016), who determined the sticking coefficient for atomic H on ASW for various incident kinetic energies and surface temperatures on the basis of classical molecular dynamics (MD) simulations. In our models, the sticking coefficient is set to 0.65 at T_s = 10 K, and 0.6 at 20 K and 30 K. These sticking coefficients are used in both cases where atomic H lands on empty sites (i.e., H₂O) and cases where H lands on either H₂S or PH₃. The sticking coefficient for H₂ is assumed to be the same as that for atomic H.

2.2.2. Binding Energy

Amiaud et al. (2006, 2015) constrained the binding energy distribution of H₂ on porous ASW on the basis of temperature-programmed desorption (TPD) experiments. They found that the binding energy distribution can be described by a polynomial function, adopting ν_des = 1.8 × 10¹² s⁻¹ (Amiaud et al. 2015):

$\begin{eqnarray}{g}_{{{\rm{H}}}_{2}}({E}_{{\rm{b}}})=\left\{\begin{array}{ll}a{\left({E}_{{\rm{b}}}^{\max }-{E}_{{\rm{b}}}\right)}^{1.6} & ({E}_{\ {\rm{b}}}^{\min }\leqslant {E}_{{\rm{b}}}\lt {E}_{{\rm{b}}}^{\max })\\ 0 & ({otherwise})\end{array}\right.\end{eqnarray} \tag{ 9 }$

where a = 4.7 × 10¹¹ sites cm⁻² meV^−2.6 and ${E}_{{\rm{b}}}^{\max }=67.5$ meV (∼780 K) (Amiaud et al. 2015). Notably, the functional form of ${g}_{{{\rm{H}}}_{2}}$ indicates that the number of adsorption sites with smaller E_b is greater. The lower bound of E_b ( ${E}_{\ {\rm{b}}}^{\min }$ ) was not constrained well by the TPD experiments because H₂ was deposited on ASW at a surface temperature of 10 K; even if adsorption sites with very low E_b exist, they would not be probed by the TPD experiments because of the efficient thermal desorption even at 10 K. For example, the timescale of thermal desorption from sites with E_b = 250 K is only ∼0.1 s at 10 K. The minimum value of E_b constrained by their experiments was ∼30 meV (∼350 K) (Amiaud et al. 2006); however, we extrapolate ${g}_{{{\rm{H}}}_{2}}({E}_{{\rm{b}}})$ to lower E_b. In the present work, ${E}_{\ {\rm{b}}}^{\min }$ is treated as a free parameter, and we test three values: ${E}_{\ {\rm{b}}}^{\min }=150$ K, 250 K, and 350 K. According to the numerical investigation of Molpeceres & Kästner (2020), who used MD simulations, approximately 75% of the adsorption sites for H₂ on ASW have an E_b in the range between 150 K and 450 K, and sites with even lower E_b exist, although the peak binding energy distribution is ∼300 K in their simulations. Because the intermolecular potential energy of a H₂–H₂O dimer is ∼120 K (Zhang et al. 1991), the model with ${E}_{\ {\rm{b}}}^{\min }=150$ K corresponds to the situation where the majority of sites have an E_b as low as the intermolecular energy of a dimer. We also tested some models in which the binding energy distributions for H₂ (and atomic H) follow a Gaussian distribution, but found no clear advantage over using a polynomial distribution (see Appendix for details). The cumulative distribution for the H₂ binding energy is shown in Figure 3.

**Figure 3.** Cumulative distribution of the H₂ binding energy adopted in the present work (black line). Note that the minimum value of E_b for H₂ is treated as a free parameter. The red dashed line and the red dotted lines represent the cumulative distribution of the binding energy for atomic H in models with f = 0.8 and f = 0.6, respectively. When f = 1.0, the binding energy distribution for atomic H is identical to that for H₂. See the text for details.
Download figure:
Standard image High-resolution image

The atomic-H binding energy is not well constrained by laboratory experiments because of the high reactivity of atomic H. In this work, we assume that the binding energy of atomic H is given by fE_b(H₂), where f is a scaling factor. Wakelam et al. (2017) have reported that the binding energy of stable molecules on ASW, as determined from the TPD experiments, is proportional to the intermolecular potential energy between a dimer of the molecule and H₂O (see their Figure 1). The intermolecular energy between a dimer of atomic H and H₂O is ∼70 K, whereas that between a dimer of H₂ and H₂O is ∼120 K (Zhang et al. 1991; i.e., the former is ∼0.6 times smaller than the latter). We varied f in the range 0.6 ≤ f ≤ 1.0 in the present work. The cumulative distribution of the atomic-H binding energy is shown by the red lines in Figure 3. Notably, the binding energy distribution for atomic H in our models is controlled by two independent parameters: ${E}_{\ {\rm{b}}}^{\min }({{\rm{H}}}_{2})$ and f. The former parameter controls the extent of E_b(H), whereas the latter parameter controls the maximum value of E_b(H), which is given by ${{fE}}_{{\rm{b}}}^{\max }({{\rm{H}}}_{2})$ . Notably, we assume that the presence of H₂S or PH₃ does not affect the binding energies of atomic H and H₂ on porous ASW.

At the beginning of each kMC simulation, E_b(H₂) at each lattice site is assigned randomly, but it follows ${g}_{{{\rm{H}}}_{2}}$ . The E_b(H) at each lattice site is set to fE_b(H₂); that is, H₂ and atomic H share adsorption sites, and shallow (deep) potential sites for H₂ are also shallow (deep) for atomic H.

According to the TPD experiments, the binding energy of H₂S on water ice is 2700 K (Collings et al. 2004; Garrod & Herbst 2006; Jiménez-Escobar & Muñoz Caro 2011). To the best of our knowledge, no experimental measurements of the binding energy of HS, PH₃, or PH₂ on water ice have been reported. However, theoretically calculated binding energies for HS, PH₃, or PH₂ on water ice are available: 2700 K for HS (Wakelam et al. 2017), 2200 K for PH₃, and 1800 K for PH₂ (mean values; Molpeceres & Kästner 2021). Sil et al. (2021) and Nguyen et al. (2021) also calculated the binding energy of PH₃ on water ice and reported a value similar to that of Molpeceres & Kästner (2021). In our model, we neglect thermal desorption and thermal hopping of H₂S, HS, PH₃, and PH₂. Given the relatively high binding energy of these molecules, this assumption should be valid under the experimental conditions (T_s ≤ 30 K).

The integration of ${g}_{{{\rm{H}}}_{2}}({E}_{{\rm{b}}})$ from ${E}_{{\rm{b}}}^{\max }$ to ${E}_{\ {\rm{b}}}^{\min }$ gives the site density (N_site) on porous ASW. The N_site on porous ASW is known to be greater than that on flat crystalline water ice (10¹⁵ molecules per cm²) (Kimmel et al. 2001; Hidaka et al. 2008). N_site is ∼6 × 10¹⁵, ∼4 × 10¹⁵, and ∼2 × 10¹⁵ sites cm⁻² when ${E}_{\ {\rm{b}}}^{\min }=150$ K, 250 K, and 350 K, respectively. We use a different number of grids to unify the resolutions of our models; n is determined to satisfy the following relation: N_site/(n × n) ≈ 2 × 10¹¹ sites cm⁻². In the models with ${E}_{\ {\rm{b}}}^{\min }=150$ K, 250 K, and 350 K, we use 161 × 161, 128 × 128, and 98 × 98 square lattices, respectively.

In the experiments performed by Oba et al. (2019), 7 × 10¹⁴ molecules cm⁻² of H₂S were deposited on ASW. In this case, the initial surface coverage of H₂S (7 × 10¹⁴/N_site) is ∼12%, ∼18%, and ∼35% in our models with ${E}_{\ {\rm{b}}}^{\min }=150$ K, 250 K, and 350 K, respectively. The initial positions of H₂S molecules on the grid are randomly determined at the beginning of each simulation.

2.2.3. Chemical Desorption Probability

The chemical desorption probabilities per reaction event for Reaction (1) (denoted as P_cd(H₂S)) and for Reaction (3) (P_cd(PH₃)) are treated as free parameters in the present work. In the experiments by Oba et al. (2018, 2019), it was not possible to distinguish between the chemical desorption of H₂S upon Reaction (1) and the chemical desorption of HS upon Reaction (2). If the ratio of the binding energy of a reaction product to the energy released by the reaction is a key parameter for determining P_cd (e.g., Garrod et al. 2007; Minissale et al. 2016), then the P_cd for the hydrogenation reactions would be greater than those for the hydrogen abstraction reactions (see Table 1). In the present work, the P_cd of HS upon Reaction (2) (P_cd(HS)) is always set to be lower than the P_cd(H₂S) by a factor of 10. Notably, our model results do not depend on the ratio of P_cd(HS) to P_cd(H₂S) but do depend on their sum because the chemical desorption of H₂S and HS occurs via the loop of Reactions (1) and (2). Similarly, the P_cd of PH₂ for Reaction (4) (P_cd(PH₂)) is set to be lower than that for Reaction (3) (P_cd(PH₃)) by a factor of 10. Throughout this work, P_cd is assumed to be independent of both T_s and the adsorption sites.

For simplicity, we assume that P_cd of H₂, which can be formed by the recombination of H atoms and the hydrogen abstraction reactions, is unity. This assumption does not affect our results because H₂ can be directly adsorbed on the surface, and because the flux of H₂ is as high as that of H atoms; the maximum formation rate of H₂ on the surface (i.e., one-half of the adsorption rate of atomic H) is comparable to the adsorption rate for H₂. Thus, the main source of H₂ on the surface is direct adsorption.

Our model has the three free parameters (i.e., ${E}_{\ {\rm{b}}}^{\min }({{\rm{H}}}_{2})$ , f, and P_cd); they are summarized in Table 2. In total, we ran 90 different models in which the three free parameters were varied. For each model, we ran the simulations five times with different random seeds and then took an average to reduce fluctuations in the calculation results. We then found the best "by-eye" fit model among the models. This model was considered to be the best fit to the experiments.

Table 2. Free Parameters for the Kinetic Monte Carlo Simulations

Parameter	Values
${E}_{\ {\rm{b}}}^{\min }({{\rm{H}}}_{2})$ ^a	150 K (161²)—250 K (128²)—350 K (98²)
f	0.6-0.7-0.8-0.85-1.0
P_cd	1%–2%–3%–5%–10%–20%

Note.

^aThe values in parentheses are the number of grids used in models with different ${E}_{\ {\rm{b}}}^{\min }({{\rm{H}}}_{2})$ ${E}_{\ {\rm{b}}}^{\min }({{\rm{H}}}_{2})$ values.

Download table as: ASCII Typeset image

3. Results

3.1. H₂S + H

3.1.1. Fiducial Model

We first consider the results from one specific model to understand how the binding energy distribution affects the surface chemistry. The left panel of Figure 4 shows the H₂S abundance on ASW normalized by the initial H₂S abundance as a function of the H-atom exposure time in the model with ${E}_{\ {\rm{b}}}^{\min }({{\rm{H}}}_{2})=150$ K, f = 0.85, and P_cd(H₂S) = 3% (hereafter referred to as the fiducial model). The results show that ∼30% of H₂S is lost within the first 1 minutes in the model at 10 K, which is clearly inconsistent with the experiments. Although the decrease in the H₂S abundance at later times (≳1 minutes) is due to chemical desorption, this rapid loss of H₂S is not due to chemical desorption; as a result of the loop of H₂S–HS interconversion via Reactions (1) and (2), the abundance of H₂S decreases, whereas the abundance of HS increases. At T_s = 10 K, the probability that Reaction (1) occurs per an encounter of atomic H with H₂S is unity (i.e., ${p}_{{{\rm{H}}}_{2}{\rm{S}}+{\rm{H}}}=1$ and the reaction is effectively barrierless) when E_hop(H) ≳ 150 K (see Figure 2). Because both Reactions (1) and (2) are (effectively) barrierless in sites with E_hop(H) ≳ 150 K, the abundances of HS and H₂S become similar as a result of the H₂S–HS interconversion in these sites (Figure 5). This H₂S–HS interconversion should occur before substantial amounts of H₂S and HS are released from the surface by chemical desorption as long as P_cd(H₂S) ≪ 1. Indeed, this rapid destruction and formation of H₂S and HS, respectively, via H₂S–HS interconversion are observed in all our models at 10 K, irrespective of the values of ${E}_{\ {\rm{b}}}^{\min }({{\rm{H}}}_{2})$ , f, and P_cd(H₂S). The degree of H₂S loss depends on the binding energy distribution of atomic H (i.e., the fraction of sites with E_hop(H) ≳ 150 K); however, more than 20% of the initial amount of H₂S is converted to HS within 1 minutes in all of our models at 10 K.

**Figure 4.** H₂S abundance (left panel) and total abundance of H₂S and HS (right panel) on ASW, normalized by the initial H₂S abundance and plotted as functions of the H-atom exposure time in the model with ${E}_{\ {\rm{b}}}^{\min }({{\rm{H}}}_{2})=150$ K, f = 0.85, and P_cd(H₂S) = 3% (lines). Square symbols represent experimental results reported by Oba et al. (2019) and based on FTIR spectroscopy. Black, red, and blue represent surface temperatures of 10 K, 20 K, and 30 K, respectively.
Download figure:
Standard image High-resolution image

**Figure 4.** H₂S abundance (left panel) and total abundance of H₂S and HS (right panel) on ASW, normalized by the initial H₂S abundance and plotted as functions of the H-atom exposure time in the model with ${E}_{\ {\rm{b}}}^{\min }({{\rm{H}}}_{2})=150$ K, f = 0.85, and P_cd(H₂S) = 3% (lines). Square symbols represent experimental results reported by Oba et al. (2019) and based on FTIR spectroscopy. Black, red, and blue represent surface temperatures of 10 K, 20 K, and 30 K, respectively.
Download figure:
Standard image High-resolution image

**Figure 5.** HS/H₂S abundance ratio in adsorption sites with different E_b(H) in the model with ${E}_{\ {\rm{b}}}^{\min }({{\rm{H}}}_{2})=150$ K, f = 0.85, and P_cd(H₂S) = 3% at t = 1 minutes (black bars). For comparison, ${p}_{{{\rm{H}}}_{2}{\rm{S}}+{\rm{H}}}$ is shown as a blue line. In the evaluation of ${p}_{{{\rm{H}}}_{2}{\rm{S}}+{\rm{H}}}$ , E_hop(H) is assumed to be 0.65E_b(H).
Download figure:
Standard image High-resolution image

We conclude that the rapid destruction and formation of H₂S and HS, respectively, via H₂S–HS interconversion are inevitable at 10 K. In principle, if all adsorption sites have E_hop(H) ≪ 150 K (E_b(H) ≪ 230 K as E_hop(H) ∼ 0.65E_b(H)), then the abundance of HS would be negligible, because ${p}_{{{\rm{H}}}_{2}{\rm{S}}+{\rm{H}}}\ll 1$ for all sites and the HS/H₂S abundance ratio would therefore be much lower than unity at 10 K (see Figure 5). In such models, however, reproducing the results of the experiments at 20 and 30 K is difficult, because efficient thermal desorption (e.g., the resident time of atomic H on the surface would be <10⁻⁷ s at 20 K and even shorter at 30 K) prevents the occurrence of Reaction (1), and the chemical desorption of H₂S and HS would be negligible in the experimental timescale. The rapid destruction and formation of H₂S and HS, respectively, observed in the model at 10 K becomes much less obvious in the models at 20 and 30 K (left panel of Figure 4), because the thermal hopping rate of atomic H increases exponentially with increasing T_s. Thus, ${p}_{{{\rm{H}}}_{2}{\rm{S}}+{\rm{H}}}$ is much less than unity for the majority of the adsorption sites; i.e., the HS/H₂S abundance ratio $\approx {p}_{{{\rm{H}}}_{2}{\rm{S}}+{\rm{H}}}\ll 1$ .

Oba et al. (2019) confirmed the chemical desorption of H₂S (and HS) by both FTIR spectroscopy in atomic-H deposition experiments and QMS measurements in the TPD experiments after H-atom deposition. During H-atom deposition, a decrease was observed in the magnitude of the band at 2570 cm⁻¹, which was assigned to an S–H stretching band of H₂S. However, distinguishing between H₂S and HS by FTIR is difficult because of overlap of the absorption bands (Jiménez-Escobar & Muñoz Caro 2011). In the TPD experiments after H-atom deposition, a decrease in the H₂S amount compared with the initial H₂S amount deposited on the ASW surface was found by QMS. The QMS measurements confirmed that no S-bearing species except H₂S were desorbed; however, both H₂S⁺ and HS⁺ were produced as fragment ions of H₂S, which could mask the existence of HS. Thus, we cannot rule out the possibility that some HS along with H₂S existed on the ASW in the experiments.

In the present work, as a working hypothesis, we interpret the 2570 cm⁻¹ band as the sum of the S–H stretching modes of H₂S and HS rather than solely as the stretching mode of H₂S. We also assume that the band intensity for H₂S and HS on ASW are the same because no literature value has been reported for HS. Under these assumptions, in the rest of this paper, we focus on the total abundance of S (i.e., H₂S + HS) on the ASW surface rather than the H₂S abundance in our simulations. In the right panel of Figure 4, the vertical axis represents the total abundance of H₂S and HS normalized by the initial abundance of H₂S at a given time. In this case, we do not observe a rapid decrease at the early times (≲1 minute) at 10 K, whereas the results at 20 and 30 K are similar between the two panels.

In the experiments, the abundance of H₂S (and HS) on ASW decreased with time and eventually became almost constant in the experimental timescale (Figure 4). This nonlinear behavior reflects the binding energy distributions of atomic H and H₂ on ASW. Figure 6 shows the number of adsorption sites occupied by either H₂S or HS at t = 0 minutes (gray), 10 minutes (red), and 150 minutes (blue) in the fiducial model at T_s = 10 K (left), 20 K (middle), and 30 K (right) as functions of E_b(H). Because the surface coverage of H₂S is ∼12% at t = 0 minutes (Section 2.2.2), the majority of adsorption sites are empty, and these empty sites are not included in the figure. At 10 K and at t = 150 minutes, H₂S and HS are populating adsorption sites with binding energies E_b(H) following a bimodal distribution (E_b(H) ≲ 150 K and E_b(H) ≳ 400 K). In the deep sites (E_b(H) ≳ 450 K, which corresponds to E_b(H₂) ≳ 530 K as E_b(H) = 0.85E_b(H₂) in the fiducial model), H₂S and HS are buried by H₂; because H₂ does not easily hop to neighboring sites, further reactions of H₂S and HS with atomic H are hindered, as described in Section 2.1. The hopping timescale for H₂ from a site with E_b(H₂) = 530 K to a site with E_b(H₂) ≲ 450 K exceeds the duration of the experiments (150 minutes). In the shallow sites (E_b(H) ≲ 150 K), ${p}_{{{\rm{H}}}_{2}{\rm{S}}+{\rm{H}}}$ is less than unity, which slows the loop of the H₂S–HS interconversion; the rate of the chemical desorption is thereby reduced. H₂S only in the sites with E_b(H) ≲ 250 K, and E_b(H) ≲ 350 K remains on the surface after 150 minutes at 20 K and 30 K, respectively. H₂ readily hops from the deepest sites to shallower sites and is thermally desorbed at 20 K and 30 K; essentially no H₂ remains on ASW. A bimodal distribution is therefore not observed at 20 and 30 K.

**Figure 6.** Number of adsorption sites occupied by either H₂S or HS at t = 0 minutes (gray), at t = 10 minutes (red), and at t = 150 minutes (blue) in the models with f = 0.85, ${E}_{\ {\rm{b}}}^{\min }=150$ K, and P_cd(H₂S) = 3% as functions of the binding energy of atomic H in the adsorption sites. The left, middle, and right panels show the model for a surface temperature of 10 K, 20 K, and 30 K, respectively. Note that the majority of adsorption sites are empty and that these empty sites are not included in the figure.
Download figure:
Standard image High-resolution image

**Figure 6.** Number of adsorption sites occupied by either H₂S or HS at t = 0 minutes (gray), at t = 10 minutes (red), and at t = 150 minutes (blue) in the models with f = 0.85, ${E}_{\ {\rm{b}}}^{\min }=150$ K, and P_cd(H₂S) = 3% as functions of the binding energy of atomic H in the adsorption sites. The left, middle, and right panels show the model for a surface temperature of 10 K, 20 K, and 30 K, respectively. Note that the majority of adsorption sites are empty and that these empty sites are not included in the figure.
Download figure:
Standard image High-resolution image

We can roughly estimate the threshold value of E_b(H) ( ${E}_{\ {\rm{b}}}^{\mathrm{thresh}}({\rm{H}})$ ); most H₂S in sites with a binding energy lower than ${E}_{\ {\rm{b}}}^{\mathrm{thresh}}({\rm{H}})$ remains on ASW after 150 minutes of H-atom deposition, avoiding chemical desorption. The flux of H atoms is 5.7 × 10¹³ cm⁻² s⁻¹ with a sticking probability of ∼0.6; thus, the fluence (i.e., the time integral of the flux) of H atoms in 150 minutes is ∼3 × 10¹⁷ molecules cm⁻². Because 7 × 10¹⁴ molecules cm⁻² of H₂S are initially present on ASW, the average number of H atoms available to desorb one H₂S molecule via chemical desorption is ∼400. Let us assume that adsorbed H atoms are consumed either by the loop of Reactions (1)–(2) in sites with E_b or by thermal desorption. Under this assumption, the probability that Reaction (1) occurs before the thermal desorption of a H atom is given by ${k}_{{{\rm{H}}}_{2}{\rm{S}}+{\rm{H}}}/({k}_{{{\rm{H}}}_{2}{\rm{S}}+{\rm{H}}}+{k}_{\mathrm{des}}({\rm{H}}))$ . Then, by solving $400[{k}_{{{\rm{H}}}_{2}{\rm{S}}+{\rm{H}}}/({k}_{{{\rm{H}}}_{2}{\rm{S}}+{\rm{H}}}+{k}_{\mathrm{des}}({\rm{H}}))]{P}_{\mathrm{cd}}({{\rm{H}}}_{2}{\rm{S}})=1$ , we obtain ${E}_{\ {\rm{b}}}^{\mathrm{thresh}}({\rm{H}})$ to be ∼100 K, 200 K, 300 K at T_s = 10 K, 20 K, and 30 K, respectively, when P_cd(H₂S) = 3%. These values are consistent with the numerical results corresponding to T_s = 20 K and 30 K (see Figure 6). In the numerical simulation at 10 K, H₂S and HS in adsorption sites even with ${E}_{{\rm{b}}}({\rm{H}})\gt {E}_{\ {\rm{b}}}^{\mathrm{thresh}}({\rm{H}})$ remain on the surface. The lack of agreement between the previously discussed simple argument and the numerical results at 10 K indicates the importance of trapping H atoms in deep-potential sites, followed by H₂ formation via recombination with another H atom. This effect can reduce the number of H atoms available for the H₂S–HS interconversion loop. Even for E_b(H) = 660 K (i.e., the deepest site), the timescale of hopping to a neighboring site is around 0.001 s at 20 K, which is comparable to the average adsorption timescale of one atomic H in our model, and the hopping timescale is even shorter at 30 K. Then the trapping of atomic H in deep sites is only important at 10 K.

Taken together, the nonlinear behavior of the total abundance of H₂S and HS on ASW can be understood as follows: At early times (≲10 minutes), chemical desorption occurs at adsorption sites, where the loop of H₂S–HS interconversion can occur most efficiently (i.e., ${p}_{{{\rm{H}}}_{2}{\rm{S}}+{\rm{H}}}\sim 1$ ). At later times, the rate of chemical desorption is reduced, because the remaining H₂S is in sites less favorable for the H₂S–HS interconversion loop; that is, either ${p}_{{{\rm{H}}}_{2}{\rm{S}}+{\rm{H}}}\lt 1$ or H₂S is buried by H₂, hindering further reactions. Eventually, H₂S only in sites with E_b(H) lower than ${E}_{\ {\rm{b}}}^{\mathrm{thresh}}({\rm{H}})$ remains on the ASW. At 10 K, H₂S (and HS) also remain in deep-potential sites, where the presence of H₂ hinders surface reactions.

The previous discussion reveals that the binding energy distribution for atomic H can be directly constrained on the basis of the 20 K and 30 K experiments, because the fraction of H₂S remaining on the ASW corresponds to the fraction of adsorption sites with E_b(H) lower than ${E}_{\ {\rm{b}}}^{\mathrm{thresh}}({\rm{H}})$ . As previously mentioned, ${E}_{\ {\rm{b}}}^{\mathrm{thresh}}({\rm{H}})$ is ∼200 K and 300 K at T_s = 20 K and 30 K, respectively. ${E}_{\ {\rm{b}}}^{\mathrm{thresh}}({\rm{H}})$ depends logarithmically on P_cd(H₂S); i.e., the dependence of ${E}_{\ {\rm{b}}}^{\mathrm{thresh}}({\rm{H}})$ on P_cd(H₂S) is weak. The abundance of H₂S becomes ∼40% and ∼70% of the initial H₂S abundance at 150 minutes in the 20 K and 30 K experiments, respectively. Thus, ∼70% of the adsorption sites on the porous ASW should have E_b(H) ≲ 300 K, whereas ∼40% of the adsorption sites should have E_b(H) ≲ 200 K. Hama et al. (2012) found in their experiments that the adsorption sites for atomic H on the ASW surface can be categorized into three groups: shallow (E_hop(H) < 18 meV ∼ 210 K), middle (E_hop(H) = 22 meV ∼ 260 K), and deep (E_hop(H) > 30 meV ∼ 350 K) potential sites. Although the fraction of each of the three sites was not constrained well in their work, the shallow sites were found to be dominant on the surface. If we assume a hopping-to-binding energy ratio of 0.65 (Asgeirsson et al. 2017), the binding energy distribution constrained as previously described indicates that ∼70% of the adsorption sites have E_hop(H) ≲ 0.65 × 300 ∼ 200 K, consistent with the diffusion barriers constrained by Hama et al. (2012).

3.1.2. Lower Limit of H₂S Desorption Probability

As previously discussed, the initial decrease of the H₂S and HS abundances at ≲10 minutes in the experiments is likely due to chemical desorption from the adsorption sites, where ${p}_{{{\rm{H}}}_{2}{\rm{S}}+{\rm{H}}}=1$ , and thus, the loop of H₂S–HS interconversion can occur most efficiently. We expect that reproducing the experimental results at ≲10 minutes using models would be easier than reproducing the experimental results at >10 minutes, where chemical desorption occurs even in sites with ${p}_{{{\rm{H}}}_{2}{\rm{S}}+{\rm{H}}}\lt 1$ . Thus, as a first step, attempting to constrain the P_cd(H₂S) while focusing only on the experimental results at ≲ 10 minutes would be worthwhile. For this purpose, we additionally ran models with a single value of E_hop(H) for all the adsorption sites. We varied E_hop(H) in the range where ${p}_{{{\rm{H}}}_{2}{\rm{S}}+{\rm{H}}}=1$ (i.e., E_hop(H) > 150 K at T_s = 10 K). In these models, we neglected thermal desorption of atomic H for simplicity, whereas we considered H₂ formation by recombination of two H atoms. The other parameters were the same as those described in Section 2.2. We here focus on the 10 K experiments because thermal desorption of atomic H should occur efficiently at 20 and 30 K, invalidating the aforementioned assumption.

Figure 7 shows the total abundance of H₂S and HS normalized by the initial H₂S abundance as functions of the H-atom exposure time at 10 K, where E_hop(H) and P_cd(H₂S) were varied. When E_hop(H) = 400 K or larger, E_hop(H) and P_cd(H₂S) degenerate (i.e., we cannot constrain P_cd(H₂S)). However, when E_hop(H) = 200 and 300 K, the results do not depend on E_hop(H) but do depend on P_cd(H₂S) because the diffusion of H atoms by thermal hopping is fast, and almost all the adsorbed H atoms are consumed by either Reaction (1) or Reaction (2) before another H atom is adsorbed onto the ASW surface. That is, Reactions (1) and (2) are adsorption-limited, and the number of H₂S and HS released from the surface is simply given by P_cd(H₂S) multiplied by the number of adsorbed H atoms. In the experiments, some adsorbed H atoms might have been trapped in deep-potential sites and/or desorbed thermally; these effects are not included in the models with a single value of E_hop but indeed occur in the models where the binding energy distribution is taken into account. Then, what we can constrain here is the lower limit of P_cd(H₂S) by comparing the models with E_hop ≤ 300 K with the experiments at t ≲ 10 minutes. If we consider the experimental results only at t ≤ 5 minutes, the lower limit of P_cd(H₂S) is ∼3%. If we consider the experimental results at t ≤ 10 minutes, the lower limit of P_cd(H₂S) is ∼2%. As a conservative choice, we consider the lower limit of P_cd(H₂S) as 2%.

3.1.3. Constraining the H₂S Desorption Probability

Figure 8 shows the total amount of H₂S and HS divided by the initial amount of H₂S at 10, 20, and 30 K as functions of the H-atom exposure time in the models with P_cd(H₂S) = 2%, where the other free parameters, f and ${E}_{\ {\rm{b}}}^{\min }({{\rm{H}}}_{2})$ , were varied. The model results are sensitive to f and ${E}_{\ {\rm{b}}}^{\min }({{\rm{H}}}_{2})$ at 20 and 30 K, whereas the impact of the two parameters is less significant at 10 K. At 20 and 30 K, H₂S and HS in adsorption sites with higher ${E}_{\ {\rm{b}}}^{\min }({\rm{H}})$ are preferentially lost from the surface (Figure 6). Lower values of f and ${E}_{\ {\rm{b}}}^{\min }({{\rm{H}}}_{2})$ lead to a decrease in the fraction of deeper sites. As a result, the chemical desorption rate is reduced with decreasing f and/or ${E}_{\ {\rm{b}}}^{\min }({{\rm{H}}}_{2})$ , as evident in Figure 8. At 10 K, H₂S and HS become populating adsorption sites with binding energies ${E}_{\ {\rm{b}}}^{\min }({\rm{H}})$ following a bimodal distribution (Figure 6). Lower values of f and ${E}_{\ {\rm{b}}}^{\min }({{\rm{H}}}_{2})$ lead to a decrease in the fraction of deep sites (i.e., the chemical desorption rate tends to increase), whereas they lead to an increase in the fraction of shallow sites (i.e., the chemical desorption rate tends to decrease). As a result, the 10 K model is less sensitive to f and ${E}_{\ {\rm{b}}}^{\min }({{\rm{H}}}_{2})$ than the 20 and 30 K models.

Because 2% is the lower limit of P_cd(H₂S), the models that overestimate the amount of H₂S and HS desorbed from the ASW surface compared to the experiments are ruled out. Examining the 20 and 30 K models, we can rule out the majority (6 of 9) of the models shown here. The three models not yet ruled out are those with ${E}_{\ {\rm{b}}}^{\min }({{\rm{H}}}_{2})=250$ K and f = 0.6, with ${E}_{\ {\rm{b}}}^{\min }({{\rm{H}}}_{2})=150$ K and f = 0.6, and with ${E}_{\ {\rm{b}}}^{\min }({{\rm{H}}}_{2})=150$ K and f = 0.8 (shown by thick lines in Figure 8). Notably, the three models share the common feature of a low-atomic-H binding energy; more than one-half of the adsorption sites have a low binding energy for atomic H (≲300 K; see Figure 3) and, thus, a low hopping barrier (≲200 K).

Among the three models, we disfavor the model with ${E}_{\ {\rm{b}}}^{\min }({{\rm{H}}}_{2})=250$ K and f = 0.6 because, although the model moderately well reproduces the experimental results at T_s = 10 K, too much H₂S and HS remain on the surface at 30 K compared to the experiments; the experiments at 10, 20, and 30 K are difficult to fit simultaneously by varying P_cd(H₂S). We thus have two remaining possibilities: (i) the model with f = 0.6 and ${E}_{\ {\rm{b}}}^{\min }({{\rm{H}}}_{2})=150$ K is more favorable, and P_cd(H₂S) is much higher than the lower limit value of 2%, because too much H₂S and HS remain on the surface compared with the experiments at all the investigated temperatures; or (ii) the model with f = 0.8 and ${E}_{\ {\rm{b}}}^{\min }({{\rm{H}}}_{2})=150$ K is more favorable, and P_cd(H₂S) is close to and slightly greater than 2%, because the model somehow underestimates chemical desorption at all the investigated temperatures.

To explore the first possibility, we tested additional models, varying f from 0.6 to 0.8 and varying P_cd(H₂S) from 5% to 20% but fixing ${E}_{\ {\rm{b}}}^{\min }({{\rm{H}}}_{2})=150$ K. The results for these models are shown in Figure 13 in the appendix. We find that reproducing the 10–30 K experiments simultaneously using these models with a relatively high P_cd(H₂S) of >5% is difficult; the models at 10 K overestimate the amount of H₂S and HS desorbed from the surface compared with the experiments, or the models at 30 K underestimate the amount of H₂S and HS desorbed from the surface compared with the experiments.

Therefore, we conclude that the second possibility is more favorable; i.e., the P_cd(H₂S) is close to the lower limit value of 2%. After some exploration, we found that the model with f = 0.85, ${E}_{\ {\rm{b}}}^{\min }({{\rm{H}}}_{2})=150$ K, and P_cd(H₂S) = 3% (i.e., the fiducial model) reasonably well reproduces the 10, 20, and 30 K experimental results simultaneously (see the right panel of Figure 4). We confirmed that the reduced χ² for this model (=1.85) is the minimum among all the models applied in the present work. We consider this model as the best-fit model.

As noted in Section 2.2.3, our model results depend on the sum of P_cd(H₂S) and P_cd(HS), but do not depend on the ratio of P_cd(HS) to P_cd(H₂S), which is assumed to be 0.1 throughout the present work. Thus, strictly speaking, the constraint applied here is that the sum of P_cd(H₂S) and P_cd(HS) is 3.3%. Nevertheless, P_cd(H₂S) can be reasonably expressed as ∼3% because P_cd(H₂S) would be greater than P_cd(HS) according to the theoretical models for chemical desorption (see Section 4.3). Finally, we also ran models in which the binding energy distributions of atomic H and H₂ are given by a Gaussian distribution rather than by a polynomial distribution; we found no clear advantage to using a Gaussian distribution (see the appendix).

Our best-fit model underestimates the total abundance of H₂S and HS remaining on the ASW surface at 10 K compared with the experiments, although it better reproduces the 20 K and 30 K experimental results. This underestimation might indicate that the presence of H₂ on the surface lowers the P_cd. The ASW surface is partly covered by H₂ at 10 K, whereas essentially no H₂ remains on the surface at 20 and 30 K. Minissale & Dulieu (2014) found in their experiments that the chemical desorption probability of O₂ upon the reaction between two O atoms decreases with increasing coverage of O₂ on oxidized graphite. The presence of an adsorbed species might enhance the dissipation of the excess energy produced by chemical reactions, resulting in a lowering of the P_cd (Minissale & Dulieu 2014). Such an effect of an adsorbed species (H₂ in our case) on the P_cd is not considered in our models.

3.2. Constraining the PH₃ Desorption Probability

We here constrain the P_cd for PH₃ upon Reaction (4). As previously described, the experiments of H₂S + H at 10, 20, and 30 K are reasonably well reproduced by the model with f = 0.85 and ${E}_{\ {\rm{b}}}^{\min }({{\rm{H}}}_{2})=150$ K. Therefore, in this subsection, we fixed the parameters f and ${E}_{\ {\rm{b}}}^{\min }$ to these values. As in Section 3.1, we discuss the total amount of PH₃ and PH₂ on the surface divided by the initial amount of PH₃, again because the rapid conversion of PH₃ to PH₂ is inevitable at T_s = 10 K in our models; the rate of Reaction (3) is even greater than that of Reaction (1) (Table 1). Figure 9 compares the experimental results with the model, where the P_cd(PH₃) was varied. Notably, the H-atom flux was approximately four times greater in the experiments involving PH₃ + H (Nguyen et al. 2021) than in those involving H₂S + H (Oba et al. 2019). We find that the model with P_cd(PH₃) = 4% reasonably well reproduced the experiments at 10 K.

**Figure 9.** Relative abundances of phosphorous (PH₃ + PH₂) on porous ASW at 10 K as functions of the H-atom exposure time in the kMC simulations (lines) and in the experiments by Nguyen et al. (2021; black squares). Solid, dashed, and dotted lines represent the models with P_cd(PH₃) of 4%, 2%, and 10%, respectively. In all the simulations, ${E}_{\ {\rm{b}}}^{\min }({{\rm{H}}}_{2})$ and f are set to 150 K and 0.85, respectively.
Download figure:
Standard image High-resolution image

4. Discussion

4.1. Uncertainties in the Flux of Atomic H

Thus far, the flux of H atoms was set to 5.7 × 10¹³ cm⁻² s⁻¹ and 2.2 × 10¹⁴ cm⁻² s⁻¹ in the H₂S + H model and in the PH₃ + H model, respectively. These fluxes were estimated by Oba et al. (2019) and Nguyen et al. (2021). The authors did not discuss the uncertainties in the estimation of the atomic-H flux in their experiments; however, empirically, the absolute uncertainties can be as large as 50% (private communication). Here, we explore how the flux uncertainty affects the constraints on P_cd.

We tested two additional models for the H₂S + H system, varying the H-atom flux. In one model, the H-atom flux was set to 3.8 × 10¹³ cm⁻² s⁻¹ (50% lower than the fiducial value); in the other model, it was 8.6 × 10¹³ cm⁻² s⁻¹ (50% higher). The other parameters were the same as in our fiducial model: ${E}_{\ {\rm{b}}}^{\min }({{\rm{H}}}_{2})=150$ K, f = 0.85, and P_cd(H₂S) = 3%. The model with the 50% lower H flux and P_cd(H₂S) = 3% and the model with the fiducial H flux and P_cd(H₂S) = 1.5% gave almost identical results, whereas the model with the 50% higher H flux and P_cd(H₂S) = 3% and the model with the fiducial H flux and P_cd(H₂S) = 4.5% gave almost identical results. Thus, the uncertainty in the H-atom flux is inversely linearly transferred to the uncertainty in the P_cd(H₂S), likely because the system is adsorption-limited rather than diffusion-limited. Considering the empirical uncertainty in the H-atom flux (50%), we conclude that the P_cd(H₂S) is 3% ± 1.5%, whereas the P_cd(PH₃) is 4% ± 2%.

4.2. Diffusion of Atomic H by Quantum Tunneling

In our fiducial models, thermal hopping is treated solely as a mechanism of surface diffusion. Asgeirsson et al. (2017) theoretically studied the diffusion of atomic H on amorphous ice, considering both thermal hopping and quantum tunneling, and showed that tunneling is important only for T_s ≲ 10 K. Kuwahata et al. (2015) experimentally showed that tunneling dominates thermal hopping for atomic H on crystalline water ice at 10 K. However, tunneling is less important for the diffusivity of atomic H on ASW, where greater heterogeneity (i.e., a greater range in E_b(H)) exists compared to the heterogeneity of crystalline water-ice surfaces. The rationale is that, even though H atoms can tunnel between shallow sites, the rate-limiting step for long-range diffusion is the transition out of deep sites, which still requires thermal activation (Smoluchowski 1983; Kuwahata et al. 2015; Asgeirsson et al. 2017). These experimental and theoretical studies indicate that tunneling has limited importance for long-range diffusion of atomic H on ASW at ≥10 K. On the other hand, tunneling diffusion might play a role in determining the reaction probability because it is determined by the competition between reaction and short-range diffusion (Equation (8)).

To check the effect of tunneling diffusion on our results, we tested additional models in which the tunneling diffusion of atomic H was considered. Lamberts et al. (2014) proposed tunneling rates for exothermic reactions based on arguments of microscopic reversibility. Because the diffusion rates should also obey microscopic reversibility, referring to Lamberts et al. (2014), we calculate the tunneling diffusion rates for atomic H between two sites with different binding energies ( ${E}_{{\rm{b}}}\lt {E}_{{\rm{b}}}^{\prime}$ ) as

$\begin{eqnarray}&&{k}_{\mathrm{diff},\mathrm{qt}}({E}_{{\rm{b}}}\to {E}_{{\rm{b}}}^{\prime} )=\nu \exp \left(-\displaystyle \frac{2d}{{\hslash }}\sqrt{2{m}_{{\rm{H}}}\alpha {E}_{{\rm{b}}}}\right),\end{eqnarray} \tag{ 10 }$

$\begin{eqnarray}&&\begin{array}{l}{k}_{\mathrm{diff},\mathrm{qt}}({E}_{{\rm{b}}}^{\prime} \to {E}_{{\rm{b}}})\\ =\ {k}_{\mathrm{diff},\mathrm{qt}}({E}_{{\rm{b}}}\to {E}_{{\rm{b}}}^{\prime} )\exp \left(-\displaystyle \frac{{E}_{{\rm{b}}}^{\prime} -{E}_{{\rm{b}}}}{{{kT}}_{s}}\right),\end{array}\end{eqnarray} \tag{ 11 }$

where surface sites are assumed to be separated by a rectangular barrier of thickness d, and m_H is the mass of atomic H. Because the barrier thickness is unknown, we here consider two values: d = 1 and 2 Å. Notably, when atomic H moves from a deeper site to a shallower site, the part of activation energy barrier that corresponds to the difference in the binding energies of the two sites should be overcome thermally (see Figure 1). This consideration is accounted for by the factor $\exp (-({E}_{{\rm{b}}}^{\prime} -{E}_{{\rm{b}}})/{{kT}}_{s})$ in Equation (11). We compared the tunneling rate and thermal hopping rate and chose the higher of the two as the diffusion rate.

Figure 10 shows the effect of the tunneling diffusion at T_s = 10 K when d = 1 Å (red solid line) and d = 2 Å (red dashed line). Compared with our best-fit model constrained as described in Section 3 (black solid line), tunneling increases the fraction of H₂S and HS released to the gas phase by chemical desorption to some extent. As a result, the model with the tunneling diffusion underestimates the amount of H₂S and HS remaining on the ASW surface after 25 minutes compared with the experimentally observed amount. Even if we adopt P_cd(H₂S) = 2%, which is the lower-limit value, the model with the tunneling underestimates the amount of H₂S and HS remaining on the ASW surface (blue lines). As discussed in Section 3.1.1, some H₂S and HS in sites with ${E}_{{\rm{b}}}({\rm{H}})\gt {E}_{\ {\rm{b}}}^{\mathrm{thresh}}({\rm{H}})\sim 100$ K remain on the surface at 10 K in our fiducial model (see Figure 6) because of the trapping of H atoms in deep-potential sites, followed by the formation of H₂ via recombination with another H atom. Because of the tunneling, the effect of trapping is reduced, and almost all of the H₂S and HS in sites with E_b(H) ∼ 150 K are released to the gas phase by chemical desorption in the model with the tunneling diffusion. A similar effect is observed in the PH₃ + H model at 10 K. At 20 and 30 K, the results obtained using the models with and without the tunneling diffusion are almost identical (not shown).

When the tunneling diffusion is considered, the model with P_cd(H₂S) = 2% better reproduces the H₂S + H experimental results at T_s = 10, 20, and 30 K rather than the model with P_cd(H₂S) = 3%, assuming ${E}_{\ {\rm{b}}}^{\min }({{\rm{H}}}_{2})=150$ K and f = 0.85. Note, however, that the reduced χ² for this model is larger than that for our best-fit model constrained in Section 3 (∼4 versus ∼2). For the PH₃ + H system, the model with P_cd(H₂S) = 3% better reproduces experimental results rather than the model with P_cd(PH₃) = 4%, when the tunneling diffusion is considered. The reduced χ² for this model is similar to that for our best-fit model constrained in Section 3. Therefore, slightly lower P_cd than that constrained in Section 3 is favorable, when the tunneling diffusion is considered.

4.3. Comparison with Theoretical Models of Chemical Desorption

Chemical desorption is caused by the energy released by reactions; after an exothermic surface reaction, some of the excess energy is dissipated into the products' translational energy in the direction perpendicular to the surface, leading to desorption (e.g., Fredon et al. 2017). The efficiency of chemical desorption depends on (i) the excess energy released by reaction (E_react), (ii) the binding energy of products to the surface, (iii) the fraction of excess energy that remains in the products and is not lost to the solid surface, and (iv) how the energy of the products is distributed among all degrees of freedom. The last two parameters are currently uncertain, which limits our quantitative understanding of chemical desorption. Several authors have suggested general formalisms of chemical desorption to be used in astrochemical models (Garrod et al. 2007; Minissale et al. 2016; Fredon et al. 2021). These formalisms include free parameters that should be determined by laboratory experiments or/and quantum chemistry calculations. Here, we compare the P_cd constrained in the present work with that predicted by the theoretical models of chemical desorption in the literature.

Garrod et al. (2007) proposed the following formalism in the case of one-product reactions, applying Rice–Ramsperger–Kessel–Marcus theory:

$\begin{eqnarray}&&{P}_{\mathrm{cd},{\rm{G}}07}=\displaystyle \frac{a\gamma }{1+a\gamma },\end{eqnarray} \tag{ 12 }$

$\begin{eqnarray}&&\gamma ={\left(1-\displaystyle \frac{{E}_{{\rm{b}}}}{{E}_{\mathrm{react}}}\right)}^{s-1},\end{eqnarray} \tag{ 13 }$

where γ is the probability for an energy E > E_b to be present in the admolecule–surface bond, and s is the number of vibrational degrees of freedom, including binding to the surface. The parameter a = ν/ν_s is the ratio of the desorption attempt frequency to the frequency at which the reaction energy is lost to the surface. The value of a is unknown, and a = 0.01 (i.e., P_cd,G07 ∼ 1%) is often assumed in astrochemical models (e.g., Garrod et al. 2007; Taquet et al. 2014; Furuya et al. 2015). In the case of two-product reactions, chemical desorption is assumed to not occur (Garrod et al. 2007).

Minissale et al. (2016) proposed a different formulation from that of Garrod et al. (2007) to explain the results of their chemical desorption experiments. They treated the energy dissipation as an elastic collision, where some of the kinetic energy is transferred to the surface from the product species, and assumed that excess energy that remains in the product species is spread equally over all degrees of freedom of the product species. Their formulation with an effective surface mass of 130 amu reproduced the efficiency of chemical desorption for the studied reactions on graphite surfaces (i.e., rigid surfaces) within the margin of error. However, the use of the formulation by Minissale et al. (2016) is questionable for ASW surfaces, where intermolecular interactions dominate, and the collisions between admolecules and the surface are inelastic (Fredon et al. 2017). Because our work is focused on the chemical desorption on ASW, hereafter we do not discuss the formalism of Minissale et al. (2016).

More recently, Fredon et al. (2021) proposed the following formulation based on classical MD simulations, where excess energy is transferred to a molecule adsorbed on ASW, and the fate of the molecule is tracked (see also Fredon et al. 2017):

$\begin{eqnarray}&&{P}_{\mathrm{cd},{\rm{F}}21}=\alpha \left(1-\exp \left(-\displaystyle \frac{\chi {E}_{\mathrm{react}}-{E}_{{\rm{b}}}}{3{E}_{{\rm{b}}}}\right)\right),\end{eqnarray} \tag{ 14 }$

where α is an empirical factor of 0.5, and the factor of 3 in the denominator arises from the assumption that only translational excitation in the direction perpendicular to the surface can lead to desorption. The free parameter χ describes the fraction of excess energy transferred to the translational excitation of the products. They proposed that χ depends on the number of product species; for the product species i, χ is given by

$\begin{eqnarray}{\chi }_{i}=\left\{\begin{array}{ll}{\chi }_{1} & (\mathrm{one}\,\mathrm{product}\,\mathrm{reaction}),\\ {\chi }_{2}\displaystyle \frac{{m}_{j}}{{m}_{i}+{m}_{j}} & (\mathrm{two}\,\mathrm{product}\,\mathrm{reaction}),\end{array}\right.\end{eqnarray} \tag{ 15 }$

where m_i is the mass of the product species i, and m_j is the mass of another species produced by the reaction. Fredon et al. (2021) conducted gas–ice astrochemical simulations using the rate-equation method in conjunction with their formalism and found that their model reasonably well reproduces the observations of gas-phase (complex) organic molecules in dark clouds when χ₁ = 0.1. They also found that the gas-phase abundances of organic molecules are not sensitive to the choice of χ₂, because the surface reactions that are efficient at low temperatures (∼10 K) involve atomic H or H₂, and thus, χ_i is small for organic molecules, which is more massive than atomic H and H₂.

Figure 11 compares P_cd for Reactions (2) and (4) constrained as described in the present work with those predicted from the formulations proposed by Fredon et al. (2021; left panel) and by Garrod et al. (2007; right panel). The binding energy of H₂S and PH₃ was set to 2700 and 2200 K, respectively (Collings et al. 2004; Molpeceres & Kästner 2021). The energy released by the reactions is listed in Table 1. The Fredon formulation predicts that P_cd(H₂S) should be lower than P_cd(PH₃), whereas the Garrod formulation predicts the opposite trend. Our results indicate that P_cd(H₂S) is slightly lower than P_cd(PH₃); thus, the Fredon formulation is more consistent with our findings. From a quantitative comparison, we find that the Fredon formulation with χ₁ ∼ 0.07 reproduces both P_cd(H₂S) and P_cd(PH₃) constrained as described in the present work; i.e., approximately 7% of the excess energy is transferred to the translational excitation of the products. The Fredon formulation predicts that, even if χ₂ is unity, the values of P_cd(HS) and P_cd(PH₂) are zero, adopting the HS and PH₂ binding energies at 2700 K and 1800 K, respectively. These results support our assumption that P_cd for the H addition reactions is higher than that for the H abstraction reactions. If we adopt the Garrod formulation to explain P_cd(H₂S) and P_cd(PH₃), the parameter a should be ∼0.05, which is larger than the typically assumed value of 0.01.

Whether the Fredon formulation with χ₁ ∼ 0.07 provides reasonable estimates of P_cd for other reaction systems is unclear. Pantaleone et al. (2020) used ab initio MD simulations to study the fate of the energy released by the reaction H + CO → HCO on crystalline water ice and found that 90% of the reaction energy is instantly injected toward the water ice. The fraction of the excess energy transferred to the translational excitation should be less than 10% in this case. However, Pantaleone et al. (2021) studied the H₂ formation on water-ice surfaces using ab initio MD simulations and found that as much as two-thirds of the reaction energy is injected toward the water ice, and the remaining energy is retained in the produced H₂. Experimental and numerical studies of additional systems are required to draw any solid conclusions.

4.4. Astrochemical Implications

Gaseous H₂S has been detected in various evolutionary stages of star and planet formation, cold dense clouds (e.g., Minh et al. 1989; Ohishi et al. 1992; Navarro-Almaida et al. 2020), envelopes around protostars (e.g., Blake et al. 1994; Wakelam et al. 2004), and protoplanetary disks (Phuong et al. 2018; Rivière-Marichalar et al. 2021). The formation of H₂S in the gas phase is inefficient because of the presence of endothermic reactions in the sequence of reaction pathways to convert S⁺ or atomic S into H₂S (e.g., Yamamoto 2017); the reaction of S⁺ with H₂ to form SH⁺ is endothermic, as is the reaction of SH⁺ (which can be formed by the reaction between atomic S and ${{{\rm{H}}}_{3}}^{+}$ ) with H₂ to form SH₂ ⁺. H₂S has been speculated to be produced on grain surfaces by the sequential hydrogenation of atomic S on grain surfaces and subsequently released to the gas phase by thermal or nonthermal desorption processes, the latter of which should dominate the former at low temperatures (∼10 K) (e.g., Garrod et al. 2007). In the molecular cloud TMC1-CP, the gas-phase H₂S abundance with respect to hydrogen nuclei is 0.6−6 × 10⁻⁹ (Navarro-Almaida et al. 2020 ). H₂S ice has not been detected in star-forming regions, and the upper limit of the H₂S/H₂O ice abundance ratio is ∼1% (Smith 1991), corresponding to the H₂S ice abundance with respect to hydrogen nuclei of ≲10⁻⁶.

Very little is known about PH₃ in star-forming regions. Neither gas-phase nor solid-phase PH₃ has been detected in star-forming regions (e.g., Turner et al. 1990; Lefloch et al. 2016). Like the gas-phase formation of H₂S, the formation of PH₃ via gas-phase reactions is inefficient (Thorne et al. 1984). Thus, the main formation pathway for gas-phase PH₃ is the formation of PH₃ ice by the sequential hydrogenation of atomic P on grain surfaces, followed by thermal or nonthermal desorption, as assumed in previous astrochemical models (e.g., Charnley & Millar 1994; Aota & Aikawa 2012; Chantzos et al. 2020; Sil et al. 2021).

To explore the effect of the chemical desorption of H₂S and PH₃ on their gas-phase abundances, we conducted gas–ice astrochemical simulations using the modified rate-equation method (Garrod 2008) under dark-cloud physical conditions; the number density of hydrogen nuclei (n_H), the temperature, and the visual extinction were set to 2 × 10⁴ cm⁻³, 10 K, and 10 mag, respectively. The cosmic-ray ionization rate of H₂ was set to 1.3 × 10⁻¹⁷ s⁻¹. In our rate-equation model, the chemistry is described by a three-phase model, where the gas phase, a surface of ice, and the chemically inert bulk ice mantle are considered (Hasegawa & Herbst 1993). As nonthermal desorption processes, which are relevant to H₂S and PH₃, chemical desorption and photodesorption were considered. The chemical desorption probability for reactions other than Reactions (2) and (4) were calculated using the method of Garrod et al. (2007), assuming a = 0.01. The photodesorption yields per incident, far-ultraviolet (FUV) photons of H₂S and PH₃, were set to 10⁻³ (see Fuente et al. 2017). Additional details can be found in Furuya et al. (2015), Furuya & Persson (2018).

Although our kMC simulations have shown that the distribution of adsorption sites with different potential-energy depths plays a role in the chemical desorption of H₂S and PH₃, we must choose single values for the binding energy and the hopping energy of atomic H in the rate-equation model. We set the hopping energy of atomic H to be 80 K, which corresponds to the low end of the distribution. The ${p}_{{{\rm{H}}}_{2}{\rm{S}}+{\rm{H}}}$ was then 3 × 10⁻³, and the ${p}_{{\mathrm{PH}}_{3}+{\rm{H}}}$ was 3 × 10⁻². The binding energy of atomic H was set to 300 K.

The amounts of elemental S and P available for gas and ice chemistry in the interstellar matter (ISM) are uncertain. In diffuse clouds, S is predominantly present in the gas phase, whereas P in the gas phase is depleted to some extent (Jenkins 2009). On the other hand, previous observations and modeling studies have suggested that the S and P abundances in star-forming regions are much lower than the values in diffuse clouds (e.g., Wakelam et al. 2004; Lefloch et al. 2016). A large fraction of elemental S and P might be incorporated into refractory compounds during the evolution from diffuse clouds to denser clouds (e.g., Bergner et al. 2019; Cazaux et al. 2022). Here, we assume that the elemental abundances of S and P available for gas and ice chemistry in the dark-cloud stage are lower than those observed in diffuse clouds by a factor 100 (Graedel et al. 1982; Wakelam & Herbst 2008). The elemental abundances of S and P with respect to H were set to 1.5 × 10⁻⁷ and 1.2 × 10⁻⁹, respectively. The elemental abundances of other elements were taken from Aikawa & Herbst (1999). Initially the species were assumed to be atoms or atomic ions except for hydrogen, which is in molecular form.

Figure 12 shows the temporal evolution of the S-bearing species (left panel) and the P-bearing species (right panel), where the P_cd(H₂S) and P_cd(PH₃) were varied. In the model with P_cd(H₂S) = P_cd(PH₃) = 0%, almost all the S and P are eventually confined in H₂S ice and PH₃ ice, respectively. The main production pathway for gas-phase H₂S and PH₃ is photodesorption of the corresponding icy molecules by cosmic-ray-induced FUV photons. In the model with P_cd(H₂S) = 3% and P_cd(PH₃) = 4%, as constrained in the present work, the gas-phase abundances of H₂S and PH₃ are higher by orders of magnitude than those in the model with P_cd = 0%. This result indicates that chemical desorption is the dominant dominant route for the supply of these gas-phase molecules. Because of chemical desorption, the abundances of H₂S ice and PH₃ ice are reduced compared with those in the model with P_cd = 0%, and the dominant reservoirs of the elemental S and P are the atomic forms.

The gas-phase H₂S abundance in the model with P_cd(H₂S) = 3% is (1–4) × 10⁻⁹ after 10⁵ yr. The predicted value is similar to that observed in the molecular cloud TMC1-CP ((0.6–6) × 10⁻⁹, depending on the position in the cloud; Navarro-Almaida et al. 2020). The gas-phase PH₃ abundance in the model with P_cd(PH₃) = 4% is (3–20) × 10⁻¹¹ after 10⁵ yr. To the best of our knowledge, PH₃ has not been detected in cold molecular clouds. Future high-sensitivity observations of the PH₃ 1₀ − 0₀ transition at 266.9445136 GHz (Müller et al. 2001, 2005) toward cold molecular clouds such as TMC-1 would provide an interesting test of the surface chemistry of P.

The model with P_cd = 1%, which is often assumed in astrochemical models, predicts the gas-phase abundances of H₂S and PH₃ similar to those in the model with P_cd(H₂S) = 3% and P_cd(PH₃) = 4%; the difference in the abundances is less than a factor of two. This weak dependence indicates that the loops of H₂S–HS and PH₃–PH₂ interconversions on dust grains are highly efficient and that the rate-limiting step of the desorption of H₂S and PH₃ in our models is the adsorption of atomic H. Notably, however, our rate-equation model, in which a single type of adsorption site is considered, may overestimate the rate of the loops of H₂S–HS and PH₃–PH₂ interconversions. In reality, deep-potential sites can trap atomic H, slowing the loops of H₂S–HS and PH₃–PH₂ interconversions, as observed in our kMC models. To accurately evaluate the effect of chemical desorption under the ISM conditions, the binding energy distribution for atomic H should be considered by solving a gas and ice chemical network.

5. Summary

Chemical desorption is caused by the energy released by reactions; after an exothermic surface reaction, some of the excess energy is dissipated into product's translational energy in the direction perpendicular to the surface, leading to desorption (e.g., Fredon et al. 2017). Chemical desorption is usually included in modern gas–ice astrochemical models after Garrod et al. (2007); however, its efficiency is poorly constrained, especially desorption from water ice. Oba et al. (2019) experimentally studied chemical desorption upon the reaction of H₂S with H atoms on porous ASW. Nguyen et al. (2021) conducted similar laboratory studies for the reactions of PH₃ with H atoms. These studies demonstrated that H₂S and PH₃ can be lost from water-ice surfaces by chemical desorption. They also estimated a chemical desorption probability for H₂S and PH₃ per reactive species (i.e., per incident H atom) of ∼1% on porous ASW. As noted by Oba et al. (2018), the desorption probability per incident of H atom corresponds to the lower limit of the desorption probability per reactive event, which the astrochemical models require as inputs, because a substantial fraction of the adsorbed H atoms on ASW surfaces would be thermally desorbed and subsequently consumed by the H₂ formation reaction on the surface. In the present work, we constrained the desorption probability of H₂S and PH₃ per reactive event on porous ASW by numerically simulating the laboratory experiments of Oba et al. (2019) and Nguyen et al. (2021). We used kinetic Monte Carlo simulations in which the position and movement of each chemical species on surfaces were tracked over time and in which the binding energy distributions of atomic H and H₂ were considered. Our findings are summarized as follows.

1.
The chemical desorption probability of H₂S and PH₃ per hydrogenation event of the precursor species on porous ASW are constrained to 3% ± 1.5% and 4% ± 2%, respectively.
2.
These probabilities are consistent with a theoretical model of chemical desorption proposed by Fredon et al. (2021), with χ ∼ 0.07% (where χ is a free parameter) describing the fraction of reaction energy transferred to the translational excitation of reaction products. Whether the Fredon et al. model with χ ∼ 0.07 provides reasonable estimates of P_cd for other reaction systems is unclear. Experimental and numerical studies of additional systems are required for better understanding of the chemical desorption process.
3.
As a byproduct, we constrained the binding energy distribution of atomic H. The abundance of H₂S became ∼40% and ∼70% of the initial H₂S abundance in the 20 and 30 K experiments by Oba et al. (2019), respectively. These results indicate that ∼70% of the adsorption sites on porous ASW should have E_b(H) ≲ 300 K, whereas ∼40% of adsorption sites should have E_b(H) ≲ 200 K (see Section 3.1.1).

Finally, we stress that, for the chemical reactions systems on which the present work is focused, basic chemical data are available in the literature (e.g., rate coefficients and binding energies). Without these basic data provided by computational chemistry and experiments, we could not have constrained the chemical desorption probability. Basic chemical data are critical not only for the astrochemical models of astrophysical objects but also for extracting the chemical parameters from laboratory experiments of surface reactions.

K.F. acknowledges Naoko Yokomura and Yuri Aikawa for their assistance in developing the kMC code used in this work. Y.O. acknowledges Naoki Watanabe, Akira Kouchi, Hiroshi Hidaka, and Nguyen Thanh for discussions of the experimental chemical desorption results. We are grateful to the anonymous referees for providing valuable comments that helped improve the manuscript. This work is partly supported by JSPS KAKENHI grant numbers 17H06087, 20H05847, 21H01145, 21H04501, and 21K13967. Numerical computations were carried out in part on a PC cluster at the Center for Computational Astrophysics, National Astronomical Observatory of Japan.

Appendix A: Additional Model Results

Figure 13 shows the results from additional models of the H₂S + H system, where ${E}_{\ {\rm{b}}}^{\min }({{\rm{H}}}_{2})$ is fixed at 150 K, and f and P_cd(H₂S) are varied between 0.6 and 0.8 and between 5% and 20%, respectively. Among the nine different models shown in the figure, only three models can reasonably well reproduce the 30 K experiments (shown by thick lines in the figure): the model with f = 0.7 and P_cd(H₂S) = 20%, the model with f = 0.8 and P_cd(H₂S) = 5%, and the model with f = 0.8 and P_cd(H₂S) = 10%. The first and third models overestimate the amounts of H₂S and HS desorbed from the surface at 10 and 20 K compared with the amounts measured in the corresponding experiments. The second model moderately well reproduces the experiments at 10, 20, and 30 K simultaneously; however, the model tends to overestimate the amount of H₂S and HS desorbed from the surface compared with the experiments at 10 and 20 K. Notably, this model is similar to our best-fit model (f = 0.8 versus 0.85 and P_cd(H₂S) = 5% versus 3%); however, the reduced χ² for this model is greater than that for the best-fit model (2.38 versus 1.85).

Appendix B: Polynomial Distribution versus Gaussian Distribution

In our models, we assumed that the binding energy distributions for H₂ are described by a polynomial function (Equation (9)) from 780 K to ${E}_{\ {\rm{b}}}^{\min }$ on the basis of the experiments by Amiaud et al. (2006, 2015). The distribution at E_b lower than ∼350 K was not constrained well in the experiments. We simply extrapolated the polynomial function to the lower E_b values. Here, we check the effect of this extrapolation on the modeling results. We tested some models in which the H₂ binding energy followed a Gaussian distribution with a mean value of 270 K and a full-width at half-maximum of 190 K. The lower boundary of E_b(H₂) was assumed to be 150 K, which roughly corresponds to the intermolecular potential energy of a H₂–H₂O dimer (∼120 K; Zhang et al. 1991). As shown in Figure 14, this Gaussian distribution matches the polynomial distribution for E_b ≳ 350 K; by contrast, for lower E_b, the two distributions deviate.

**Figure 14.** Binding energy distribution for H₂ assuming a Gaussian distribution (red) and a polynomial function (black). See the text for additional details.
Download figure:
Standard image High-resolution image

Figure 15 shows the results of a small grid of models, where f and P_cd(H₂S) are varied between 0.6 and 0.8, and between 3% and 10%, respectively. As shown, no model can reproduce the experiments at 10, 20, and 30 K simultaneously. We conclude that there is no clear advantage to using the Gaussian distribution instead of the polynomial distribution.

Quantifying the Chemical Desorption of H₂S and PH₃ from Amorphous Water-ice Surfaces

Article metrics

Author e-mails

Author affiliations

ORCID iDs

Dates

Abstract

1. Introduction

2. Methods

2.1. Numerical Setup

2.2. Chemical Parameters

2.2.1. Sticking Coefficient

2.2.2. Binding Energy

2.2.3. Chemical Desorption Probability

3. Results

3.1. H₂S + H

3.1.1. Fiducial Model

3.1.2. Lower Limit of H₂S Desorption Probability

3.1.3. Constraining the H₂S Desorption Probability

3.2. Constraining the PH₃ Desorption Probability

4. Discussion

4.1. Uncertainties in the Flux of Atomic H

4.2. Diffusion of Atomic H by Quantum Tunneling

4.3. Comparison with Theoretical Models of Chemical Desorption

4.4. Astrochemical Implications

5. Summary

Appendix A: Additional Model Results

Appendix B: Polynomial Distribution versus Gaussian Distribution

Quantifying the Chemical Desorption of H2S and PH3 from Amorphous Water-ice Surfaces

Article metrics

Share this article

Author e-mails

Author affiliations

ORCID iDs

Dates

Abstract

1. Introduction

2. Methods

2.1. Numerical Setup

2.2. Chemical Parameters

2.2.1. Sticking Coefficient

2.2.2. Binding Energy

2.2.3. Chemical Desorption Probability

3. Results

3.1. H2S + H

3.1.1. Fiducial Model

3.1.2. Lower Limit of H2S Desorption Probability

3.1.3. Constraining the H2S Desorption Probability

3.2. Constraining the PH3 Desorption Probability

4. Discussion

4.1. Uncertainties in the Flux of Atomic H

4.2. Diffusion of Atomic H by Quantum Tunneling

4.3. Comparison with Theoretical Models of Chemical Desorption

4.4. Astrochemical Implications

5. Summary

Appendix A: Additional Model Results

Appendix B: Polynomial Distribution versus Gaussian Distribution

Quantifying the Chemical Desorption of H₂S and PH₃ from Amorphous Water-ice Surfaces

3.1. H₂S + H

3.1.2. Lower Limit of H₂S Desorption Probability

3.1.3. Constraining the H₂S Desorption Probability

3.2. Constraining the PH₃ Desorption Probability