Modelling the Radical Chemistry on Ice Surfaces: An Integrated Quantum Chemical and Experimental Approach

Heterogeneous radical processes on ice surfaces play a vital role in the formation of building blocks of the biologically relevant molecules in space. Therefore, quantitative mechanistic details of the radical binding and radical reactions on ices are crucial in rationalizing the chemical evolution in the Universe. The radical chemistry on ice surfaces was explored at low temperatures by combining quantum chemical calculations and laboratory experiments. A range of binding energies was observed for OH, HCO, CH3, and CH3O radicals binding on ices. Computed reaction paths of the radical reactions on ices, OCS + H and PH3 + D, explained the experimentally observed products. In both radical reactions, quantum tunnelling plays a key role in achieving the reactions at low temperatures. Our findings give quantitative insights into radical chemistry on ice surfaces in interstellar space and the planetary atmospheres.


INTRODUCTION
The mechanisms of the heterogeneous chemical reactions on dust and ice surfaces in space are essential in understanding the molecular evolution in the Universe that leads to the origin of life. A number of molecular species have been detected in interstellar space (Herbst and van Dishoeck, 2009;Guélin and Cernicharo, 2022;Herbst and Garrod, 2022) and the planetary atmospheres. (Zapata Trujillo et al., 2021) However, their origin and chemical reactions, giving rise to complex organic molecules (COM), are not fully understood. The most plausible theory for COMs formation in interstellar circumstances is the radical-driven chemistry on ice dust particles. (Garrod and Herbst, 2006;Garrod et al., 2009;Chuang et al., 2016) At the early stage of the COMs formation, the hydrogen atoms addition reactions on the primordial species, such as C, N, O, and CO is occurred on ices at very low-temperatures (~10 K). The hydrogen addition processes give rise to various small radical species (e.g., OH, HCO, CH 3 O, NH, NH 2 , CH, CH 2 , and CH 3 ) (Herbst and van Dishoeck, 2009) and stable molecules (e.g., methanol, formaldehyde, ammonia, methane). As the temperature rises, i.e., during the warming-up stage, small molecules or radicals begin to diffuse on ice to form COMs when they encounter. Thus, the elementary chemical processes on ices, such as radical adsorption, radical diffusion, and radical reaction mechanisms are required to rationalize the mechanisms of the COMs formation. However, radical processes on ices are challenging to characterize from laboratory experiments alone. For example, radical species are highly reactive and consume quickly, and therefore the number density of the radical species is too low for detection. Also, the temperature-programmed desorption method cannot be used to determine the binding energies of radicals on ices. Therefore, computational methods employing ab-initio calculations become critical in quantitative determination of the radical processes on ices. (Sameera et al., 2017;Enrique-Romero et al., 2019;Lamberts et al., 2019;Duflot et al., 2021) The quantum chemical methods have been the most successful approaches in calculating molecule structures and their properties. High-level ab-initio calculations, configuration interaction (CI), multireference CI, coupled-cluster singles doubles perturbative triples CCSD(T), (Pople and Raghavachari, 1987) for instance, can be performed for a small molecular system with few atoms. For large molecular systems, density functional theory (DFT) would be the method of choice. The quantum mechanics/molecular mechanics (QM/MM) methods can be used for modelling relatively large molecular systems. (Sameera and Maseras, 2012;Chung et al., 2015;Sameera et al., 2017;Sameera and Maseras, 2018) In a QM/ MM calculation, the electronically important part of the molecular system is calculated using a QM method, and the remaining part of the molecule is calculated by a MM method. Therefore, QM/MM calculations give accurate results at a low computational cost. Chemical processes on ices can be computed using DFT or DFT/MM methods, employing ice cluster models or periodic ice structures.
The first step of a computational study is to calculate the potential energy surface (PES) of the molecular system, which is a geometric hypersurface representing the potential energy as a function of the coordinates of the molecular system. The PES consists of stationary points, specifically local minima (LM, i.e., reactant, intermediates, product) or first-order saddle points (i.e., transition states, TS). A simplified potential energy surface of a single-step reaction is shown in Figure 1. In this PES, the LM1 represents the reactant(s), and LM2 is the product(s) of the reaction. The LM1 and LM2 are connected through the TS1. The molecular structures of LM1, LM2, and TS1 can be optimized using the ab-initio computations. (Schlegel, 2011) Also, the connectivity between the LM can be confirmed, for instance, by employing the intrinsic reaction coordinate (IRC) calculations. (Fukui, 1981) Then, the reaction barrier (i.e., the relative energy difference between LM1 and TS1) and the reaction rate can be calculated.
Multistep chemical reactions, consisting of several chemical steps, can proceed on ices in interstellar space and planetary atmospheres. PES of a multistep reaction is challenging to calculate, as it goes through several LM and TS ( Figure 1B). There are two methods for calculating reaction paths of multistep chemical reactions, specifically conventional or systematic approaches. The conventional method uses chemical knowledge and intuition to sketch the possible reaction mechanisms. (Sameera and Maseras, 2012) Then, PESs of the guessed reaction paths are calculated. After that, the lowest energy path connecting the reactant and the product, the socalled reaction mechanism, can be determined. The conventional approach is the most common method in calculating reaction mechanisms. However, the accuracy of the computed reaction mechanism depends on the guessed reaction paths. On the other hand, systematic methods determine the PESs of all possible reaction paths (known, unknown, or unexpected) in an automated way. (Sameera et al., 2016) Therefore, systematic approaches are vital for investigating complex chemical networks. As the systematic methods explored reaction paths in a broad sense, the computational cost of a systematic reaction path survey is substantial.
Once the reaction paths are determined from conventional or systematic methods, PESs of the reaction mechanisms can be prepared. Then, reaction rates can be calculated using the transition state theory (TST). (Eyring, 1934) Quantum tunnelling is very important for the reactions if the system temperature is very low ( Figure 1C). The classical TST does not include tunneling effects. However, a semiclassical approach can be used for including tunneling effects. Nyman, 2021) After determining the PESs for the radical process on ices from quantum chemical methods, computed data must be compared with the available experimental results. Then, an accurate picture of the radical chemistry on ices can be achieved. This mini-review focuses on our recent progress on radical binding and radical reaction mechanisms on ices. The main focus of the article is the computational results and their implications. Relevant experimental results were described briefly to compare with the computational results.

