O₂ formation in cold environments

Abstract

The diffusional dynamics of atomic oxygen in and on amorphous solid water (ASW) to form molecular oxygen is characterized. Reactive molecular dynamics simulations to study bond breaking and bond formation show that vibrational relaxation of the highly excited diatomic occurs on the 10 ns to 100 ns time scale. The relaxation process is highly nonexponential and can be characterized by a stretched exponential decay reminiscent of the dynamics of glasses. The stretched exponents range from β = 0.15 for relaxation on the surface to β = 0.21 for the dynamics in bulk. It is also found that coupling of the O₂ relaxation to the internal water modes occurs which speeds up the vibrational relaxation by a factor of 4. Extrapolation of the stretched exponential decay to 1 μs yields a final vibrational quantum number v = 2 for O₂(X³Σ⁻_g), consistent with experimental results from photolysis of SO₂ on ASW at 193 nm which find v ≤ 3. Desorption energies of water from the surface range from 1.5 to 2.0 kcal mol⁻¹ compared with 1.8 kcal mol⁻¹ found from experiment, depending on whether the water molecules are flexible or not.

---

1 Introduction

The mobility of atomic and molecular species on solid surfaces is a primary driver for surface-induced chemistry at low temperatures. It has been found from experiments^1,2 and atomistic simulations^3,4 that H and O atoms can diffuse on amorphous water surfaces down to low temperatures (∼10 K). This is relevant for the possibility to form diatomic or larger molecules (such as O₂, O₃ or CO₂) in laboratory experiments² and may also play a role for the conditions in interstellar space.^5–7 If other atomic species than hydrogen (such as N, C, or O) are able to diffuse on cold surfaces, the range of molecules that can be formed by diffusion and aggregation is considerably expanded. Given the very long processing times (∼years or longer) that are available in translucent, diffuse and dense molecular clouds,^5,8 the accessible chemical space from such reaction pathways may be affected appreciably.

Molecular oxygen is an important tracer for characterizing the physical properties of star forming regions and molecular clouds and is implicated in various chemical processes leading to diverse chemical compositions in the interstellar space. Molecular oxygen in interstellar clouds was directly observed in the ρOphiuchi cloud using submillimeter wave spectroscopy.⁹ Understanding of the formation and ensuing concentration¹⁰ of O₂ is relevant for better characterizing the evolution of stellar clouds and their chemical composition.

In interstellar clouds, dust provides a suitable substrate for deposition and chemical synthesis of a variety of small molecules in a bottom-up fashion. Due to the low temperatures (10 K to 50 K),^5,8 residence times are long despite the small barriers (∼1 to 2 kcal mol⁻¹)^11,12 for desorption and migration which also allows reactions to occur. In laboratory experiments with amorphous solid water as the substrate, deposition of oxygen (O) and O₂ leads to ozone formation at 10 K.^13,14 These experiments suggest that using pure oxygen species, only two reactions are relevant: O + O → O₂ and O + O₂ → O₃.² In the experiments atomic oxygen (O(³P)) was generated from dissociating O₂ through a microwave discharge and it was verified that all participating species are in their ground electronic state. In modeling interstellar processes, the O + O → O₂ channel is sometimes included⁵ and sometimes not.¹⁵ Another, and probably the predominant,^16,17 O₂-formation reaction is O + OH → O₂ + H. A more exhaustive list of O₂-formation reactions, including O + OH → O₂ + H in the gas phase¹⁸ or on the surface of grains,¹⁹ can be found in the KIDA database.²⁰

It has been previously shown by experiment and simulations that migration of atomic species, such as oxygen, in ASW is readily possible on the sub-microsecond time scale over distances of several ten Angstroms.^2,4 Similarly, oxygen atom diffusion in Neon matrices at temperatures as low as 4 K has recently been established.²¹ Diffusing oxygen atoms have also been implicated in the formation of CO₂ through O + CO → CO₂.^22–24 This pathway for the formation of molecular species from adsorbates is generally known as the Langmuir–Hinshelwood mechanism.²⁵

In laboratory experiments using amorphous solid water, formation of ground state ³Σ⁻_g (from OH + O(³P) with ΔH = −68 kJ mol⁻¹) and excited state ¹Δ_g O₂ (from OH + O(³P) with ΔH = 26 kJ mol⁻¹) was found following electronic excitation of the ASW with 157 nm irradiation at 90 K.¹⁷ This process for O₂ formation follows a non-equilibrium route, involves higher electronic states and differs from the processes considered in temperature-programmed desorption experiments.² The study of SO₂ photolysis at 193 nm adsorbed on ASW yields O₂(³Σ⁻_g) from recombination of O(³P) + O(³P) with ΔH = −498 kJ mol⁻¹.¹⁷ Furthermore, a model for formation of O₂ following electronic excitation of ice by fast ions, low-energy electrons, and UV photons at low temperatures was presented.²⁶ In this model trapped O(³P) is the principal species considered for formation of O₂ at low doses. This route for O₂ formation is also implicated in the atmospheres of the icy moons of Jupiter and in the rings of Saturn, where molecular oxygen has been detected as a result of the decomposition of ice by energetic ions and electrons or by UV photon-induced production from ice.

In the present work the dynamics of O₂ formation in and on amorphous solid water (ASW) from collision of two O(³P) is studied. The main questions considered are the possibility for formation and stabilization, the likely formation time scales, the vibrational relaxation of highly excited O₂ following recombination, the desorption of molecular oxygen following a Langmuir–Hinshelwood mechanism and the energy transfer between the energized molecule and the surrounding water matrix. After presenting the methods, atomic oxygen recombination within and on the surface of an ice grain is considered and analyzed. Then, the desorption dynamics from the oxygen is considered.

