Diﬀusion of molecules in the bulk of a low density amorphous ice from molecular dynamics simulations †

The diﬀusion of molecules in interstellar ice is a fundamental phenomenon to take into account while studying the formation of complex molecules in this ice. This work presents a theoretical study on the diﬀusion of H 2 O, NH 3 , CO 2 , CO, and H 2 CO in the bulk of a low density amorphous (LDA) ice, while taking into account the physical conditions prevailing in space, i.e. temperatures below 150 K and extremely low pressure. This study was undertaken by means of molecular dynamics simulations. For CO 2 for which no experimental data were available we conducted our own experiments. From our calculations we show that, at low temperatures, the diﬀusion of molecules in the bulk of a LDA ice is driven by the self-diﬀusion of water molecules in the ice. With this study we demonstrate that molecular dynamics allows the calculation of diﬀusion coeﬃcients for small molecules in LDA ice that are convincingly comparable to experimentally measured diﬀusion coeﬃcients. We also provide diﬀusion coeﬃcients for a series of molecules of astrochemical interest.


I Introduction
Astrophysical observations of cold interstellar clouds, comets and icy planetary bodies have revealed the existence of ices mainly composed of water but also containing a few simple molecular species like CO 2 , CO, NH 3 , and CH 4 . 1,2 These ices form during the collapse of dense regions within interstellar clouds. 3,4 During this process, the density increases while the temperature decreases, allowing the condensation or formation of simple molecules on the dust grains present in these dense regions, also called star-forming regions. It is believed that later in the life of these regions of star formation, when the ices are warmed up until sublimation, these ices could catalyze the formation of complex organic molecules (COMs). 5,6 However, even at these temperatures, the scarcity of the molecules present in the ices and the resulting need for the reactants to diffuse and meet each other strongly limit ice chemistry. Whether or not long interstellar time scales can counterbalance the diffusion-limited solid-state reactivity in the diffusion-reaction equation derivation depends on both the reaction rate constants and the reactants diffusion coefficients. Diffusion of molecules in ices is therefore a key phenomenon to take into account when studying the formation of complex molecules in interstellar ices.
Diffusion-limited reactivity is also encountered in terrestrial ices, such as polar ice or atmospheric ice, 7,8 and drives the capacity of ice particles to modify atmospheric chemistry. An experimental study 9 devoted to the structure of interstellar ices showed that at very low temperature (T = 10 K) the condensation phase leads to the formation of a high density porous amorphous ice, which transforms into a low density amorphous (LDA) ice when the temperature increases (from 20 K to 120 K). 10 Molecular dynamics (MD) simulations have been applied to understand the evolution of a LDA ice to a high-density amorphous ice 11,12 while considering temperature and pressure effects or to study the high-density amorphous ice/crystalline ice I h compressiondecompression at a fixed low temperature. 13 Diffusion of small molecules, in crystalline ice, [14][15][16][17][18] clathrates 18,19 and hydrates, 20 has also been studied with molecular dynamics. The previously established interstitial mechanism for He atoms in ice I h 16 has been found inappropriate for the N 2 molecules. 17 From the comparison of diffusion coefficients of N 2 , O 2 , CO 2 , and CH 4 in hexagonal ices I h 14 the bond-breaking diffusion mechanism has been suggested with diffusion velocity being several orders of magnitude larger than that found in the case of an interstitial mechanism. In the interstitial mechanism He atoms have been found to migrate from a stable interstitial site to an adjacent site without a distortion of the lattice. In the bondbreaking mechanism, the hydrogen bonds in the lattice are broken and the molecule jumps between stable sites. The bondbreaking mechanism has also been proposed for the diffusion of formaldehyde in hexagonal I h ice. 15 The most recent work on CO 2 diffusion in clathrates has revealed that CO 2 diffusion is possible in these solids if water vacancy exists. 18 Diffusion has also been quantified in structure materials like H 2 , tetrahydrofuran clathrate hydrates 19 or for water in cyclodextrin hydrates. 20 For a more complete insight into molecular diffusion in clathrates we recommend the very recent review by English and MacElroy 21 on molecular simulation of clathrate hydrates. All the studies cited above have demonstrated the applicability of several different effective potentials such as TIP4P, 22 TIP4P-ice, 23 SPC, 24 MCY, 25 ST2, 26 SPC/E, 27 and Kawamura potential model to describe various types of ices and to study ice phase changes as well as the diffusivity of molecules in them. A fully atomistic, off-lattice kinetic Monte Carlo technique was applied to compute the diffusion to desorption barrier ratios of CO and CO 2 at crystalline ice surfaces at temperatures equal to 25 and 70 K. 28 The CO mobility, studied by the same technique, has been reported to be very sensitive to the amorphous ice surface morphology, with the CO molecules trapped by surface nanopores for low CO coverage. 29 To the best of our knowledge, diffusivity of molecules in the bulk of LDA ices has not been considered by means of molecular dynamics. The present work aims at bringing additional knowledge on diffusion by studying theoretically the bulk diffusion of H 2 O, NH 3 , CO 2 , CO, and H 2 CO in a LDA ice as a relevant model for an interstellar ice, at very low temperatures and no pressure.
Experimental studies of bulk diffusion in LDA ices are scarce too. 30 This paper is organized as follows. First the computational approach and protocols of the MD simulations are detailed (Section II). The experimental set-up for the CO 2 bulk diffusion measurements is presented (Section III). Then the computational results are discussed with respect to available experimental data from the literature and from our own experiments (Section IV).