COMPUTATIONAL METHODS
LM and TS on the ground state PESs were optimized using QM or QM/MM methods as implemented in the Gaussian16 program. (Frisch et al., 2016) The two-layer Our own N-layered Integrated molecular Orbital and molecular Mechanics (ONIOM) (Chung et al., 2015;Sameera and Maseras, 2018) method was employed for QM/MM calculations. We have used density functional theory (DFT) and standard double-zeta (Ditchfield et al., 1971;Hehre et al., 1972) or triple-zeta basis sets (Schäfer et al., 1992;Schäfer et al., 1994) for QM calculations. The ωB97X-D (Chai and Head-Gordon, 2008) or M062X (Zhao and Truhlar, 2008) density functionals were employed for QM calculations. The AMOEBA09 polarizable force field was used for MM computations, employing the SICTWO interface. (Sameera and Maseras, 2018) Vibrational frequency calculations were performed to confirm the nature of the stationary points (i.e., no imaginary frequency for LM and one imaginary frequency for TS) and to calculate zero-point energies. Connectivity between the LM was confirmed by performing IRC calculations. (Fukui, 1981) Excited-state PESs were calculated using the time-dependent DFT method. The energy decomposition analysis (EDA) (Morokuma, 1971;Ziegler Rauk, 1977) was performed using the ADF program. (te Velde et al., 2001) The excited states seam crossing calculations were performed using the GRRM program. (Maeda et al., 2018) We have used the ice structure models of Andersson et al. (Andersson et al., 2006) to make ice cluster models for crystalline hexagonal water ice (I h ) and amorphous solid water (ASW). In ONIOM(QM:MM) structure optimizations, the H 2 O molecules in the MM region were frozen (i.e., kept as in the ice structure models of Andersson et al. (Andersson et al., 2006)).

