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₃O, NH, NH₂, CH, CH₂, and CH₃) (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 so-called 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. (Senevirathne et al., 2017; 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.

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₂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)).

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₃ radicals on I_h (Figure 1E) was calculated. (Sameera et al., 2017) Computed binding energies of OH (0.20–0.67 eV, 2321–7775 K), HCO (0.12–0.47 eV, 1393–5454 K), and CH₃ (0.11–0.26 eV, 1277–3017 K) on I_h indicated a range of binding energies. According to the calculated averaged binding energies, binding preference follows the order of OH (0.46 eV, 5338 K) > HCO (0.30 eV, 3481 K) > CH₃ (0.20 eV, 2321 K). Recently, both I_h and ASW models (Figure 1F) were employed to calculate the binding energy of OH on ices. (Miyazaki et al., 2020) Moreover, a broader range of binding energies (0.06–0.74 eV, 696–8587 K) was observed relative to the previous work that employed only I_h models (0.20–0.67 eV, 2321–7775 K). The computed binding energies of the other radical species are also in a broad range, CH₃O: 0.10–0.50 eV (1160–5802 K), (Sameera et al., 2021) PH₂: 0.13–0.21 eV (1503–2406 K), (Nguyen et al., 2021a) PH: 0.10–0.19 eV (1105–2205 K), (Nguyen et al., 2021a) P: 0.07–0.16 eV (831–1800 K), (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₃O radical binding on ASW. (Sameera et al., 2021) The energy difference between ice-CH₃O complex and the energy sum of the isolated ice structure (A) and CH₃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₃O structures at the ice-CH₃O complex and the energy sum of the fully optimized ice and CH₃O structures (i.e., A₀ and B₀). According to EDA, the deformation energy of the CH₃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₃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₂O molecules. Computed ground and excited-state PESs of the OH radical on ices suggested that one-photon absorption of the OH-(H₂O) _n complex (Figure 1G), having strong hydrogen bonds between OH and surrounding H₂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 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-by-side to rationalize the radical reactions on ices. The following section focuses on the characterization of two radical reaction mechanisms on ices.

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₃ (1119–1654 K)] and Enrique-Romero et al.(Enrique-Romero et al., 2022) [CH₃O (3139–4582 K)] 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₂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₃ and D yields PD₃, 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.