II-1 The molecular dynamics protocol
In order to investigate the diffusion of small molecules in a low density amorphous ice mimicking an interstellar ice, classical molecular dynamics (MD) simulations using the OPLS-AA (optimized potentials for liquid simulations -all atom) [35][36][37] force field were performed. The s and e parameters of the Lennard-Jones (LJ) potential, implemented in the OPLS-AA force field, are summarized in Table 1. The geometrical parameters of the molecules are reported in Table 2. It is worth noting that for a correct description of the charges in the carbon monoxide, a dummy atom D was introduced between C and O, bearing a positive charge in order to account for the nucleophilic character of both C and O atoms. 38 The water molecules of the ice were described by the four point rigid TIP4P model 22,39 (transferable intermolecular potential four points) with the water molecules described by three sites corresponding to the atoms and an atom M located on the bisector of the HOH angle bearing the negative charge of the oxygen atom (Table 1). TIP4P was chosen after test calculations with the more recent TIP4PQ/ 2005 potential model developed by McBride et al. 40 and Noya et al. 41 that will be discussed later.
The diffusive molecules were described using for each atom the LJ parameters and a charge without any adjustment or extra site (except for CO). An equivalent potential model for a water molecule would be TIP3P 22 (where water is represented with three sites corresponding to the atoms with their respective charge). To validate our diffusion calculations for NH 3 , CO 2 , CO, and H 2 CO in the LDA ice we have also calculated the diffusion in the same LDA ice of a single TIP3P water molecule. As will be recalled later, such calculated diffusion will be designated in our paper as ''water diffusion'' to be distinguished from the ''water self-diffusion'' calculated for all the water molecules described by the TIP4P potential.
The cross-interaction between the solute molecule and the water molecule was defined by using the usual Lorentz-Berthelot combining rules (geometric mean of the C (6) and C (12) van der Waals coefficients) for each site-site interaction pair. Periodic boundary conditions were applied in all three directions. The cut-off for Lennard-Jones and real-space coulombic interactions was set to 0.9 nm. The particle mesh Ewald method was used to compute reciprocal electrostatic forces and energies. The equations of motion were integrated using the leap-frog algorithm with a time The initial arrangement of the water molecules is that of their liquid state. To achieve the correct LDA organization, an annealing procedure was carried out, following the methodology suggested by Martoňák et al. 11,12 First, the system was relaxed using a steepest descent algorithm. Second, it was placed in a Berendsen thermostat at 15 K with a coupling time constant of 0.1 ps, for a 100 ps simulation time. Third, the Parrinello-Rahman barostat at 0 bar was introduced with a coupling time constant of 0.5 ps and finally the equilibration of the system was continued for another 100 ps. This system was further equilibrated in the NPT ensemble, where the pressure and the temperature were fixed by coupling the system to the barostat and the thermostat. The temperature was raised to 170 K by steps of 15-20 K with an equilibration of the system for 10 ns, at each step. At 170 K, the system was equilibrated for another 1 ms simulation in order to reach the LDA organization. The cooling was done by decreasing the temperature from 170 K to 60 K by steps of 15-20 K and equilibrating the system for 100 ns at each step. With this procedure, thermalized ice structures were obtained for each temperature. Finally, performing an extra 500 ns simulation equilibrated each one of these structures. To study the diffusion of the solute molecules in the LDA ice, 5 ms MD trajectories were calculated starting from the equilibrated systems at 170 K, 150 K, 135 K, 120 K, 105 K, 90 K, 75 K, and 60 K.
To obtain equilibrated LDA structures at T = 200, 225, 250, and 275 K a similar procedure was performed by starting from the ice equilibrated at T = 170 K. However, as it will be discussed later in Section IV-1, for these high temperatures the ice was losing its LDA organization.
The ice densities were calculated all along the simulated annealing procedure. The initial density was that of the liquid state created by the initial GROMACS arrangement. It decreased during the heating phase and stabilized during the long equilibration at 170 K to a value of 0.958. At each temperature along the ice thermalization, the density increased again as the ice was cooled, but remained below 1 as required for a LDA ice. 11,12,[43][44][45] The evolution during our simulations of the densities can be viewed in Fig. S1 of the ESI. † As the size of the simulation box was previously reported to influence the diffusion coefficient calculations, 46 we have tested a larger box of 2.57 nm containing 540 water molecules. The procedure described above was followed. Densities calculated with the larger simulation box of 540 water molecules, in the range of 0.972 at 60 K to 0.949 at 170 K, showed no difference with respect to densities calculated for the 350-380 water cell.
Densities for the equilibrated structures were also calculated using the TIP4PQ/2005 potential. As can be viewed in Fig. S4 (ESI †), they are higher than the densities calculated with the TIP4P model with values above 1. These high densities do not characterize LDA ices but rather crystalline or high-density amorphous ices. Indeed TIP4PQ/2005 was improved by increasing the charges on the hydrogen and the M atoms in order to reproduce the densities and the structures of the crystalline ices.

