Reconstruction of water ice: the neglected process OH + OH → H 2 O + O

Context. Although H 2 O is the most important molecular material found in the solid state in the interstellar medium, the chemical routes leading to ice through surface reactions are still a matter of discussion. Three reaction pathways proposed in the past are at the heart of current research: hydrogenation of atomic oxygen, molecular oxygen, and ozone. The reaction network ﬁnally leads to a small number of processes giving H 2 O: H + OH, H 2 + OH, and H + H 2 O 2 . To these processes, OH + OH should be added. It is known to be efﬁcient in atmospheric chemistry and takes the irradiations of the interstellar grains into account that, directly or indirectly, create a number of OH radicals on and in the icy mantles. Aims. We study the role of the existing ice in its own reconstruction after it is destroyed by the constant irradiation of interstellar grains and focus on the OH+OH reaction in the triplet state. Methods. We used numerical simulations with a high level of coupled cluster ab initio calculations for small water aggregates and methods relevant to density functional theory for extended systems, including a periodic description in the case of solid water of inﬁnite dimensions. Results. OH + OH → H 2 O + O reaction proﬁles are reported that take the involvement of an increasing number of H 2 O support molecules into account. It is found that the top of the barrier opposing the reaction gradually decreases with the number of supporting H 2 O and falls below the level of the reactants for H 2 O layers or solid water. Conclusions. In contrast to the gas phase, the reaction is barrierless on water ice. By adding a reconstructed H 2 O molecule and a free oxygen atom at the surface of the remaining ice, this reaction leaves open the possibility of the ice reconstruction.