2 Methods

2.1 Molecular dynamics simulations

All reactive MD simulations are carried out using CHARMM²⁷ with MS-ARMD^28–30 and the Leap Frog integrator with a time step of Δt = 0.1 fs. Such a short time step is required due to the high vibrational excitation upon oxygen–oxygen recombination and also in simulations with flexible water (vide infra). The simulation system is a cubic box with dimensions (31 Å)³ with periodic boundary conditions and the nonbonded cutoff is 14 Å. The box size along the z-axis is elongated by 20 Å to allow for surface desorption of the O and O₂ adsorbate. Initial velocities are drawn from a Maxwell–Boltzmann distribution at 50 K. All simulations are initialized as follows: a minimization of 100 steps using the steepest descent method is followed by 10 ps heating, where the temperature is increased by 10 K every ps and then relaxed for the last 5 ps, and 50 ps equilibration dynamics with rescaling of the temperatures every 2.5 ps. Then, the following different types of simulations are considered: the (rigid) TIP3P³¹ or the (flexible) KKY^32,33 water model are combined with either a reproducing kernel Hilbert space (RKHS) or a Morse representation of the reactive O–O potential (vide infra) to follow the recombination and relaxation dynamics.

Oxygen (O₂) formation from collision of two oxygen atoms is studied by using reactive MD which requires two potential energy surfaces, one for the unbound state (O + O), and a second one for bound O₂. MS-ARMD then handles the transition between the two states.²⁸ Only the state with the lowest energy is propagated in the simulations except for configurations for which the energy difference between the lowest and next higher-energy states is within ΔV for which mixing of the two (or multiple) states occurs. In those cases the total energy for each configuration x is the weighted sum of the involved states:

where V_MRMD is the total energy surface, w_i(x) is the weight of the ith state with energy V_i. The weightis calculated by re-normalizing the raw weights w_i,0, using an exponential decay function of the energy difference between the minimum v_min and the ith surface.In this work two states are used for the dynamics: the unbound state, described by two isolated oxygen atoms which interact through nonbonded interactions with a charge of zero and van der Waals parameters σ = 1.7 Å and ε = 0.12 kcal mol⁻¹, as parametrized in the CHARMM36 force field³⁴ and the bound O₂ in its ground electronic state (³Σ⁻_g) for which accurate MRCISD + Q calculations are available.³⁵ These calculations employed the cc-pCVQZ basis set and the 12 electrons in the last 8 orbitals, including core-valence and relativistic contributions from the lower electronic states, were correlated.³⁵ The energies are represented using a reproducing kernel Hilbert space (RKHS) approach^36,37 or they were fitted to a Morse potential V(r) = D_e(1 − e^{−β(r−r_e)})². The parameters from a least-squares fit are D_e = 123.1 kcal mol⁻¹, r_eq = 1.2075 Å, and β = 2.68673 Å⁻¹ for the well depth, equilibrium separation and steepness parameters, respectively. For the simulations using a RKHS representation the MRMD module in CHARMM was suitably modified.

For water two models are considered in the present work. One is the standard transferable interaction potential 3 points (TIP3P) water model³¹ which is rigid. In order to study energy transfer between the highly excited oxygen molecule and the internal water degrees of freedom a flexible water model, a modified Kumagai, Kawamura and Yokokawa (KKY) model was used as well.^32,33 This model has also been used for studies of ASW in the past.^3,33 The functional form of the KKY model involves a stretching term

E_str = D_e(1 − e^{−β(r−r_eq)}) − D_e (4)

with D_e = 75.037 kcal mol⁻¹, r_eq = 0.957 Å, and β = 2.74 Å⁻¹. The bending potential iswithwhere θ_eq = 99.50°, r_i is one of the O–H bond distances and g_r = 6.99 Å⁻¹ and r_m = 1.402 Å are original force field parameters and f_k = 0.11 J is a re-parametrized bending parameter, fitted to match the experimental bending frequency (see Table S1, ESI†).³³

The ASW is prepared by annealing liquid water from 300 K to 50 K for a total simulation time of 10 ns. The final density (1.01 g cm⁻³) of the system is compatible with what is expected for high density amorphous solid water.³⁸ To study O₂ formation within ASW, a cavity with volume ∼1000 ų is created by removing 15 water molecules. Then, the system is minimized and equilibrated and further relaxed by running a 1 ns equilibrium simulation at 50 K. This leads to one main cavity with a volume of 583 ų and several smaller cavities with volumes ranging from 30 to 100 ų. The cavity volumes were determined from GHECOM.³⁹ No further deformation is observed after this preparation.

3 Results

In the following, the results are discussed separately for formation of O₂ in an inside ASW-cavity and for recombination at the surface.

3.1 Recombination in the cavity

Initial conditions for O₂-formation were generated from a 500 ps equilibrium MD simulation at 50 K with the two oxygen atoms in the same internal cavity as illustrated in Fig. 1. From this simulation, 1000 structures were chosen to initiate MS-ARMD simulations. The initial distribution of oxygen–oxygen separations (see Fig. 2) has distinct maxima at 3.4 Å, 7.0 Å, and 10.0 Å covering a wide range of initial conditions. The oxygen atoms are almost thermalized, with a kinetic temperature of 70 K compared with 50 K of the surrounding water lattice, although the meaning of temperature for a diatomic molecule is somewhat debatable.