II-3 The diffusion coefficient calculations
The diffusion coefficient, D T , of a target molecule T can be derived from the mean square displacement (MSD) as a function of time t by fitting r i (t), the position of a particle i (from the ensemble of T diffusing molecules) as a function of time obtained against the Einstein relation: 47 The temperature dependence of the diffusion coefficients is obtained by fitting the diffusion coefficient curve against an Arrhenius law:

RT
where E diff is the activation energy of the diffusion process and D 0 is the pre-exponential factor.
Diffusion coefficients were calculated from the MSD following the methodology described above and implemented in GROMACS. 42

III-1 The ice film preparation
The experiments were performed using the RING experimental set-up as described elsewhere. 48 A gold plated copper surface is maintained at low temperature using a closed-cycle helium cryostat (ARS Cryo, model DE-204 SB, 4 K cryogenerator) within Table 2 Geometry parameters for all molecules considered in the MD simulations a high-vacuum chamber at a few 10 À9 hPa. The sample temperature is measured with a DTGS 670 silicon diode with a 0.3 K uncertainty. The temperature is controlled using a Lakeshore Model 336 temperature controller and a heating resistance. Infrared absorption spectra are recorded in the reflection absorption mode by means of Fourier transform reflection absorption infrared (FTIR-RAIRS) spectroscopy using a Vertex 70 spectrometer with either a DTGS detector or a liquid N 2 cooled MCT detector. A typical spectrum has a 1 cm À1 resolution and is averaged over a few tens of interferograms. Water vapor is obtained from deionized water, which was purified by several freeze-pump-thaw cycles, carried out under primary vacuum. Gas-phase CO 2 is commercially available in the form of 99.9995% pure gas from Air Liquide. Gas-phase H 2 O and CO 2 are mixed together at room temperature in a primary pumped vacuum line using standard manometric techniques, with CO 2 mixing ratios of a few percent. Then the homogeneously mixed gas-phase mixture is sprayed onto the cold gold plated copper surface at a normal incidence at 80 K to give a homogeneously mixed CO 2 : H 2 O ice mixture, with a large excess of amorphous solid water (ASW), as displayed in Fig. 1. The 80 K deposition temperature is chosen to obtain a compact morphology of the amorphous ice (c-ASW) and prevent the diffusion on the open pores surface. At 80 K the pore collapse was measured at 1.7 Â 10 À4 s À1 , 33 i.e. a characteristic time of approximately an hour and a half. Fig. 2 shows the characteristic spectra we get for CO 2 : H 2 O mixtures.
The column density N (molecules cm À2 ) of each molecular species is derived from the IR absorption spectrum right after deposition from the film IR spectra, as seen in Fig. 2, using the expression where the optical depth t u is equal to ln(10) times the integrated absorbance, and A is the band strength in cm molecule À1 . Carbon dioxide is identified by its asymmetric stretching mode at 2342 cm À1 and its bending mode band at 667 cm À1 . The band strength for the CO 2 asymmetric stretching band was measured to be 7.6 Â 10 À17 cm molecule À1 for the pure solid, 49,50 while in water ice a value of 1.4 Â 10 À17 cm molecule À1 was found. 51 Water ice has three characteristic bands at 3280, 1660 and 760 cm À1 corresponding to the OH stretching, HOH bending and libration modes respectively. The corresponding band strengths are 2.1 Â 10 À16 , 3.1 Â 10 À17 and 3.1 Â 10 À17 cm molecule À1 , respectively. 50 There is an approximate 30% uncertainty on the band strengths and so on the calculated column densities.
The ASW ice film thickness l the CO 2 molecules need to diffuse through is an important parameter in our experiment. The method we used to measure it is based on the quantity of matter as determined from the IR absorption bands, using the H 2 O OH stretching band. The ASW thickness l is derived from the measured column density N (molecule cm À2 ) using r = 0.94 g cm À3 as the amorphous ice density, and using with N A being the Avogadro number 18 g mol À1 is the molar mass for H 2 O. The cos(181) comes from the 181 incidence angle between the FTIR beam and the ice film normal angle. The factor a half comes from the reflection mode, which probes double the column density. The uncertainty on the ASW thickness l is therefore mainly given by the uncertainty on the band strengths, around 30%. For thin ice films (few hundreds of nm thick) it is more precise to determine the thickness from IR spectra than from He-Ne laser interference.
Since the ice sample can be modeled as a cylinder of a few centimeters diameter and a few hundreds of nanometers thick, and since the IR beam has a smaller diameter, we can reasonably assume that the CO 2 molecules are mainly diffusing along the x direction within the ASW layer and that a negligible amount of them escape from the cylinder side. The diffusion Fig. 1 Scheme of the ASW film of thickness d, where the CO 2 is homogeneously mixed in the dominant ASW at the initial time. The onedimension diffusion of the CO 2 molecules along the x direction is monitored at a fixed temperature by recording the evolution of its characteristic IR absorption bands as a function of time using FTIR spectroscopy. The diffusion boundary conditions at the surface are a null flow at x = 0 and a flow set by the desorption rate at x = d.  Table 3.