RADICAL BINDING ON ICES
The binding energy of the radical species on ices gives some insights into their desorption and diffusion. These properties are vital in understanding radical reactions on ices. Radical species interact with the dangling-H (d-H) and dangling-O (d-O) atoms on ice surfaces ( Figure 1D). To understand the roles of the dangling atoms on the binding energy, binding energy of OH, HCO, and CH 3 radicals on I h ( Figure 1E (Nguyen et al., 2021a) and OCSH: 0.19-0.46 eV (2205-5338 K), (Nguyen et al., 2021b).
Understanding the origin of the strong or weak radical binding energies on ices is essential. During the radical binding on ices, the radical and ice structures can be deformed that cost some energy. Also, the interactions between the radical and the ice surface are vital for the binding energy. We have performed an energy decomposition analysis (EDA, Figure 1H) (Morokuma, 1971;Ziegler Rauk, 1977) to rationalize the deformation energy and interaction energy of CH 3 O radical binding on ASW. (Sameera et al., 2021) The energy difference between ice-CH 3 O complex and the energy sum of the isolated ice structure (A) and CH 3 O (B), without changing their internal geometries, is defined as the interaction energy, INT. The INT term was further decomposed into electrostatic interactions (ΔE elstat ), Pauli repulsion (ΔE Pauli ), and orbital interactions (ΔE Oi ). The deformation energy (DEF) is defined as the energy sum of the isolated ice and CH 3 O structures at the ice-CH 3 O complex and the energy sum of the fully optimized ice and CH 3 O structures (i.e., A 0 and B 0 ). According to EDA, the deformation energy of the CH 3 O radical or ice is almost zero. (Sameera et al., 2021) Thus, the radical bonding process does not require geometric or electronic preparation. On the other hand, the interaction energy between CH 3 O and ASW is significant, and therefore the interaction energy controls the strength of the binding energy. According to the computed data, the ΔE elstat and ΔE Oi overcome the ΔE Pauli , giving rise to strong or weak binding energy ( Figure 1I).
The physicochemical behaviour of OH radicals on ices depends on the adsorption sites. (Miyazaki et al., 2020;Tsuge and Watanabe, 2021) Recent experiments observed photodesorption of OH radicals on ASW upon absorption of 532 nm photons. However, 532 nm photons cannot be absorbed by isolated OH and H 2 O molecules. Computed ground and excited-state PESs of the OH radical on ices suggested that one-photon absorption of the OH-(H 2 O) n complex ( Figure 1G), having strong hydrogen bonds between OH and surrounding H 2 O molecules, leads to OH desorption. Further, the excited OH radical enters a dissociation channel through the conical intersection (CI) near the excited state minima ( Figure 1G), leading to the OH desorption from the ASW with the kinetic energy (KE) of 0.24 eV, which is in agreement with the experimentally determined KE.

RADICAL REACTIONS ON ICES
Radical reactions in interstellar space and the planetary atmospheres produce small molecules or radical species. Thus, reaction mechanisms of the radical reactions on ices are crucial in understanding the chemical evolution in the Universe. We have performed experimental and quantum chemical studies side-byside to rationalize the radical reactions on ices. The following section focuses on the characterization of two radical reaction mechanisms on ices.

Carbonyl Sulphide + H
Chemistry of the sulphur (S)-bearing species in space has been an active topic in the astrochemical community. (Oba et al., 2018;Laas and Caselli, 2019;Cabezas et al., 2021) A number of S-bearing species have been detected in space. Carbonyl sulphide (OCS), a S-bearing molecule, was detected in the gas (Goldsmith and Linke, 1981) and solid phases of the ISM. (Palumbo et al., 1995) According to the laboratory experiments under the ISM conditions, OCS can be produced through the photon and cosmic ray-induced processes of H 2 S-containing ices. (Palumbo et al., 1995;Ferrante et al., 2008;Jiménez-Escobar et al., 2014) The reactivity of OCS is, however, not well established. Thus, we have performed experiments on the surface reaction of OCS solid with H atoms on icy surfaces at 10 K. (Nguyen et al., 2021b) We found that the addition of H atoms to OCS solid proceeded via quantum tunnelling. The main products were H 2 S, CO, H 2 CO, and CH 3 OH. Thioformic acid (HC(O)SH) was also observed as a minor product. Reaction paths for the formation of H 2 S, CO, and HC(O)SH are unknown. The H 2 CO and CH 3 OH molecules may be formed by successive H addition to CO. (Watanabe and Kouchi, 2002) Quantum chemical calculations were performed to rationalize the reaction mechanisms. Calculated potential energy surfaces are shown in Figure 2A.

Phosphine + D
Phosphine (PH 3 ) is an important phosphorus-bearing molecule in space. PH 3 was detected in the atmosphere of Jupiter, Saturn, and Venus. (Sousa-Silva et al., 2020;Greaves et al., 2021) The formation of PH 3 in the gas phase of the ISM is a matter of controversy. (Thorne et al., 1984;Millar, 1991) Therefore, the chemistry of PH 3 in space is very important in understanding its origin and the chemical evolution of the phosphorus-bearing species in space. We have performed a combined experimental and computational study to rationalize the reactions between PH 3 and D on ASW at low temperatures. According to our experimental data, PD 3 is formed on ice and partly released into the gas phase by chemical desorption. We have also detected PH 2 D and PHD 2 , but their concentration is relatively low compared to PD 3 since their deuteration is fast on icy surfaces. (Nguyen et al., 2021a) However, the mechanism of the reaction between PH 3 and D is unknown. Thus, we have performed quantum chemical calculations to rationalize the mechanism.
Computed reaction pathways for the reaction between PH 3 and D are shown in Figure 2B. The computed reaction barrier for the PH 3 + D → PH 2 + HD reaction is 0.13 eV (1509 K). After the formation of PH 2 radical, two reaction paths are possible; 1) PH 2 + D → PH 2 D and 2) PH 2 + D → PH + HD. The former reaction is barrierless, while the latter reaction has a barrier of 0.11 eV (1277 K). Thus, PH 2 D is the main product of the reaction. Similar vein, PH 2 D + D → PHD + HD and PHD 2 + D → PD 2 + HD reactions occur through the reaction barriers of 0.13 eV (1509 K). Finally, a barrierless reaction between PD 2 radical and D leads to PD 3 . It is important to note that we have performed electronic structure calculations, where zeropoint energy of the systems was calculated introducing the mass of H and D. Computed reaction paths explained the mechanism for the formation of PD 3 that goes through the experimentally characterized PH 2 D and PHD 2 intermediates. Our findings give implications to understand chemical networks that include phosphine in interstellar space and planetary atmospheres.