Fig. 1 Two physical conditions in which O₂ formation is studied: panel A shows two oxygen atoms (red) trapped inside a cavity of the ASW (grey); in panel B shows the two on top of the water surface.

Fig. 2 Initial distribution of the oxygen–oxygen separation P(r) at the beginning of the rebinding simulations using [TIP3P/RKHS]. P(r) has three maxima at ∼3.4 Å (the minimum of the O–O van der Waals complex), at 7.0 Å and at 10.0 Å. The total distribution (black) is compared with the one for initial conditions that lead to rebinding (red, 793 events on the 100 ps time scale). The inset shows the time series for the oxygen–oxygen distance r(t) for one sample trajectory and confirms that all distances in P(r) can be sampled even for one single trajectory.

All 1000 rebinding simulations are run for 100 ps and 16 of them are extended out to 5 ns to better characterize the relaxation process. The fraction of O₂-forming simulations is 79.3% within 100 ps with half of them forming O₂ within 20 ps. The whole range of initial separations (between 3.4 Å and ∼13 Å) leads to formation of O₂, see red distribution in Fig. 2. Typical O₂-formation curves (black, red, blue, indigo) are reported in Fig. 3 whereas for the green trace no molecular oxygen is formed on the 100 ps time scale. The red trace shows that O₂ can also be formed transiently (at 40 ps), break up, reform and stabilize at a later stage. The overall distribution of the O₂-formation times is strongly peaked at short times with an extended tail, see top panel of Fig. 3.

Fig. 3 Top: Histogram P(τ) of the formation times for the 793 reactive simulations; 40% react within the first 20 ps and the cumulative fraction is shown as the red line. Several time series r_OO(t) are shown in the inset. In black and indigo two examples of a reaction within the first 20 ps, a late reaction (around 40 ps, blue), and a simulation where dissociation occurs soon after O₂ recombination (red) and in green a simulation where the atoms remain in the unbound state. Bottom: Interatomic distance (panel A), O₂ kinetic energy (panel B) and temperature profiles for water (red trace) and the overall system (black trace, panel C) for one recombination process. At t = 23 ps a nonproductive collision occurs and stabilization of O₂ starts after t > 55 ps. In panel B the right hand side y-axis reports the temperature associated with the kinetic energy of O₂.

Because rebinding of the two oxygen atoms leads to the liberation of the O₂ binding energy, it is expected that the kinetic temperature in the system increases. This is illustrated for one particular trajectory in Fig. 3 bottom panel. The lower panel in Fig. 3A reports the O₂ distance as a function of simulation time for one specific trajectory. At t = 23 ps a short encounter of the two oxygen atoms occurs which, however, does not stabilize O₂. Only after 55 ps molecular O₂ forms. Concomitantly, the kinetic energy and the temperature increase significantly, as is shown in Fig. 3B. This also leads to a marked temperature increase in the overall system from 50 K to >60 K.

Due to the weak coupling between the O₂ stretch and the phonons of the solid water matrix, vibrational relaxation of the diatomic is slow. In order to quantify this, 16 trajectories using [TIP3P/RKHS] were continued out to ∼5 ns, see Fig. S2A (ESI†). Relaxation times were estimated by fitting the time dependent oxygen–oxygen separation r_OO(t) to a sum of exponentials

or to a stretched exponential

r_OO(t) = a exp(−(t/τ)^β) + r₀ (8)

with r₀ = 1.2075 Å, which is the O₂ equilibrium separation. Using stretched exponentials is probably more meaningful in the present context as such a dependence has been found to be particularly suitable to describe the relaxation of glasses⁴⁰ and for ligand rebinding to proteins⁴¹ both of which are characterized by distributed barriers. In the present case (see Table 1) it is also found that the quality of the fit improves when one stretched exponential function is used compared to a fit to two exponentials. Considering one particular trajectory, the exponential decays involve a rapid decay τ₁ = 0.5 ns and a longer relaxation τ₂ = 14.6 ns whereas one stretched exponential leads to τ^SE = 1.4 ns with β = 0.21. The average long-time relaxation for the 16 trajectories is τ₂ = 15 ± 0.5 ns, see Fig. S2A (ESI†).

Table 1 Parameters and quality for fitting the O₂ vibrational relaxation to exponential (first line), one stretched exponential (second line) using eqn (7) and (8), respectively. The force fields (O₂ potential and water) with which the simulations were run are also indicated

Region	FF	a 1 (Å)	a 2 (Å)	τ 1 (ns)	τ 2 (ns)	β	RMSD (Å)
Inside	[TIP3P/RKHS]	0.28	0.38	0.5	14.6	—	3.638 × 10^−3
		1.38	—	1.4	—	0.207	6.546 × 10^−5
Surface	[TIP3P/RKHS]	0.35	0.31	1.1	83.4	—	9.154 × 10^−3
		1.08	—	1.1	—	0.144	2.032 × 10^−4
Surface	[KKY/RKHS]	0.27	0.34	0.7	20.7	—	2.156 × 10^−3
		1.12	—	0.9	—	0.182	3.082 × 10^−4
Surface	[TIP3P/Morse]	0.29	0.25	3.2	206.9	—	3.223 × 10^−2
		1.08	—	1.3	—	0.139	4.680 × 10^−4

---