III-2 The kinetics of the ASW film pore collapse
The morphology of the ice depends on the temperature of the gold surface the water vapor is deposited on. Once deposited the ASW (porous p-ASW or compact c-ASW for low-and highdeposition temperatures respectively) is metastable and will tend to reach a thermodynamically equilibrium state. The pore structure will collapse and the amorphous structure will reorganize into an organized crystalline network. This way the equilibrium can be decomposed by phase transitions, from an initial high-density amorphous ice Ia h (temperature deposition at 15 K) to a lowdensity amorphous ice Ia l between 38 K and 68 K, and possibly a third amorphous form Ia r preceding the crystallization into cubic ice Ic 1. Because of the deposition temperature we choose, our ice does not undergo the high density to low density transition. Ice porosity is visible from the small OH dangling bands at 3720 and 3696 cm À1 , indicating a large effective surface. 52,53 Once slightly warmed up these dangling bands quickly disappear, indicating the start of the pore collapse and the decrease of the effective surface. The ice samples we prepared have a small surface to volume ratio as indicated by the absence of OH dangling bands in our spectra as seen in Fig. 2. The pore collapse also induces modifications of the OH stretching band. The pore collapse kinetics can be estimated from the changes in the OH stretching band. 33,54 Fitting the experimental kinetic rates measured in ref. 33 between 40 K and 120 K with an Arrhenius law gives an approximate 0.8 AE 0.1 kJ mol À1 reorganization energy, with a 6.5 AE 1.5 Â 10 À4 s À1 pre-exponential factor. This reorganization is fast. For example at the onset of crystallization, at 150 K, it is around 18 minutes. 33 In the interstellar medium, at 10 K, it is approximately 9 months, which is nothing on an interstellar timescale. Above 120-140 K the crystallization kinetics takes over the reorganization process which precedes it. Both the deposition temperature we have chosen, 80 K, and the temperature interval we are working in ensure that the ASW ice is compact, with a small surface to volume ratio, and that we are measuring volume diffusion mainly.

III-3 The isothermal kinetic experiments on CO 2
In isothermal kinetic (IK) experiments, right after deposition, the ice film is heated as fast as possible to a fixed target temperature T, in typically a few tens of seconds. Once the target temperature T is reached, we set the initial time t = 0 s of our isothermal kinetics. Assigning a ''time zero'' can be difficult because of the time it takes to reach the isothermal temperature.
If the deposition temperature is close to the target temperature the uncertainties are small, ca. a few seconds, but the error can be significant for low temperature deposition, ca. a few tens of seconds, which is nevertheless small compared to hours long kinetics. In an IK experiment, the CO 2 molecules diffuse within the water ice film, up to the top surface of the ice, and then desorb. The diffusion-desorption of CO 2 along the x direction at the fixed temperature T is monitored by recording its abundance decay from its characteristic absorption band at 2342 cm À1 as a function of time, until it reaches a plateau, as shown in Fig. 3 for several temperatures. The IR decay curve is directly related to the CO 2 molecules diffusion in the compact ASW ice. The decay curves are normalized to the initial CO 2 abundance. The temperature range that can be accessed to is limited. The lowest temperature is set by the time duration and the long-term stability of the experiment, i.e. by the deposition of residual water on top of the ASW film and the consequent increase in the ASW ice thickness. Within our vacuum conditions, ca. at a few 10 À9 hPa, the time to grow ca. 10% of the initial thickness corresponds to approximately three days. The highest accessible temperature is set by the ASW substrate desorption compared to the time needed to acquire one spectrum. The H 2 O desorption rate constant is k(H 2 O) = 10 15 s À1 Â exp(46.6 kJ mol À1 /RT), 55 which corresponds to a residence time of 240 s at 140 K and 17 s at 150 K. At 150 K, the ASW desorption rate is too high and prohibits IK studies. Fig. 4 shows that at T = 150 K the CO 2 decay curve is not due to diffusion but due to the water ice substrate desorption. Therefore we will consider IK experiments between 95 K and 140 K only.  Fig. 3 Decay curve of CO 2 at temperatures ranging from 100 to 140 K using two different deposition temperatures (80 K and 100 K) as presented in the legend. An example of the fit, using the formula introduced in Section III-3, of the CO 2 decay curve at 135 K is plotted with a dashed line.
The reproducibility and the dispersion of the experiments are estimated by measuring the decay of CO 2 at 140 K (experiments 8, 9, 10 and 11 in Table 3). The four experiments show decay rates of the same order of magnitude (8.4 Â 10 À4 s À1 , 4.2 Â 10 À4 s À1 , 3.6 Â10 À4 s À1 , 1.6 Â 10 À4 s À1 , respectively).
The CO 2 abundance decay curves are fitted using a onedimension diffusion equation in a plane sheet, where the boundary condition at the x = 0 surface is given by a null flow and where the boundary condition at the free surface is set by the CO 2 desorption rate. The initial concentration is C 0 at t = 0. We chose temperatures for which the CO 2 residence time on top of the ice surface (x = 1) is much smaller than the diffusion time, so that we have an infinite desorption rate and therefore a null concentration at x = 1. Since we never totally deplete the initial CO 2 reservoir, we can make the assumption that the concentration at x = 0 is kept constant at C 0 . In that case the total amount of CO 2 present at time t in the plane ice sheet, normalized to the initial quantity, is 56 The quantity M(t)/M(t = 0) is the experimentally measured quantity as shown in Fig. 3.
For some of the temperatures, it was not possible to derive a diffusion coefficient, but rather an upper limit for the diffusion coefficient for a given thickness and temperature. The experimental results are presented in Table 4 and can be viewed in Fig. 9b. They are discussed with respect to the calculated values in Section IV-4. One must emphasize that the decay curve fitting depends on the estimated final value of the curves, which is taken arbitrarily, and not set to zero as we might theoretically expect. This remaining CO 2 problem has been previously observed in Mispelaer et al. 33 and will be discussed in Section IV-4.