Introduction
Although water ice has been found throughout the cold interstellar medium (ISM; Merrill et al. 1976;Whittet et al. 1988;Dartois et al. 1998;Gibb et al. 2004;Boogert et al. 2008), the chemical origin and longevity of the icy grain mantles in which H 2 O is the dominant component is still debated. The observed abundance of water ice cannot be explained by a simple freezeout of H 2 O onto the grain surface, and the persistent survival of the mantles under an intense radiation bombardment raises a question. Because the molecule formation in the gas phase is not very efficient, on-grain processes have to play a major role in the survival of the ice mantles. These on-grain mechanisms have been proposed in chemical models in the 1980s and are still discussed today (Tielens & Hagen 1982;d'Hendecourt et al. 1985;Hasegawa & Herbst 1993;Cuppen & Herbst 2007). The three reaction pathways proposed by Tielens & Hagen (1982), that is, hydrogenation of atomic oxygen, molecular oxygen, and ozone, are not only at the origin of the chemical models but also inspired most of the laboratory experiments. The last steps of the reaction network, which are those that effectively yield H 2 O, are Reactions (1) and (2) appear as the final steps in the hydrogenation of the oxygen atom; they also appear in the hydrogenation of ozone after OH has been formed by Reaction (3) is effective in the hydrogenation of molecular oxygen after H 2 O 2 has been obtained following H + O 2 → HO 2 H + HO 2 → H 2 O 2 .
All these mechanisms have been investigated in the laboratory. For instance, experimental evidence for water formation through ozone hydrogenation (4) onto dust grain analogs has been obtained by Mokrane et al. (2009). The water formation by H hydrogenation of an O 2 surface has been probed by Ioppolo et al. (2010) and Cuppen et al. (2010). These studies also confirmed the implication of reaction (3) along the way (Hiraoka et al. 1998;Kouchi et al. 2008). However, although presumed to exist in the icy mantle of interstellar grains, neither O 2 nor O 3 has been identified so far.
On the other hand, the OH radical appears at one stage or another in all mechanisms of H 2 O formation, and we recall that a number of these radicals exist on the ice cover or around it as the incident cosmic irradiation breaks the H 2 O bonds, which A&A 638, A125 (2020) produces a copious amount of radicals such as O, H, and OH. Here we focus on the OH + OH → H 2 O + O (7) reaction that has not been considered properly before, probably because of an activation in the gas phase (∼2 kcal mol −1 , according Karkach & Osherov 1999) and processes along a triplet surface reaction, which is not the reaction with the lowest energy. The aim of this research is to assess the possible role of the ice in the formation, or more precisely, the re-formation, of the ice mantle. In order to lead to the formation of a water molecule from two interacting OH radicals, reaction (7) has to proceed in a triplet state, regardless of the surface underneath (the end products, H 2 O( 1 A 1 ) and O( 3 P), are both in their ground states). The other possible reaction starts in the singlet state of the [HO...OH] complex, which implies the coupling of the two OH radicals and leads to hydrogen peroxyde H 2 O 2 . This is only one more way to start reaction (3), which is taken into account in existing chemical models. Although the oxygen atom does not participate as a reactant in the reaction OH + OH → H 2 O + O and appears only as a final product when the reaction is completed, it might be wondered whether the presence of water molecules would reverse the triplet/singlet stability ordering of atomic oxygen (Molpeceres et al. 2019). In the gas phase, atomic spectroscopy has shown that the singlet state is less stable than the triplet by 1.97 eV (Kramida et al. 2019). A test of atomic oxygen on a layer of 24 H 2 O molecules cut in solid water ice showed a singlet-triplet gap of ∼2 eV at the MPWB1K level (Zao & Truhlar 2004). Our density functional theory (DFT) periodic calculations (often referred to as "first principle" treatments) on a solid surface of infinite dimension give ∼3 eV. This great difference in energy means that a triplet-to-singlet swap is unlikely at the end of a reaction that occurs entirely on the triplet potential surface.
Here we present an extensive survey of reaction (7). This survey starts from a gas-phase situation intended to serve as reference, continues in the presence of an increasing number of H 2 O molecules; then a bilayer is considered that mimics crystalline ice, and we conclude with a real solid surface of infinite dimension. We summarize our conclusions in Sect. 7.

Theoretical background
We investigate the possible role of ice in its own reconstruction through reaction (7). To achieve this, it is necessary to accurately describe the chemical process at the molecular level for simple species such as OH radicals or small aggregates and simultaneously at the level of the solid water that supports the reaction. These are two different domains for which elaborate and conceptually different theories have been developed to account for changing electronic correlation/dispersion effects in the course of chemical reactions. It is particularly important when weak interactions are involved. In both cluster and solid-state approaches, the reconstruction of water ice has been considered to proceed in an overall triplet state. In other words, all structures, from the pre-reactant complex [OH. . .OH] and transition state to the final products H 2 O + O( 3 P), are calculated in a way that satisfies the triplet multiplicity.

Cluster approach
In ab initio post-Hartree-Fock (HF) methods designed for small-dimension systems (here simple models where the ice is represented by small H 2 O aggregates), the effects of electronic correlation are treated using a systematic progression in the number of electronic excitations from the HF reference state by means of perturbative treatments for double excitations (MP2) and by high-level coupled cluster calculations for double and triple excitations (CCSD and CCSD(T); Bartlett & Shavitt 1977;Raghachavari et al. 1989). Atomic basis sets of increasing levels of flexibility, from 3-21+G(d,p) to correlation consistent cc-pVnZ of double, triple, and quadruple (n = 2-4) quality, possibly augmented with diffuse functions, were employed to perform the calculations.
Calculations for much larger systems (extended models where the ice is represented by aggregates of dozens of H 2 O molecules) are the domain of the DFT. We considered the hybrid formalisms that are used most frequently, B3LYP and BHandHLYP, as well as the functionals that were recently improved to account for dispersion effects, PBE1PBE, MPW1K, and MPWB1K (Lynch et al. 2000;Zao & Truhlar 2004). The basis sets were the same in ab initio and DFT calculations.
We note a last technical point: the zero-point vibrational energies (ZPE) were included in every building of the potential energy surfaces. Because analytical frequencies cannot be calculated at the CCSD(T) level in the code we employed, we instead used the force constants that are available at this level by double numerical differentiation. Hence, harmonic vibrational frequencies were calculated at this level numerically by the standard procedure. For all other single points, the vibrational frequencies were calculated analytically at the level at which the geometries were optimized. We performed these calculations using levels of theory and basis sets as implemented in the GAUSSIAN09 package (Frish et al. 2009).

Solid-state approach
The approach presented above might suffer, a priori, from inherent dimension constraints that might introduce some distortion in the results. To avoid these limiting factors, a solution is to set up the icy environment with a proper solid surface of infinite dimension. Based on our preceding works on molecular adsorption on water ices (Lattelais et al. 2011(Lattelais et al. , 2015Bertin et al. 2017), we took a periodic representation of solid hexagonal ice I h , composed of bilayers of water molecules. The ice surface was described by the (0001) Miller indices. This surface was shown to be the most stable (Pan et al. 2010). This surface also satisfies the Bernal & Fowler (1933) constraint: exactly two hydrogens of neighboring water molecules form hydrogen bonds with oxygen lone pairs of the adjacent water molecules. More details can also be found in previous work on the structures of ice (Casassa & Pisani 2002;Calatayud et al. 2003;Casassa et al. 2005).
The electronic treatment relies on the generalized gradient approximation (GGA) in the form of the PBE functional. Long-range van der Waals and hydrogen-bonding interactions are better described within the corrective scheme (PBE+D2 ) of Grimme et al. (2010). The core electrons were kept frozen and were replaced by pseudo-potentials generated by the planeaugmented wave (PAW) method.
Cluster calculations show that the functional must include the exact HF exchange. We performed solid-state calculations of the molecular two-body OH + OH → H 2 O + O case system using a cell large enough to prevent any interaction between the molecules in adjacent cells. To make the solid-state description comparable to that of the MPWB1K hybrid DFT approaches we employed for cluster calculations, a 50% percentage of the exact HF exchange was added to the PBE+D2 functional.  Notes. E e is the electronic energy, and E 0 is the electronic + zero-point vibrational energy.
These solid-state periodic DFT treatments, also known as first-principle calculations, were carried out using the Vienna ab initio simulation package (VASP; Kresse & Hafner 1994;Kresse & Furthmüller 1996;Kresse & Joubert 1999). More details on the computational procedure using VASP are given in the appendix.

Two-body OH + OH reaction as a case study
An extensive study of the reaction profile of reaction (7) was carried out. We anticipated calculations on much larger systems and therefore used multiple levels of wave functions, from highly correlated post-HF to DFT treatments using a large panel of functionals.
The reaction pathway illustrated in Fig. 1 from the values listed in Table 1 shows the same stationary points at all levels of theory. A shallow complex [OH...OH] forms in the entrance channel, followed by a transition state whose position with respect to the starting species is largely dependent on the level of theory, and another shallow complex in the exit channel [O...H 2 O] that is close in energy to the isolated products H 2 O and O( 3 P).
We found that the entire reaction proceeds in the same plane, in agreement with a previous study (Karkach & Osherov 1999) at the QCISD(T)/6-311(d,p) level. The reaction is exothermic with a value of approximately −19 kcal mol −1 at the CCSD(T) and QCISD(T) levels. This value is reduced to approximately −16 kcal mol −1 for coupled cluster calculations when zero-point energies are taken into account. At the same level of theory, the entrance channel complex is stabilized by approximately −2.4 kcal mol −1 and the exit channel complex is found to be approximately −1.2 kcal mol −1 below the separated products; the transition state is ∼1.6 kcal mol −1 above the two separated OH radicals, that is, ∼4 kcal mol −1 above the [OH...OH] complex.
Compared to the highest level of theory (CCSD(T)/augcc-pVTZ), the exothermicity of the reaction is overestimated at MP2 and underestimated at the DFT level, except for the MPW1K and MPWB1K functionals. The complexes on the reaction pathway are also best represented by the MPW1K and MPWB1K functionals. The energies given by the MP2 calculations, even with augmented basis sets, are too high. This suggests that higher excitations are needed to determine reliable activation energies. However, we expect that the geometry obtained at this level is correct; when we used the MP2 geometries, we found about the same activation energy in the CCSD(T)/aug-cc-pVTZ//MP2/aug-cc-pVDZ as in the fully optimized CCSD(T)/aug-cc-pVTZ treatment 1 .
Although the list of functionals we tested is far from being complete, the DFT approach to be used in forthcoming studies appears essentially limited to MPW1K and MPWB1K functionals (augmented with diffuse functions). We verified and confirmed this point in the following developments by comparing it to an increasing complexity of the reactive environment. It is worth mentioning that B3LYP, PBE1PBE, and MPWB1K without diffuse functions are inappropriate enough to lead to artificial hidden barriers (i.e., below the reactants).

Test evolution when H 2 O molecules are added
After reconsidering the gas-phase reaction to determine what level of theory is appropriate for the calculations on larger systems, we tested the potential energy surface of reaction (7) in the presence of an increasing number of H 2 O molecules at the CCSD(T) and DFT levels. We study the influence of the icy support here, therefore we use as starting point reference the A&A 638, A125 (2020) Notes. E e is the electronic energy, and E 0 is the electronic + zero-point vibrational energy. situation where one OH is attached to the support (from one H 2 O molecule to the full solid). In this context of reactivity on the surface, we consider as final point the complex in which the oxygen is still attached to the support.

Potential surface at CCSD(T) level
Starting with one H 2 O support, we obtained the reaction profile illustrated in Fig. 2 from the CCSD(T)/aug-cc-pVTZ//MP2/augcc-pVDZ energies reported in Taking a support of two H 2 O in the form of the dimer obtained above, we obtained the reaction profile illustrated in Fig. 3 from the CCSD(T)/aug-cc-pVTZ//MP2/aug-cc-pVDZ energies reported in Table 3. The starting point of the reaction was a triangular complex, [H 2 Odimer...OH + OH], in which one OH radical makes two hydrogen bonds with the two H 2 O molecules forming the water dimer. The reaction progresses by accreting the second OH radical in the form of a quadrangular complex that is better represented by the OH...OH complex found for the gas phase that interacts at both ends with the water -4.9 dimer. The next stationary point, that is, the transition state, follows from the latter complex with a structure similar to that of the gas-phase transition state (Sect. 3) that interacts with the supporting water dimer. It is found 1.2 kcal mol −1 below the starting point and 3.7 kcal mol −1 above the intermediate complex. We confirmed that this former triangular complex also correlates with another transition state that is much higher in energy (∼13.2 kcal mol −1 ) above the reference energy, which makes it a highly unfavorable pathway (not shown). The final product is a complex between the oxygen atom and the water trimer that is formed as before.
With a support of three H 2 O added in the form of the optimized trimer obtained above, we found a reaction profile different from those previously obtained, that is, the trend that the activation barrier is lowered while the number of water molecules increases is no longer true for n = 3. We found a real activation barrier of 1.7 kcal mol −1 on the way to a cubic structure that is the most stable form. This agrees with the result obtained by Wales & Hodges (1998) for n = 4 in a systematic study of the structure of the ice cluster (n = 2−21) H 2 O molecules. This structure does not reproduce even a partial geometry of solid ice, however. We therefore discontinued this type of modeling for ice, that is, optimized clusters. However, the indication given by the former very small clusters remains valid as a trend.

Potential surface at the DFT level
Knowing that CCSD(T) is not applicable to larger models, we reconsidered the above potential surfaces at the MPW1K and MPBW1K that provided the best approaches to the reference in the gas phase (Table 1).
With one H 2 O support (see Table 2), the MPW1K/ MPWB1K functionals provide a fair approach to the CCSD(T) A125, page 4 of 8 P. Redondo et al.: Reconstruction of water ice Table 3. Relative energies (kcal mol −1 ) for the reaction profile given in Fig. 3  Notes. E e is the electronic energy, and E 0 is the electronic + zero-point vibrational energy. exothermicity and activation energy, with a slight advantage for MPWB1K.
In the case of two H 2 O support (see Table 3), as in the studies we referred to above, the results obtained with MPW1K/ MPWB1K are similar to those of the CCSD(T) calculations for both the exothermicity and the activation energy. The trend is identical: the result for MPWB1K is slightly better.
The same drawbacks as in CCSD(T) appear when the H 2 O trimer is the starting complex that reacts with the OH radicals.

Modeling ice with small clusters
We found roughly the same energy profile with similar structures for each of the stationary points. The two approaches we used, CCSD(T) and DFT (especially MPWB1K) give similar values for the exothermicities and activation barriers.
The parameter that controls the facility for the OH + OH reaction to proceed is the activation energy. It decreases from its gas phase value with the addition of one or two molecules of water. However, the fact that this tendency, which is valid for CCSD(T) and DFT, is not followed when more H 2 O is added to the support is consistent with the results of Wales & Hodges (1998), who reported a systematic study of the structure of ice clusters with increasing dimension (up to 21 H 2 O molecules). From the H 2 O trimer on, none of them can be taken as a substructure of the most stable I h compact ice. We add that numerically extrapolating the interaction energy between one single H 2 O and a free-flyer molecule to an adsorption energy on ice can lead to erroneous conclusions, especially when more than one adsorption site is possible (Vakelam et al. 2017). Consequently, the approach based on the successive growth of simple water clusters can only give an indication of trends that has to be confirmed by a more elaborate model.

Reaction on water layers
To approach the ice situation while keeping a tractable model support, an H 2 O aggregate of reasonable dimension can be obtained by an appropriate cut of 24 H 2 O molecules (24-1L) assembled in the form of six fused cyclohexane-like structures that might represent a significant part of the surface bilayer of the crystalline ice (I h ). An even more realistic description can be obtained by superimposing two of these bilayers (48 H 2 O molecules) to model the role of the underlying bulk. Both clusters satisfy the Bernal & Fowler (1933) constraint.
With the MPWB1K functional we previously selected, we can study the reconstruction of the ice surface following the reaction path in additon to the ice model. This is not possible at the highly correlated level we used previously for the treatment of the small aggregates. We anticipated calculations on (48-2L) clusters and therefore employed two basis sets of different dimensions for the (24-1L) cluster: aug-cc-pVTZ as in the small cluster calculations, and 6-311++G(3df, 2p) to investigate alternative solutions of lower cost. Because the reaction profiles and stationary points are quite similar for the 24(1-L) and 48(2-L) ice supports, we present only one figure as illustration with snapshots to show the structural details of each critical point (Fig. 4).
(i) The first test was calculating OH + OH on the single bilayer (24-1L). We were aware of the problems of an eventual flexibility of the structure (as in the 3 H 2 O case) and therefore kept the bilayer frozen in order to keep an ice-like structure. Only the reactants were optimized along the reaction path.
The reaction starts from an OH radical on the ice surface in the form of a bridged complex in which OH is attached to the 1, 3 positions of the hexagonal surface in a chair conformation and was taken as energy reference for the process we considered. The next stationary point is also a bridged complex in which the incoming OH forms an H 2 O dimer with the OH that was previously adsorbed on the surface. This complex is stabilized by ∼12 kcal mol −1 below the energy of the starting products. The next step on the reaction path is a transition state that is ∼10 kcal mol −1 above the bridge complex, but still ∼2 kcal mol −1 below the reactants energy. This confirms the trend we previously found with the clusters, that is, the presence of a hidden barrier. This means that there is no activation energy to oppose the continuation of the reaction. The final product is the complex composed of the original bilayer increased by one H 2 O molecule with the oxygen atom attached to the surface. In this position, the molecule is well placed to react with any diffusing hydrogen atom to form a new OH radical and reinitiate the whole mechanism. The two basis sets lead to practically the same values within the uncertainty margin (Table 4). This justifies the use of 6-311++G(3df, 2p) for the next larger calculations.
(ii) A second calculation was performed with the two-bilayer structure (48-2L). The subsurface frozen bilayer is intended to maintain the structure that would be imposed by the bulk of real ice. By imposing this constraint, we leave the possibility for relaxation in the upper bilayer by optimizing the six nearest H 2 O molecules that form the cyclohexane-like environment that supports the reaction without distorting the whole structure dramatically. The reaction profile of OH + OH supported by a two-bilayer of water is very similar to the profile that is obtained with one bilayer. The energy values reported in Fig. 4 are those where the ice support is optimized.
The reaction starts with an incoming second OH radical that leads to the formation of a bridged complex (∼7.5 kcal mol −1 below the reference energy) in which the OH...OH dimer is  Notes. E e is the lectronic energy, and E 0 is the lectronic + zero-point vibrational energy. chair conformation. The next stationary point is a transition state whose structure is shown in Fig. 4, situated between 1 and 2 kcal mol −1 below the reactants level. It shows again that we are in presence of a hidden barrier on the reaction path and that there is no activation energy to oppose the continuation of the reaction. As with one bilayer, the final product is a structure composed of the original layer augmented by a H 2 O molecule to which the oxygen atom is attached at the surface. This system can either react with any diffusing hydrogen atom to form a new OH radical and then reinitiate the reconstruction mechanism of the ice mantle, or it can detach from the ice.

Reaction on solid ice
We finally configured the icy environment with a solid surface of infinite dimension. In this periodic approach, the model unit cell that is cut into apolar hexagonal ice I h contains two H 2 O bilayers. It is an elementary volume whose replication in three spatial directions creates a series of ice surfaces of infinite dimensions. This volume is defined by three structural parameters a, b, and c. Here, the basal plane dimensions (corresponding to the surface on which the reaction proceeds) are a = 7.145 Å and b = 8.73 Å. The third dimension of the cell was determined so that there is no interaction between molecules generated by a translation perpendicular to the ice surface. The same reaction path as for clusters (see Fig. 4) was considered, starting from one OH radical adsorbed on the ice model surface of infinite dimensions, taken as reference. As in the cluster case, the second OH is a free flyer, and we follow the bridge conformation throughout, as for the 48(2-L) simpler model. The trend for stabilization shown by the interaction energy of the transition state with the ice support when n(H 2 O) increases is even emphasized. It shows a hidden barrier of −3.8 kcal mol −1 below the reactants, that is, deeper than that found in the cluster case. This confirms the hypothesis of a possible ice reconstruction through this mechanism.   Snapshots of the energy profile of the solid-state reaction mechanism are given in Fig. 5. The critical energies are reported in Table 5 together with the results of the preceding calculations for comparison.

Summary and concluding remarks
The irradiation of the interstellar ices is known to break the water molecules into three fragments: atomic hydrogens, oxygens, and OH radicals. The main goal of this work was to increase the chemical information on the reconstruction mechanisms of ice mantles of the interstellar grains at the molecular level. The direct recombination of H and OH fragments (reaction (1)) was originally proposed for the formation of the ice. However, the hydrogen atoms produced by the diffusing irradiation easily give mostly molecular H 2 . These hydrogen molecules could also react with the OH radicals (reaction (2)), and it has been shown that this reaction could proceed through the tunneling effect (Oba et al. 2012;Meisner et al. 2017).
The reconstruction of the ice under the hypothesis that the OH radicals produced by destructive irradiations might recombine through reaction (7) at the top of the remaining ice was considered. Two types of numerical simulations based on conceptually different approaches were employed. The first approach relies on individualizing the ice support down to the molecular level using high-level ab initio calculations (CCSD(T) and recently developed DFT functionals. The second approach relies on globalizing the ice support as a solid of infinite dimensions using a first-principle periodic calculation. The two approaches give a similar qualitative answer: the results for large H 2 O molecular clusters are similar to the results we obtain for the solid surface that has infinite dimensions.
We found that in both approaches, no activation energy prevents the pursuit of the H 2 O synthesis at the surface of water ice. In the present case, where we constrained the top bilayer to respect the structure of the underlying bulk, the new H 2 O molecule that formed is well positioned to help in the reconstruction of the next bilayer. The oxygen atom released in the process is also accessible to diffusing H radicals or to other H atoms that are present in the environment. The first-principle periodic calculation carried out finally agreed with the large cluster treatments: the reconstruction of the ice mantles by reaction OH + OH → H 2 O + O is a barrierless process when it is catalyzed by ice support.