Upon O₂ formation the binding energy released will be partially converted into heat and dissipates through the ASW. To quantify this, the averaged kinetic energy of all water molecules projected along the z-direction is averaged on a (x, y)-grid with voxels 1 × 1 × 30 ų (i.e. a column of height 30 Å and (x, y)-base 1 × 1 Ų). The average kinetic energy of all water atoms in a single column is calculated, and averaged over 8 independent trajectories for times 10 ps before the recombination (at t = 0), i.e. for −10 < t < 0. Recombination was considered to have occurred for r_OO < 2.5 Å. The corresponding heat map is shown in Fig. 4A. The low energy region (blue) in the center of the chart corresponds to the cavity region where fewer water molecules are located along the z-direction.

Fig. 4 Averaged temperature distribution in the (x, y)-plane averaged over 8 simulations using the RKHS PES and TIP3P water before (panel A: −10 ≤ t ≤ 0 ps) and after (panel B: 0 ≤ t ≤ 10 ps) O₂-formation. Here, t = 0 is defined as the instance along each of the trajectories for which r_OO < 2.5 Å for the first time and all time series are aligned with respect to this reference. The color scale reports the kinetic energy (in kcal mol⁻¹). The averaged position of the center of mass of O₂ is indicated by filled circles. The low-energy region in the center of the simulation box is due to the cavity which contain fewer water molecules along the z-axis that can contribute to the kinetic energy. Panel C shows the difference between panels A and B which highlights the direction and magnitude of energy transfer. The region “a”, delimited by the continuous line ellipse, represents the region where the recombination takes places, while the second region “b” is a qualitative delimitation of the energy flow after the reaction. Panel D shows the energy distribution for two not reacting oxygen atoms.

A similar heat map was generated for the time 0 < t < 10 ps, i.e. after recombination, see Fig. 4B and the difference map is shown in Fig. 4C. Heating of the ASW molecules around the recombination region is found to occur in a range of ∼8 Å on the 10 ps time scale. Furthermore, energy flow is anisotropic due to the amorphous nature of the water matrix. Fig. 4D shows an ensemble average over 10 ps and 8 trajectories for a situation in which the two oxygen atoms are separated by an initial distance of 5 Å and are not allowed to recombine, i.e. MS-ARMD was disabled and the two oxygen atoms only interact through van der Waals interactions. The distribution of the kinetic energy in Fig. 4A and D are comparable and differ substantially from panel B.

3.2 O₂ formation: surface reaction

Atomic oxygen recombination on the surface of ASW is probably even more relevant under astrophysical conditions. There are two typical scenarios: the Langmuir–Hinshelwood (LH)^42,43 mechanism whereby the two reacting species are in thermal equilibrium with the environment before they react or the Eley–Rideal (ER)⁴⁴ mechanism for which a thermalized, surface-adsorbed reactant collides with an impinging, non-thermalized partner. The scenario explored here is the LH mechanism.

For these simulations, one of the oxygen atoms is initially localized (but unconstrained) in a surface cavity sufficiently deep to avoid spontaneous desorption for the present simulations and the second oxygen atom diffuses freely on the surface. Altogether, one 110 ns [Morse/TIP3], one 25 ns [KKY/RKHS], and one 10.5 ns [RKHS/TIP3] simulation were initially run and analyzed. The relaxation of the O–O bond length for the [Morse/TIP3] simulation is shown in Fig. 5 and Fig. S2, S3 (ESI†). Relaxation is characterized by a rapid initial phase during which ∼50% or the binding energy relax (relaxation to v = 15) within a few nanoseconds. This is followed by a second, considerably slower phase extending over many 10 s to 100 s of nanoseconds during which collisions with the surrounding water lattice can lead to de-excitation and re-excitation. Similarly, the [TIP3P/RKHS] simulation (see Fig. S2, ESI†) shows a first phase (with time constant τ^DE₁ = 1.1 ns from a double exponential (DE) fit) during which the average O–O distance decreases from 〈r_OO〉 = 2.1 Å to 1.5 Å within 3.6 ns, followed by a much slower relaxation (τ^DE₂ = 83.4 ns) during which the O–O separation only decreases by another 0.1 Å. Using a single stretched exponential yields τ^SE₁ = 1.1 ns and β₁ = 0.144, see Table 1 and for two stretched exponentials the parameters are β₁ = 0.51 for the short-time process (τ^SE₂ = 0.9 ns) and β₂ = 0.02 for the long-time process (τ^SE₂ = 1.2 ns). At the end of this simulation the water temperature has increased by ∼15 K, see Fig. 6.

Fig. 5 Time evolution of the O₂ vibrational relaxation and its fit to a stretched exponential function (eqn (8)) for a [TIP3P/Morse] simulation with recombination on the ASW surface. Scattering of the O₂ from the surrounding water lattice leads to the peaks in the relaxation curve. The right-hand y-axis label reports the vibrational quantum number from solving the 1d-Schrödinger equation using LEVEL.⁴⁵ After 110 ns of simulation time the O₂ molecule is still in a state corresponding to v = 15.

Fig. 6 The kinetic temperature of the O₂ (top panel) and the water (bottom panel, red) and the total system (green) after formation at t = 0 from a [TIP3P/RKHS] surface simulation. After 20 ns the average temperature of the system has increased by ∼15 K. The spikes in the black trace are due to collisions of O₂ with the surrounding water matrix.

Simulations with the Morse Potential for the O–O coordinate and the TIP3P water model establishes that relaxation is slow. This model is computationally less demanding than the RKHS representation and longer simulations are possible. The relaxation times from the 110 ns simulation are τ₁ = 3.2 ns and τ₂ = 207 ns which can only be an estimate because the data from which this value was determined does not probe dynamics on this time scale. Nevertheless, these results confirm that coupling between the oxygen stretching mode and the water phonon modes is small.