IV-1 The LDA ice structure
The LDA ice structure equilibrated at various temperatures and obtained by employing the TIP4P effective potential, including for the H 2 O solute molecule, is defined by us as a ''reference system'', and this notation will be used later in the text.
In order to ensure that the ice structures, resulting from the methodology described in Section II-2, have the required local structure of LDA ices, we computed for the ''reference system'' the O-O radial distribution function (RDF), g OO (r), at 170 K, 105 K and 60 K. Fig. 5 shows that our calculated g OO (r) behaves similarly for the three temperatures. The curves exhibit two peaks corresponding to the two first solvation shells at 2.75 Å and 4.5 Å respectively, and a deep minimum at 3.1 Å. These peaks (positions and relative intensities) are very similar to those found experimentally by Finney et al. 43 and Bowron et al. 45 using neutron diffraction, and from molecular dynamics simulations by Martoňák et al. 11,12 The existence of these peaks is an indication of a local structural organization in the LDA ice, whereas the amorphous character of the ice model can be concluded from the g OO (r) at r 4 5 Å, where the fast averaging of the oxygen density evidences the non-existence of crystallinity as expected for a LDA ice. The peak intensities vary with the temperature. This is more pronounced for the first peak at 2.75 Å: a lower intensity being associated with a decrease of the local structural organization. Our calculated g OO (r) indicates, as expected, a higher structural organization of the ice for lower temperature. The positions of the peaks are the same for the three temperatures demonstrating that the LDA ice structure is preserved in the interval of T = 60 to 170 K.
Our calculated RDF, at 170 K, using the TIP4PQ/2005 and the TIP4P potential models are compared in Fig. S2 (ESI †).

View Article Online
One can see from this comparison that the two peaks corresponding to the first solvation shells are shifted toward smaller distances (2.65 Å and 4.35 Å) in the case of TIP4PQ/2005, confirming the nonadequacy of this potential to describe LDA ices. We note here that the radial distribution functions for the larger cell of 540 water molecules were also computed. There is no difference with the reference system as follows from the results reported in Fig. S3 (ESI †).
In order to evaluate the porosity of the modeled ice, a cavity distribution calculation was carried out using the SURFNET program. 57 This program allows for the calculation of the position and volume of the cavities present in the ice. Fig. 6 is a snapshot of our LDA ice where water and cavities are shown.
The cavity distribution was then built by calculating the size and the volume of the cavities for each step and averaging for 1000 steps of the molecular dynamics. The distributions of the cavity radii and volumes are shown in Fig. 7 at 170 K and 105 K. Fig. 7 shows that the majority of the cavities of our ice model are very small with radii below 2 Å. They are much smaller than a ''vacancy'' representing the volume leftover when a water molecule is absent and for which, from our calculated distance of the first solvation layer (Fig. 5), we expected a radius of 2.8 Å minimum. The cavities we calculate are likely due to the statistical motion of the water molecules of our ice water box creating small interstices. Our LDA ice is evidently not porous since pores would appear with at least the radius of a vacancy. We have also calculated the densities of the ice at various temperatures. They range between 0.98 g cm À3 (T = 60 K) and 0.96 g cm À3 (T = 170 K). They are very close to the density of 0.97 g cm À3 determined by Martoňák et al. 11,12 at 80 K and zero pressure for the same type of ice, and to the experimental density of 0.94 g cm À3 measured at 117 K. 44 The size of the cavities of our LDA ice model and its density confirm its compactness and therefore its relevance for modeling astrophysical ices for T = 60-170 K.