SUMMARY
Modern computational methods, employing quantum chemical computations, allow quantitative determination of the heterogeneous radical processes on ice surfaces, such as the binding energies of radicals on ices and radical reaction mechanisms. Thus, theory and computations play a crucial role in modern astrochemistry. Computed binding energies of the radical species are summarized in Table 1. In all cases, a range of binding energies were found. Thus, a distribution of binding energies would be used in developing astrochemical models.
Reported binding energies of Ferrero et al. (Ferrero et al., 2020) [OH (1551-5321 K), HCO (1315-3081 K), CH 3 (1119-1654 K)] and Enrique-Romero et al. (Enrique-Romero et al., 2022)  ] are also in a range, where a few binding sites of a large water ice cluster model were taken into account for calculating binding energies. Wakelam et al. (Wakelam et al., 2015) also reported a range of binding energies for OH (3300-5300 K) and HCO (2300-2700 K) using small water ice cluster models. Compared to their computed binding energies, we have found a broader range of binding energies due to the fact that our studies employed a number of binding sites in large water ice clusters. We have rationalized the radical reaction mechanisms on ice surfaces by combining laboratory experiments and quantum chemical calculations. According to the calculated potential energy surfaces for the reaction between OCS and H on ices, the first hydrogenation gives rise to OCS-H, where the reaction barrier is relatively high. Therefore, quantum tunnelling is critical. The second hydrogenation yields CO and H 2 S as the major products and HC(O)SH as the minor product. A recent astronomical observation indicated the presence of HC(O)SH in space.
Computed reaction mechanisms suggested that the reaction between PH 3 and D yields PD 3 , where step-wise dehydrogenation is occurred with a relatively high reaction barrier. Thus, quantum tunnelling is critical for the mechanism to operate at very low temperatures. Deuterated isotopologues of phosphine have not been detected in astronomical observations. However, based on our experimental and computational studies, we argue that the deuterated isotopologues of phosphine can be formed on ices. Thus, we suspect that astronomical observations would detect the deuterated isotopologues of phosphine in the near future.

AUTHOR CONTRIBUTIONS
WMCS made the first draft of the manuscript. All authors revised the manuscript and were approved the final version of the manuscript.

FUNDING
This work was partly supported by JSPS KAKENHI Grant Numbers JP19K03940 and JP21H05416 (to WS), JP17H06087 (to NW), JP21H04501 (to YO), and JP21K13974 (to TN).