To evaluate how water flexibility influences the relaxation dynamics of molecular oxygen, a 10.5 ns [KKY/RKHS] simulation was carried out. Such simulations are even more computational demanding because the KKY water is flexible and a time step of Δt = 0.1 fs is required in the MD simulation for energy conservation. Vibrational relaxation times decreases by a factor of two to four, with τ₁ = 0.7 ns and τ₂ = 20.7 ns for a double exponential fit. For the stretched exponential fits, the values of β increase by ∼30% compared to simulations with rigid waters (see Table 1) and the two relaxation times decrease to τ^SE₁ = 0.5 ns and τ^SE₂ = 1.0 ns. This suggests that the O₂ relaxation dynamics is somewhat less nonexponential probably due to coupling to the water matrix which also speeds up the relaxation process.

Stretched exponential decays have been employed in the past to discuss relaxation within a diffusion-trap model and to characterize the dynamics in disordered molecular systems.⁴⁶ The case β = 1 in eqn (8) corresponds to a conventional exponential decay and for all other cases 0 < β < 1. The value of β can be related to the effective dimensionality of the problem through β = d*/(d* + 2) where d* = f × d, d is the dimensionality of the problem (here d = 3) and f is the fraction of active channels for relaxation.⁴⁶ Discussing the present results (see Table 1) for one stretched exponential, it is found that for surface relaxation β ∼ 0.15 whereas for relaxation within the bulk β ∼ 0.2. For relaxation in the bulk and d = 3 this yields f = 1/6 whereas for relaxation on the surface the value for d is not obvious a priori. For d = 3 the number of relaxation channels is f = 0.12 whereas for d = 2 the value increases to f = 0.18. It is reasonable to expect that more rapid relaxation is caused by a larger number of relaxation channels (i.e. a larger number of collision partners). Because τ^bulk > τ^surface for [TIP3P/RKHS] one finds d = 2 for the relaxation on the surface.

In order to better characterize the formation probability for the surface reaction an additional 534 simulations, using [TIP3/RKHS], were performed. The two oxygen atoms were placed in two neighboring cavities on the surface with an initial distance of 8 Å and velocities randomly chosen from a Maxwell–Boltzmann distribution at T = 50 K. Within 25 ps, 25% of the simulations lead to recombination with an average formation time of 〈t_formation〉 = 10.8 ps, see Fig. 7 This establishes that not only a specifically prepared initial state leads to O₂ formation but also initial conditions which require diffusion of the two oxygen atoms before recombination can take place.

Fig. 7 Top: Time series of the O–O distance during 180 ps for eight different simulations, with one oxygen in a deep crevice: for three trajectories no recombination is observed, because the two atoms diffuse away from the cavity where the first atom is trapped; five reactions instead leads to the molecule formation. The colored arrows indicate the initial distance, r_OO, of the two atoms. In the central panel, the moving average of those last simulation is shown: as from the previous simulations, a slow relaxation of the O–O distance is observed. The bottom panel shows the maximum height of the ice surface, i.e. the maximum z-coordinate of any water molecule, compared with the z-coordinate of the center of mass of O₂. The molecules remain, in all simulations, on the surface, regardless of whether they rebind or not.

3.3 Desorption of O₂

Despite energy in excess of 100 kcal mol⁻¹ is available upon oxygen recombination, which compares with desorption energies of 2 to 3 kcal mol⁻¹,^4,47 O₂ remains localized on the surface for the simulations described so far. In order to explicitly probe desorption, a set of 10 simulations, 50 ps each, is run with velocities from a Maxwell–Boltzmann distribution at 50 K but with the z-component (perpendicular to the ASW surface) scaled according to

v_z = v_z0(1 + λa) (9)

with λ ∈ [0,…,50] in increments of 1, and a = 0.1. With increasing λ the velocity orthogonal to the surface increases and desorption becomes more probable. For each value of λ the average probability for desorption was determined.

For a rigid water model, desorption is observed for an average perpendicular kinetic energy of 2.0 kcal mol⁻¹ (see Fig. 8 black curve). If desorption occurs it takes place within 0.5 ps. Similar results are obtained when using a flexible water model but the energy for desorption is reduced to 1.5 kcal mol⁻¹ (see Fig. 8) because the O₂ stretch can now couple to the water deformation modes. Such energies for desorption are consistent with those determined from laboratory experiments of 1.8 kcal mol⁻¹.⁴⁷

Fig. 8 Probability for desorption from an average over 10 independent simulations with different z components perpendicular to the ASW surface. Black trace: simulations with rigid (TIP3P) water; red trace: simulations with the flexible KKY model. Desorption occurs for a vertical component of the kinetic energy of ∼1.7 kcal mol⁻¹ and for ∼1.2 kcal mol⁻¹, for TIP3P and KKY, respectively. The process is illustrated in the schematic.

4 Discussion and conclusion

In this work formation of molecular oxygen on cold ASW surfaces is examined both inside amorphous solid water and on its surface. Oxygen (³P) migration in and on ASW has been found from simulations⁴ and in experiments using neon matrices.²¹ If two oxygen atoms O(³P) approach one another sufficiently closely the process is rapid as it is barrierless. Once O₂ has formed, relaxation of the highly excited molecule takes place on the 10 to 100 ns time scale with full relaxation requiring considerably longer, though, well beyond the current simulation times. The rate limiting step for vibrational relaxation is coupling to the water phonon modes. For O₂ recombination on ASW it is found that even the highly vibrationally excited molecule does usually not desorb which also occurs due to inefficient coupling of the O₂ stretch to the water-phonon modes. With a flexible water model the desorption energy decreases by 25% compared to rigid water and vibrational energy relaxation speeds up by a factor of ∼4.