IV-2 Water self-diffusion and diffusion of a water molecule in a LDA ice
In this section, we are comparing the calculated ''water selfdiffusion'' in a LDA ice with ''water diffusion'' in the same ice. We recall that in the former case, all molecules of the simulation box, all described with the TIP4P effective potential (i.e. the ''reference system''), are considered as diffusive molecules. In the later case, only the motion of one solute H 2 O molecule in the LDA ice is considered. In this last case, the diffusing water molecule is described using the TIP3P effective potential while the remaining water molecules composing the ice are described using the TIP4P effective potential. Such a comparison is undertaken to validate our approach for the diffusion coefficient calculations of a solute molecule in a LDA ice (Section IV-3). Molecular diffusion trajectories of 1 ms were considered to calculate the water self-diffusion coefficients while these trajectories were extended to 5 ms for the ''water diffusion'' coefficients calculations.
Our calculated O-O radial distribution functions for the ''reference system'' and the LDA ice containing one solute  molecule (Fig. S6 of ESI †) are identical, both in very good agreement with available data from the literature. 11,12,[43][44][45] Fig. 8 reveals that the calculated diffusion coefficients follow different regimes: a ''high temperature'' regime above 170 K, a ''medium temperature'' regime between 170 K and 90 K, and a ''low temperature'' regime below 90 K. Below 90 K, an asymptotic behavior is observed. These plateaus are an indication that at low temperatures, the mean square displacements from the MD simulations are too small to determine the diffusion coefficients with a good precision. Indeed the MSD derived from the trajectories at 170 K and 60 K and reported in Fig. S7 and S8 (ESI †), respectively, show at 170 K a linear MSD behavior with time, while at 60 K the MSD curve is flat and difficult to extract from the noise. We conclude that the present MD calculations are not able to give diffusion coefficients for temperature below 90 K, as are the IK experiments.
Our simulations between 170 and 90 K are of interest for the present study. The ice model is of LDA type as confirmed from the calculated g OO (r) and r LDA at these temperatures (as well as calculated cavity distributions at 170 K and 150 K). These temperatures correspond to those where diffusion and thermal reactivity could occur at reasonable timescales in the interstellar medium. 6 Indeed, it has been shown that in the laboratory for conditions mimicking those in the interstellar medium, NH 3 and CO 2 react to form ammonium carbamate in 3 hours at 90 K. 58 In our limited temperature interval, the calculated diffusion coefficients exhibit an Arrhenius behavior (linear from Fig. 8), decreasing with temperature as expected, very similarly for both the self-diffusion and the water molecule diffusion. The calculated water self-diffusion coefficients and the calculated diffusion coefficients of a water molecule in the same LDA ice are reported in Table 5. They compare very satisfyingly regarding the uncertainty of our calculations (of an order of magnitude at most from our calculations).
The self-diffusion coefficients derived from molecular dynamics trajectories computed for the larger 2.57 nm water cell are compared to the self-diffusion coefficients calculated for the 2.25 nm reference system (Fig. S9 of ESI †). Their similarity as well as the very good agreement of radial distribution functions and densities, computed for both water box sizes, justified our choice to use the smaller box with 350-380 water molecules to study the diffusion of molecules in a LDA ice.
Using the Arrhenius formula defining the temperature dependence of the diffusion coefficients, the activation energies were estimated from the self-diffusion and water molecule diffusion calculations for T = [170-90 K]. They are 15 AE 5 kJ mol À1 and 13 AE 5 kJ mol À1 with pre-exponential factors of (7 AE 1) 10 À6 cm 2 s À1 and (9 AE 1) 10 À7 cm 2 s À1 , respectively, and are reported in Table 7.
Following the Vogel-Fulcher-Tammann (VFT) law for glass, as given below, Smith et al. 32 obtained an activation energy, a pre-exponential factor and a reference temperature of 7.2 AE 0.8 kJ mol À1 , 2.8 AE 2 Â 10 À3 cm 2 s À1 , and 119 AE 3 K, respectively.
The experimental self-diffusion coefficients determined by Smith et al. 32 in an amorphous solid water ice are in orange in Fig. 8. They range from 10 À15 to 10 À12 cm 2 s À1 for temperature between 150 K and 160 K. The variation with temperature of these self-diffusion coefficients is attributed to an amorphous solid ice glass transition into a quasi-liquid phase prior to crystallization. 59,60 However at 150 K, crystallization kinetics have been proved to be a minute-scale phase transition, 61 so that the diffusion coefficient of 10 À15 cm 2 s À1 determined experimentally at 150 K might account for an ice in a crystalline phase explaining its discrepancy with our calculated one at the  same temperature in the LDA ice but its agreement with our calculated one in a crystalline ice (10 À14 cm 2 s À1 at 150 K These values are much lower than our calculated amorphous water ice diffusion coefficient of 10 À12 cm 2 s À1 at 150 K while closer to our calculated one in a crystalline I h ice of 10 À14 cm 2 s À1 at 150 K (see also Fig. 8). As expected the self-diffusion of water in an amorphous ice is a couple of orders of magnitude faster than in a crystalline ice.
To allow for further comparison of our theoretical calculations with data from the literature 16,63 we have calculated the water ice ''self-diffusion'' coefficients for higher temperatures. Between 200 and 270 K they range from 10 À7 to 10 À5 cm 2 s À1 comparing satisfyingly with the experimental data of Goto et al. in crystalline ice 63 of 10 À6 -10 À5 cm 2 s À1 measured between 230 and 260 K (see Fig. 8). Ikeda et al. 16 calculated water diffusion coefficients in ice I h to be 10 À6 -10 À5 cm 2 s À1 at T = 200-270 K. The radial distribution functions for these high-temperature simulations are shown in Fig. S5 (ESI †). One can see that the RDFs lose their structures, as temperature increases. The densities also increase as shown in Fig. S4 (ESI †): a tendency expected for an ice to liquid water transition.
We have carried out a supplementary fit using our water selfdiffusion coefficients including the supercooled liquid points and considering temperatures from 90 K to 275 K. Comparing the resulting fit with the fitted values of Smith et al. 32 following the VFT law (in orange in Fig. 8), we can state that as expected diffusion in amorphous ice, which is generally understood as diffusion in liquid water, is better described by an Arrhenius law. 64 The successful inclusion of the supercooled liquid zone in the fit argues for the same diffusion mechanism for both the LDA ice and the supercooled water; this certainly characterizes a viscous-like diffusion.
The similarity of our results (considering our one order of magnitude uncertainty) for both the self-diffusion and water diffusion approaches and the agreement of our results with measured values argue for the validity of our protocol to calculate diffusion coefficients of molecules by considering a molecule diffusing in a LDA ice of 350 to 380 water molecules for MD trajectory of 5 ms. This protocol was used to calculate the diffusion coefficients of CO, CO 2 , NH 3 , and H 2 CO in the same LDA ice.

IV-3 Diffusion of CO, CO 2 , NH 3 , and H 2 CO in a LDA ice
We have calculated diffusion coefficients for CO, CO 2 , NH 3 and H 2 CO, in the LDA ice model. Their diffusion coefficients are reported in Fig. 9a-d, along with water self-diffusion coefficients, for temperature ranging between 170 K and 90 K. We did not report our calculated diffusion coefficients below 90 K because they are not reliable, as discussed above. Fig. 9a-d show that in the 90-170 K temperature range, the diffusion coefficients of CO, CO 2 , NH 3 and H 2 CO present the same dependence on temperature as the water self-diffusion despite their difference in mass, geometry, polarizability, dipole moment and ability to make hydrogen bonding. A deeper analysis of these diffusion coefficients is possible from Table 6 where they are reported together with the available experimental diffusion coefficients of NH 3 , CO CO 2 and H 2 CO measured in LDA ices.
Using the above Arrhenius formula, activation energies were estimated from the diffusion calculations for T = [170-90 K]. They are reported in Table 7.
From Table 6a-d one can see that the calculated diffusion coefficients are rather similar while corresponding to molecules with different physical properties (mass, dipole moment, polarizability, hydrogen bonding). In between 120 K and 135 K they are of the order of 10 À12 cm 2 s À1 for H 2 O self-diffusion, 10 À13 -10 À12 cm 2 s À1 for the diffusion of CO 2 , 10 À12 cm 2 s À1 for NH 3 , 10 À11 -10 À12 cm 2 s À1 for H 2 CO and 10 À13 cm 2 s À1 for CO. At 150 K they are of the order of 10 À11 cm 2 s À1 for H 2 O selfdiffusion, 10 À12 cm 2 s À1 for the diffusion of CO 2 , 10 À11 cm 2 s À1 for NH 3 and H 2 CO and 10 À13 cm 2 s À1 for CO. Taking into account the one order of magnitude uncertainty of our calculations, we confirm the trend already observed from Fig. 9a to d that the diffusion in these LDA ices is driven by the diffusion of the water molecules of the ice, regardless of the solute molecule diffusing in. The classification suggested by Collings et al. 34 for desorption can obviously not be transposed to bulk diffusion in amorphous water ice.

IV-4 Discussion
Let us now compare our calculated bulk diffusion coefficients of CO, CO 2 , NH 3 and H 2 CO with the experimentally measured ones (or tentatively measured in the case of CO) in the bulk of LDA ices. Table 6 shows that within the experimental and theoretical uncertainties, i.e. one order of magnitude for both, we can argue for a rather satisfying agreement between the calculated diffusion coefficients and the experimental one for CO 2 , NH 3 and H 2 CO, validating our molecular dynamics simulations and therefore our suggestion above of a diffusion mechanism driven by the diffusion of the water ice molecules. The case of CO 2 is important since its diffusion does not involve hydrogen bond breaking. We do not observe the same agreement for CO for which the experimentally measured diffusion is expected to be dominated by surface diffusion even at low temperature because of its low desorption energy (9.8 kJ mol À1 ) and the experimental difficulties to get compact ices at such low temperatures. Indeed, at these low temperatures, the ices are very porous, allowing for the diffusion of CO through the pores towards the surface. The diffusion coefficients measured for CO, much higher than what would be expected for bulk diffusion, are consistent with the surface diffusion coefficients calculated by Karssemeijer for CO and CO 2 28,29 (see Fig. 9a and c).
Our calculated diffusion coefficients for NH 3 in a crystalline ice (Table 6) are of the order of 10 À14 cm 2 s À1 between 90 K and 150 K. At 140 K Livingston 30 measured in a crystalline ice a bulk diffusion coefficient for NH 3 of 4.0 Â 10 À10 cm 2 s À1 . This value is higher than our diffusion coefficient calculated in the crystalline ice (but closer to the one in the LDA ice). Two explanations are possible. (i) As discussed in Livingston et al. 30 the high diffusion coefficient measured might be related to a crystalline ice lattice disruption. Breaking the lattice would generate vacancies causing the ammonia to move via a faster surface diffusion. (ii) Taking into account that during our molecular dynamics no destruction of the crystalline lattice is observed, it is clear that our calculations cannot account for a bulk diffusion mechanism mediated by vacancy formation.
We have derived activation energies for the diffusion of CO, CO 2 , NH 3 and H 2 CO from our Arrhenius fits. They are given in Table 7. One can notice that the activation energies might be related to the polarity of the molecule with CO 2 presenting the lowest activation energy and H 2 O and NH 3 the highest one. The height of the activation energies might also be related to the ability of the molecules to build hydrogen bonds. Our activation energies compare rather well with the experimental values of Mispelaer et al.: 33 for instance, for CO and H 2 CO the measured activation energies are respectively 1 and 12 kJ mol À1 which are in the same range as ours (8 and 9 kJ mol À1 , respectively). A larger discrepancy is observed for ammonia with an experimental value of 71 kJ mol À1 and a theoretical one of 17 kJ mol À1 . However the agreement between corresponding measured and calculated diffusion coefficients suggests that this discrepancy for the activation energies might be due to fitting artifact and due to the limited number of experimental values used for the fitting.
Our derived activation energies are totally correlated to the corresponding derived pre-exponential factors. To get around this bias we have derived another set of activation energies for CO, CO 2 , NH 3 and H 2 CO using a unique value for the preexponential factor, the value of 0.22 cm 2 s À1 derived from fitting the water self-diffusion coefficients including the selfdiffusion of the supercooled water. The new set of activation energies obtained for CO, CO 2 , NH 3 , and H 2 CO as well as for H 2 O, for T = 90-170 K, are given in Table 7. These activation energies are all of the order of 25 kJ mol À1 , arguing for a diffusion mechanism equivalent for all species and therefore for a solvent-driven mechanism. is similar to what was found for other molecules, strengthening the argument of a water solvent-driven diffusion. The non-total depletion of desorbing molecules has been previously ascribed to a trapping phenomenon. This trapping can be viewed as the creation of a local structure around the diffusive molecule similarly in principle to the clathrate example. Thus, this over-structuration drastically slows down the diffusion, ''entrapping'' a fraction of the molecule. It can also be viewed as a change of the diffusion coefficient (couple of order of magnitude lower) during the phase change. 32 The astrophysical implication of such trapping raises a real hot question, because it could strongly influence the ice molecular compositions: volatiles can be trapped in the ice above their desorption temperatures, and they can possibly be involved in high barrier reactions. It is therefore important to understand how far this ''trapping'' slows down the diffusion process. Indeed, when a week-long experiment cannot see the full depletion, can we expect a similar ''trapping'' in molecular cloud ices warmed up slowly over thousands of years? A theoretical molecular investigation might bring more answers and enlighten this challenging issue.
We can relate the experimental diffusion coefficient to a characteristic time using the Einstein-Smoluchowski formulation of the one-dimension diffusion equation for Brownian motion given below: If we relate the mean squared displacement to the average spacing between two binding sites, ca. 2.8 Å, we can estimate a characteristic thermal hopping time for a given temperature. At 90 K, our calculated diffusion coefficients for almost all molecules including H 2 O are about 10 À14 cm 2 s À1 . This gives a 100 millisecond characteristic thermal hopping time at 90 K. A 100 nm diameter grain, with a typical site density of 10 15 sites cm À2 , has approximately 10 6 sites per monolayer. Thus, scanning azimuthally a whole grain monolayer takes around 10 5 seconds or around 30 hours. Crossing radially 100 ML (ca. 30 nm) thin ice to reach the surface takes around 17 minutes. These rough estimations show that above 90 K, the diffusion of neutral molecules is significant at long timescales, especially for the star formation timescales (10 5 to 10 6 years depending on the mass of the star). For the sake of comparison, at 10 K, a hydrogen atom would scan an ASW surface within a few days. 65

V Conclusion
We have demonstrated that molecular dynamics enables one to calculate diffusion coefficients of small molecules such as CO, CO 2 , NH 3 and H 2 CO in a LDA ice and that the calculated coefficients compare satisfyingly with experimentally measured coefficients, within the experimental and theoretical uncertainties. From these calculations we are able to suggest a bulk diffusion mechanism at low temperature driven by the diffusion of the water molecules in the ice. The validation of the molecular dynamics approach from experimental measurements is of prime importance if we want to extend calculations to other molecules and to lower temperatures. Theory and experiments are complementary as experiments measure macroscopic diffusion while molecular dynamics calculations give a microscopic insight into the diffusion. However, they both suffer from the same limitations at very low temperature. Experiments at low temperature are limited by the IK experiment maximum duration time, as discussed in the experimental section, while theory is limited by computational time; yet, calculations enable us to investigate slightly lower temperatures. Moreover, it has been possible to calculate CO bulk diffusion while it was not possible to measure it experimentally.
The knowledge of bulk diffusion coefficients, which we have demonstrated to be obtainable from calculations, should now allow astrochemical models to introduce diffusion kinetic limitation for reactivity and desorption in multilayered ice. The volume of the ice has long been considered chemically inert and considered independently. Comparing our CO and CO 2 values with those measured or calculated 29,33 at lower temperature, we have outlined the large difference between bulk and surface diffusion. Our work provides the tools to refine this statement as a function of temperature and scarcity of the reactants.