Experimentally, vibrational energy transfer has been found to be very inefficient for NO near the surface of an insulator. For NO(v = 1) the survival probability of scattering from LiF is 70% to 100%, with error bars up to 100% at all energies.⁴⁸ More recently, high survival probabilities (small probabilities for vibrational relaxation) have also been found for NO(v = 12) on LiF(001).⁴⁹ Similarly, the probability for vibrational excitation of NO scattered from diamond was found to be only 5 × 10⁻³ and the vibrational relaxation probability of CO(v = 1) on LiF is of comparable magnitude (10⁻³).⁵⁰ Finally, on NaCl(100), CO vibrational energy relaxation occurs on the ms time scale.⁵¹ Relaxation on insulators (mainly by phonons) is inefficient due to the large mismatch in the frequencies of the diatomic and the phonons which leads to small coupling between these degrees of freedom. Conversely, for a conductive solid, coupling of intramolecular vibrations with phonons and electron–hole pairs accelerates energy transfer such as for vibrational relaxation of CO on Cu(111) which occurs on the ps timescale.⁵²

Photolysis of adsorbed SO₂ onto ASW and recombination of O(³P), which is the process considered here, leads exclusively to O₂(³Σ⁻_g), predominantly in vibrationally excited states v ≤ 3. The time of flight in these experiments is on the order of μs which agrees quite well with results from the present simulations. Extrapolating the stretched exponential decay from the [KKY/RKHS] simulations from recombination on the ASW surface (see Table 1) to 1 μs yields r_OO = 1.2386 Å which corresponds to v = 2, consistent with experiment.¹⁷ Formation of electronically excited O₂(¹Δ_g and ¹Σ⁺_g) has been found from OH + O(³P) involving a barrierless and an activated process, respectively.¹⁷ The fraction of O₂ in its ¹Δ_g state was 2.5%.⁵³ Given that the ³Σ⁻_g, ¹Δ_g, and ¹Σ⁺_g of O₂ dissociate to the same asymptote (O(³P) + O(³P)) it will be interesting to include electronically excited states in future work. For these states PESs of similar quality are available.³⁵ The recombination dynamics can be treated with trajectory surface hopping (TSH)⁵⁴ including nonadiabatic transitions within the Landau–Zener^55,56 formalism. This has recently been done using multidimensional RKHS PESs for the reactive dynamics of the [CNO]-system.⁵⁷

The long vibrational relaxation times (ns to ∼μs or probably even longer) found here for O₂ on ASW are consistent with vibrational relaxation of a high-frequency diatomic on an insulating surface such as water. If the water molecules lack internal degrees of freedom (rigid, TIP3P) the relaxation times are longer by about a factor of 4 compared with simulations in which the water-bending mode is available for coupling (flexible, KKY). For more quantitative results on the vibrational relaxation times much longer MD simulations will be required, though.

The present findings suggest that diffusion of atomic species followed by collisions in and on ASW and subsequent vibrational relaxation can lead to stable molecules at low temperatures. This complements and is in agreement with laboratory experiments that find generation of O₂ and O₃.^{2,13,14,17,21,58–60} Such a mechanism for O₂ formation may also be operative at conditions prevalent in the interstellar medium although there, alternative routes such as the reactive O + OH collision are considered to be important, too.^16,17

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

This work was supported by the Swiss National Science Foundation through grants 200021-117810, and the NCCR MUST.

References

1. E. Matar, E. Congiu, F. Dulieu, A. Momeni and J. L. Lemaire, Mobility of D atoms on porous amorphous water ice surfaces under interstellar conditions, Astron. Astrophys., 2008, 492, L17–L20.
2. M. Minissale, E. Congiu, S. Baouche, H. Chaabouni, A. Moudens, F. Dulieu, M. Accolla, S. Cazaux, G. Manico and V. Pirronello, Quantum Tunneling of Oxygen Atoms on Very Cold Surfaces, Phys. Rev. Lett., 2013, 111, 053201.
3. M. W. Lee and M. Meuwly, Diffusion of atomic oxygen relevant to water formation in amorphous interstellar ices, Faraday Discuss., 2014, 168, 205–222.
4. M. Pezzella, O. T. Unke and M. Meuwly, Molecular Oxygen Formation in Interstellar Ices Does Not Require Tunneling, J. Phys. Chem. Lett., 2018, 9, 1822–1826.
5. A. G. G. M. Tielens and W. Hagen, Model Calculations of the Molecular Composition of Interstellar Grain Mantles, Astron. Astrophys., 1982, 114, 245–260.
6. D. Ruffle and E. Herbst, New models of interstellar gas-grain chemistry – III. Solid CO₂, Mon. Not. R. Astron. Soc., 2001, 324, 1054–1062.
7. S. B. Charnley and S. D. Rodgers, in Astrochemistry: Recent Successes and Current Challenges, ed. D. C. Lis, G. A. Blake, E. Herbst, International Astronomical Union, 2006, pp. 237–246.
8. H. M. Cuppen and E. Herbst, Simulation of the formation and morphology of ice mantles on interstellar grains, Astrophys. J., 2007, 668, 294–309.
9. B. Larsson, et al., Molecular oxygen in the ρ Ophiuchi cloud, Astron. Astrophys., 2007, 466, 999–1003.
10. R. Liseau, et al., Multi-line detection of O₂ toward ρ Opiuchi A, Astron. Astrophys., 2012, 541, A73.
11. J. He, K. Acharyya and G. Vidali, Binding Energy of Molecules on Water Ice: Laboratory Measurements and Modeling, Astrophys. J., 2016, 825, 89.
12. J. He, S. Emtiaz and G. Vidali, Measurements of Diffusion of Volatiles in Amorphous Solid Water: Application to Interstellar Medium Environments, Astrophys. J., 2018, 863, 156.
13. C. J. Bennett and R. I. Kaiser, Laboratory Studies on the Formation of Ozone (O₃) on Icy Satellites and on Interstellar and Cometary Ices, Astrophys. J., 2005, 635, 1362.
14. B. Sivaraman, C. S. Jamieson, N. J. Mason and R. I. Kaiser, Temperature-dependent Formation of Ozone in Solid Oxygen by 5 keV Electron Irradiation and Implications for Solar System Ices, Astrophys. J., 2007, 669, 1414.
15. C. Ceccarelli, D. J. Hollenbach and A. G. G. M. Tielens, Far-Infrared Line Emission from Collapsing Protostellar Envelopes, Astrophys. J., 2009, 471, 400.
16. S. Viti, E. Roueff, T. W. Hartquist, G. Pineau des Forêts and D. A. Williams, Interstellar oxygen chemistry, Astron. Astrophys., 2001, 370, 557–569.
17. T. Hama, M. Yokoyama, A. Yabushita and M. Kawasaki, Role of OH radicals in the formation of oxygen molecules following vacuum ultraviolet photodissociation of amorphous solid water, J. Chem. Phys., 2010, 133, 104504.
18. R. Atkinson, D. L. Baulch, R. A. Cox, J. N. Crowley, R. F. Hampson, R. G. Hynes, M. E. Jenkin, M. J. Rossi and J. Troe, Evaluated kinetic and photochemical data for atmospheric chemistry: volume – gas phase reactions of O_x, HO_x, NO_x and SO_x species, Atmos. Chem. Phys., 2004, 4, 1461–1738.
19. D. Quan, E. Herbst, T. J. Millar, G. E. Hassel, S. Y. Lin, H. Guo, P. Honvault and D. Xie, New Theoretical Results Concerning the Interstellar Abundance of Molecular Oxygen, Astrophys. J., 2008, 681, 1318.
20. V. Wakelam, et al., A KInetic Database for Astrochemistry (KIDA), Astrophys. J., Suppl. Ser., 2012, 199, 21.
21. J.-I. Lo, S.-L. Chou, Y.-C. Peng, H.-C. Lu, J. F. Ogilvie and B.-M. Cheng, Thresholds of photolysis of O₂ and of formation of O₃ from O₂ dispersed in solid neon, Phys. Chem. Chem. Phys., 2018, 20, 13113–13117.
22. M. Minissale, E. Congiu, G. Manicò, V. Pirronello and F. Dulieu, CO₂ formation on interstellar dust grains: a detailed study of the barrier of the CO channel, Astron. Astrophys., 2013, 559, A49.
23. E. Hebrard, M. Dobrijevic, P. Pernot, N. Carrasco, A. Bergeat, K. M. Hickson, A. Canosa, S. D. Le Picard and I. R. Sims, How measurements of rate coefficients at low temperature increase the predictivity of photochemical models of Titan's atmosphere, J. Phys. Chem. A, 2009, 113, 11227–11237.
24. R. T. Garrod and T. Pauly, On the Formation of CO₂ and Other Interstellar Ices, Astrophys. J., 2011, 735, 15.
25. D. A. Williams and C. Cecchi-Pestellini, The Chemistry of Cosmic Dust, The Royal Society of Chemistry, 2016, pp. P001–304.
26. R. Johnson, P. Cooper, T. Quickenden, G. Grieves and T. Orlando, Production of oxygen by electronically induced dissociations in ice, J. Chem. Phys., 2005, 123, 184715.
27. B. Brooks, et al., CHARMM: the biomolecular simulation program, J. Comput. Chem., 2009, 30, 1545–1614.
28. T. Nagy, J. Y. Reyes and M. Meuwly, Multi-Surface Adiabatic Reactive Molecular Dynamics, J. Chem. Theory Comput., 2014, 10, 1366–1375.
29. J. Danielsson and M. Meuwly, Atomistic Simulation of Adiabatic Reactive Processes Based on Multi-State Potential Energy Surfaces, J. Chem. Theory Comput., 2008, 4, 1083.
30. D. R. Nutt and M. Meuwly, Studying Reactive Processes with Classical Dynamics: Rebinding Dynamics in MbNO, Biophys. J., 2006, 90, 1191–1201.
31. W. L. Jorgensen, J. Chandrasekhar, J. D. Madura, R. W. Impey and M. L. Klein, Comparison of simple potential functions for simulating liquid water, J. Chem. Phys., 1983, 79, 926–935.
32. N. Kumagai, K. Kawamura and T. Yokokawa, An Interatomic Potential Model for H₂O: Applications to Water and Ice Polymorphs, Mol. Simul., 1994, 12, 177–186.
33. N. Plattner and M. Meuwly, Atomistic Simulations of CO Vibrations in Ices Relevant to Astrochemistry, ChemPhysChem, 2008, 9, 1271–1277.
34. J. Huang and A. D. MacKerell, CHARMM36 all-atom additive protein force field: validation based on comparison to NMR data, J. Comput. Chem., 2013, 34, 2135–2145.
35. L. Bytautas, N. Matsunaga and K. Ruedenberg, Accurate ab initio potential energy curve of O₂. II. Core–valence correlations and relativistic contributions and vibration–rotation spectrum, J. Chem. Phys., 2010, 132, 1–15.
36. T. Hollebeek, T.-S. Ho and H. Rabitz, Constructing Multidimensional Molecular Potential Energy Surfaces from Ab Initio Data, Annu. Rev. Phys. Chem., 1999, 50, 537–570.
37. O. T. Unke and M. Meuwly, Toolkit for the Construction of Reproducing Kernel-Based Representations of Data: Application to Multidimensional Potential Energy Surfaces, J. Chem. Inf. Model., 2017, 57, 1923–1931.
38. T. Hama and N. Watanabe, Surface Processes on Interstellar Amorphous Solid Water: Adsorption, Diffusion, Tunneling Reactions, and Nuclear-Spin Conversion, Chem. Rev., 2013, 113, 8783–8839.
39. T. Kawabata, Detection of multiscale pockets on protein surfaces using mathematical morphology, Proteins, 2010, 78, 1195–1211.
40. J. Phillips, Microscopic aspects of Stretched Exponential Relaxation (SER) in homogeneous molecular and network glasses and polymers, J. Non-Cryst. Solids, 2011, 357, 3853–3865.
41. P. J. Steinbach, et al., Ligand binding to heme proteins: connection between dynamics and function, Biochemistry, 1991, 30, 3988–4001.
42. I. Langmuir, The mechanism of the catalytic action of platinum in the reactions 2CO + O₂ → 2CO₂ and 2H₂ + O₂ → 2H₂O, Trans. Faraday Soc., 1922, 17, 0621–0654.
43. C. N. Hinshelwood and N. V. Sidgwick, On the theory of unimolecular reactions, Proc. R. Soc. A, 1926, 113, 230–233.
44. D. D. Eley, On the solubility of gases. Part I. The inert gases in water, Trans. Faraday Soc., 1939, 35, 1281–1293.
45. R. J. L. Roy, LEVEL: a computer program for solving the radial Schrödinger equation for bound and quasibound levels, J. Quant. Spectrosc. Radiat. Transfer, 2017, 186, 167–178.
46. J. Phillips, Stretched exponential relaxation in molecular and electronic glasses, Rep. Prog. Phys., 1996, 59, 1133–1207.
47. J. A. Noble, E. Congiu, F. Dulieu and H. J. Fraser, Thermal desorption characteristics of CO, O₂ and CO₂ on non-porous water, crystalline water and silicate surfaces at submonolayer and multilayer coverages, Mon. Not. R. Astron. Soc., 2012, 421, 768–779.
48. J. Misewich, H. Zacharias and M. Loy, State-to-State Molecular-Beam Scattering of Vibrationally Excited NO from Cleaved LiF(100) Surfaces, Phys. Rev. Lett., 1985, 55, 1919–1922.
49. A. Wodtke, Y. Huang and D. Auerbach, Interaction of NO(v = 12) with LiF(001): evidence for anomalously large vibrational relaxation rates, J. Chem. Phys., 2003, 118, 8033–8041.
50. R. Jongma, G. Berden, T. Rasing, H. Zacharias and G. Meijer, Scattering of vibrationally and electronically excited CO molecules from a LiF(100) surface, Chem. Phys. Lett., 1997, 273, 147–152.
51. H.-C. Chang and G. E. Ewing, Infrared fluorescence from a monolayer of CO on NaCl(100), Phys. Rev. Lett., 1990, 65, 2125–2128.
52. D. J. Auerbach and A. M. Wodtke, in Dynamics of Gas-Surface Interactions: Atomic-level Understanding of Scattering Processes at Surfaces, ed. R. Díez Muiño and H. F. Busnengo, Springer Berlin Heidelberg, Berlin, Heidelberg, 2013, pp. 267–297.
53. S. T. Lunt, G. Marston and R. P. Wayne, Formation of O₂(a¹Δ) and vibrationally excited OH in the reaction between O atoms and HO_x species, J. Chem. Soc., Faraday Trans. 2, 1988, 84, 899–912.
54. J. R. Stine and J. T. Muckerman, J. Chem. Phys., 1976, 65, 3975.
55. L. D. Landau, Zur Theorie der Energiübertragung. II, Phys. Z., 1932, 2, 46.
56. C. Zener, Non-Adiabatic Crossing of Energy Levels, Proc. R. Soc. London, Ser. A, 1932, 137, 696.
57. D. Koner, R. J. Bemish and M. Meuwly, The C(³P) + NO(X²Π) → O(³P) + CN(X²Σ⁺), N(²D)/N(⁴S) + CO(X¹Σ⁺) reaction: Rates, branching ratios, and final states from 15 K to 20 000 K, J. Chem. Phys., 2018, 149, 094305.
58. B. Sivaraman, A. M. Mebel, N. J. Mason, D. Babikov and R. I. Kaiser, On the electron-induced isotope fractionation in low temperature ³²O₂/³⁶O₂ ices ozone as a case study, Phys. Chem. Chem. Phys., 2011, 13, 421–427.
59. D. Jing, J. He, J. R. Brucato, G. Vidali, L. Tozzetti and A. D. Sio, Formation of Molecular Oxygen and Ozone on Amorphous Silicates, Astrophys. J., 2012, 756, 98.
60. H. M. Cuppen, S. Ioppolo, C. Romanzin and H. Linnartz, Water formation at low temperatures by surface O₂ hydrogenation II: the reaction network, Phys. Chem. Chem. Phys., 2010, 12, 12077–12088.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: 10.1039/c8cp07474g

---