Kinetic Monte Carlo simulations of water ice porosity: extrapolations of deposition parameters from the laboratory to interstellar space

Aspen R. Clements *a, Brandon Berk a, Ilsa R. Cooke a and Robin T. Garrod ab
aDepartment of Chemistry, University of Virginia, Charlottesville, VA 22904-4319, USA. E-mail:
bDepartment of Astronomy, University of Virginia, Charlottesville, VA 22904-4319, USA

Received 31st August 2017 , Accepted 23rd January 2018

First published on 1st February 2018

Dust grains in cold, dense interstellar clouds build up appreciable ice mantles through the accretion and subsequent surface chemistry of atoms and molecules from the gas. These mantles, of thicknesses on the order of 100 monolayers, are primarily composed of H2O, CO, and CO2. Laboratory experiments using interstellar ice analogues have shown that porosity could be present and can facilitate diffusion of molecules along the inner pore surfaces. However, the movement of molecules within and upon the ice is poorly described by current chemical kinetics models, making it difficult either to reproduce the formation of experimental porous ice structures or to extrapolate generalized laboratory results to interstellar conditions. Here we use the off-lattice Monte Carlo kinetics model MIMICK to investigate the effects that various deposition parameters have on laboratory ice structures. The model treats molecules as isotropic spheres of a uniform size, using a Lennard-Jones potential. We reproduce experimental trends in the density of amorphous solid water (ASW) for varied deposition angle, rate and surface temperature; ice density decreases when the incident angle or deposition rate is increased, while increasing temperature results in a more-compact water ice. The models indicate that the density behaviour at higher temperatures (≥80 K) is dependent on molecular rearrangement resulting from thermal diffusion. To reproduce trends at lower temperatures, it is necessary to take account of non-thermal diffusion by newly-adsorbed molecules, which bring kinetic energy both from the gas phase and from their acceleration into a surface binding site. Extrapolation of the model to conditions appropriate to protoplanetary disks, in which direct accretion of water from the gas-phase may be the dominant ice formation mechanism, indicate that these ices may be less porous than laboratory ices.

1 Introduction

Amorphous solid water (ASW) is a major component of astronomically important ice structures – including the ices found on the surfaces of interstellar dust-grains,1–5 as well as ices in comets,6,7 planetary rings, and satellites.8 Infrared observations of protostellar cores and quiescent (i.e. non-star-forming) clouds indicate that water is the most abundant component of solid ice in the interstellar medium;2,3,9 its structure will therefore play a significant role in regulating the surface chemistry in interstellar regions.

The dust grains (of canonical radius ∼0.1 μm) upon which interstellar ices reside comprise around 1% of the mass of interstellar clouds in our galaxy, and are composed primarily of carbonaceous and silicaceous compounds.2,10–12 At the low temperatures in dense interstellar clouds (on the order of 10 K), atoms and molecules from the surrounding gas can stick to the grain surfaces. The resulting physisorbed species may further desorb, diffuse, and/or react with other species. In this way, the accretion of gas-phase material onto the grains leads to the build up of ice mantles of thicknesses up to around 100 monolayers (ML).13 Besides water, much of the rest of the ice is composed of CO and CO2, as well less abundant species such as CH3OH, NH3 and CH4.

ASW has been shown experimentally to be capable of incorporating and trapping other chemical components, e.g. H, CO, CO2, NH3, CH3OH,14 and potentially more complex organic molecules.15 The trapping efficiency is dependent on the ice morphology; for instance, porous amorphous water ice has been shown to absorb 20–50 times more gas (such as N2, CO, Ar, O2, and CH4) than compact water ice.16 The rates of chemical reactions involving the adsorbed molecules may also be affected by the underlying structure. For example, hydrogenation of both CO and C2H4 occurs more frequently because of an increased sticking coefficient of hydrogen atoms.17,18 These structural effects will be important to simulations of interstellar ices, in which chemical reactions are critical to both the composition and structure of the ice.19 It is therefore necessary that structural treatments become available for interstellar grain-surface kinetics simulations, and that these treatments also be consistent with the available experimental data.

Laboratory ice experiments are one means to gain an understanding of the physical structure and chemistry that can occur in astronomically-relevant ices (for further discussion on ASW surface formation, characteristics, and processes, see Hama and Watanabe 200820). However, it is yet unclear whether the ices that are grown under laboratory conditions are closely representative of interstellar ices. Various conditions must be considered when forming an interstellar ice analog in the laboratory, such as the deposition rate, temperature, deposition angle, and thickness of the deposited ice. All of these parameters can alter the ice morphology, and may play a significant role in interstellar ice structure. Typical laboratory experiments involve a high- or ultra-high vacuum system (UHV), i.e. 10−7 to 10−10 Torr, with temperatures usually ranging from ∼10 to 300 K. The typical gas density corresponding to a UHV system is <107 molecules cm−3 – higher than the values that obtain in dense interstellar clouds, which range from around 106 down to 104 molecules cm−3. Furthermore, the production of water ice under interstellar conditions typically occurs through the deposition of atomic precursors, which then form water through chemical reactions (although in certain types of interstellar object, such as protoplanetary disks, direct water deposition may indeed dominate). Due to the different gas densities and formation mechanisms operating in typical laboratory and interstellar regimes, to reach equivalent coverages or thicknesses in either case requires vastly different timescales (minutes versus many thousands of years).

Monte Carlo chemical kinetics models by Garrod (2013)19 and a laboratory study by Oba et al. (2009)21 found direct deposition should produce much greater porosity than the chemical formation of water ice. In the models, diffusion of oxygen atoms at 10 K, prior to their hydrogenation to water by atomic H, allowed microporous structures to be filled in more effectively over long interstellar timescales, whereas the water molecules directly accreted from the gas phase were essentially immobile, and therefore formed highly-porous structures. A prevalence of porous structures within interstellar ices will also significantly increase the time required for any arbitrary diffusing particle to visit all surface binding sites, increasing the meeting and reaction timescales in cases where only a few reactants are present on the surface at any moment.

The importance of the trapping of other molecules within interstellar water-ice mantles becomes yet more important during the process of star formation; during the so-called hot-core phase, the ices are warmed and sublimated, resulting in dramatic, detectable changes in gas-phase chemistry. More volatile species are released first, but much of this material may remain to higher temperatures as a direct result of trapping within the water-ice structure. These effects are demonstrated in temperature-programmed desorption (TPD) experiments for a range of molecules.14 Here, mixed molecular ices – formed through successive or concurrent low-temperature deposition of water and some other molecule – are heated linearly until the mixture has sublimated, while the rates of desorption are measured over time/temperature through mass spectrometry. Different components of the ice appear to desorb from the surface at different temperatures, depending on the strength of binding between the different molecular components and the local ice structure. For example, CO desorbing from a compact amorphous water ice will produce a different TPD spectrum than from a pure CO ice or from a porous amorphous water ice.22 Indeed, the desorption of CO from water ice demonstrates up to four distinct desorption events at distinct temperatures. Understanding this temperature-dependent behavior requires kinetics models that can treat structure explicitly. The accurate extrapolation of the laboratory data to interstellar timescales also requires such kinetics models. In both cases, the simulations must be efficient enough to treat the entire heating process kinetically, while also being accurate enough to simulate the ice structure that determines the kinetics of diffusion and desorption. The first step in modeling this process is therefore the reproduction of the porosities and structures produced in the laboratory for interstellar ice analogs at fixed, low temperatures.

A number of techniques, including X-ray diffraction, optical measurements of refractive index, infrared spectroscopy, and chemical adsorption, have been used to quantify the effects of experimental deposition parameters on ASW ice structure and density. Diffraction experiments are used to find distribution functions because atom arrangements yield unique diffraction patterns, which illuminate the underlying structure.23–25 In a chemical adsorption experiment, the ice is exposed to a gas species e.g. N2, CH4. The uptake of the molecules increases with porosity due to the increased surface area at the ice-vacuum interface. Infrared spectroscopy can also be used to qualitatively study the ice porosity via features that occur at the surface of the ice. In the case of ASW ice, dangling bond features from surface H and OH groups around 2.7 μm have been used to monitor the ice porosity; however, a number of studies have found these features may not be present even for a significantly porous ice.26–28 This technique has the advantage that it can be directly compared to infrared observations of astrophysical regions.29

The formation of laboratory water ices has been treated computationally using a number of simulation methods, such as ballistic, molecular dynamics, and kinetic Monte Carlo models.2,30–33 Kimmel et al. (2001) utilized a simple ballistic model where deposition angle was considered. Temperature was not explicitly treated, but the incoming particle was allowed to hop a fixed number of times before the next particle was deposited.34 Treating the temperature and diffusion appears to be important for correctly characterizing structure. Another condition that should be considered is the kinetic energy of each incoming particle, because it may be converted to translation energy. This relaxation process has been modeled extensively by Yang et al. (2001) and Zhou et al. (1997) for high-energy collisions (1 to 20 eV) for thin nickel films deposited from ionized vapor.35,36 In the case of laboratory ice experiments, as well as under interstellar conditions of interest, typical energies will be much lower (well under 1 eV). However, the adsorption of a water molecule into a local surface potential minimum will imbue that particle with some further kinetic energy as it is accelerated toward the surface at short range. This combined kinetic energy may allow the newly-adsorbed particle to undergo one or more surface diffusion events before thermalizing with the surface. Such occurrences have not yet been considered in kinetic simulations of water-ice production. Some recent work has investigated non-thermal hopping, but the source of energy was the exothermicity of reactions.30,31 The work described here investigates non-thermal hopping activated by the energy gained from deposition.

In this paper we investigate the formation of porous structures within ices formed in a laboratory setting and seek to reproduce the trends in porosity associated with the main experimental parameters during the ice-formation process. We modify the interstellar off-lattice kinetic Monte Carlo model MIMICK (Model for Interstellar Monte Carlo Ice Chemical Kinetics),19 a version of which was first presented by Garrod (2013), to allow application to laboratory conditions; here, periodic boundary conditions are applied, allow a small section of laboratory ice to be modeled self-consistently. The use of this off-lattice model allows a fully three-dimensional picture of amorphous ice structure to be traced, as it is formed kinetically through deposition and diffusion, and allows the conditions used in published experimental studies to be replicated. The tuning of model parameters, achieved through the successful reproduction of experimental results, allows ice production under the more extreme conditions of the interstellar medium – currently inaccessible in the laboratory – to be explored through simulations, and the fidelity of interstellar ice analogues to real interstellar ices to be assessed. The full exploitation of experimental ice data to improve our understanding of interstellar ice formation makes kinetics models like the one presented here essential.

2 Aims and methods

We use an off-lattice microscopic kinetic Monte Carlo model called MIMICK to simulate laboratory deposition of amorphous water ice. With this model, we examine how amorphous H2O ice structure is affected by various deposition parameters, including the ice temperature, the rate of deposition, and the angle at which the incoming gas molecules approach the surface. Past computational studies have varied the deposition rate, matching the rates used in laboratory experiments, with the sole intention of replicating laboratory results and have not been extended to replicate the ISM. The presented work is the first application of a model to the direct comparison of laboratory and interstellar ices, using respective deposition parameters.

The main feature of the experimental ices that we aim to reproduce is the porosity. The advantage of using MIMICK is that the positions of all particles are known at all times, so that porosity may be explicitly traced, and can be directly observed using three-dimensional representations of the simulated ice. The degree of porosity produced in the models can be quantified in various ways, but only methods that can be compared directly with laboratory values are useful in constraining the model parameters. We therefore use the ice density as a proxy for porosity and the structure, as is done by various experimental groups. The densities of simulated ices are relatively simple to calculate (see Section 2.2).

The predecessor of this model is discussed in detail by Garrod (2013),19 while the latest version is described by Garrod (in prep.). The code has been adjusted to use a surface with periodic boundary conditions, and now includes a relaxation process discussed later in this section. Important features of the code used here are described below.

The MIMICK model is based around the need to simulate surface diffusion kinetics, which is governed by rates of hopping between surface potential wells. The model therefore traces the positions of entities that could plausibly engage in surface diffusion – either atoms or molecules (including radicals). Thus, where a molecule exists, it is treated as a self-contained entity, rather than separate atoms. Here, only water molecules are considered (excluding the immobile atoms that comprise the substrate surface). All binding between entities on the surface is of a van der Waals-type nature; no chemical bonds are considered (except implicitly, within molecules). While the model is capable of tracing chemical reactions between surface species, the choice of system considered here (water only) leaves this feature redundant. Furthermore, all physical interactions between surface entities are described by isotropic pairwise potentials, making all surface species effectively spherical. We therefore refer to all independent entities as “particles”, whether atom or molecule.

Fig. 1 shows a simple representation of the flow of computational processes used in MIMICK. The main loop of the model involves: the evaluation of rates for every possible physical/chemical process that could happen; the random choice of a single process to occur as the next event (with probabilities weighted according to the magnitudes of each rate); the carrying out of the chosen process; and the advance of the time parameter. All physical/chemical events occur sequentially, with rates dependent on the outcomes of previous events. The model ends once the final time or the desired coverage has been reached; the desired coverage is achieved when a number of water molecules has been deposited that is equal to the desired number of layers times the number of surface sites on the bare substrate (see Section 2.2). The model traces the positions of all particles on the surface at all times, as various deposition, diffusion and desorption processes occur.

image file: c7cp05966c-f1.tif
Fig. 1 Basic flow chart of the operations carried out within the kinetic Monte Carlo code MIMICK.

Because MIMICK is an off-lattice model, there is no pre-determined grid of allowed positions for particles; all particle positions are determined through the minimization of the sum of pairwise potentials that they experience due to interactions with other surface particles. Since the constant global minimization of potentials would be prohibitively expensive in computational time, only the position of a particle that has just accreted or diffused is optimized. Low temperature ices with amorphous structures are best treated using an off-lattice model, where disorder is inherent. The off-lattice treatment allows amorphous surface structures to form containing pores.

The model starts with certain pre-defined initial conditions, which consist of a bare substrate surface at some temperature. This surface is constructed of two layers of atoms fixed in a body-centered cubic arrangement; the substrate surface atoms are not allowed to move during the simulation (although water molecules are), but each one produces an isotropic pairwise potential with any deposited particle that comes close enough to interact. The highly-ordered substrate surface employed here is not very representative of interstellar grain surfaces, which are expected to have irregular or fractal geometries, but it is a good basis to compare to laboratory experiments, which often use regular surfaces such as well-defined metal or ionic crystalline surfaces. Periodic boundary conditions in the lateral directions are utilized, to allow a small section of surface to represent a large experimental surface. A series of surface sizes ranging from 80 × 80 to 160 × 160 atoms is used to ensure that structural features remain smaller than the scale size of the ice.

The only process which may occur at the very beginning of the simulation is the deposition of a water particle from the gas phase. All deposition events consist of the placement of a new particle at a starting position at a fixed height far above the surface, whose lateral position is randomized in the x and y directions. If so-called background deposition is assumed, then the velocity vector of the incoming particle is randomized over the solid angle of a hemisphere centered on the starting position. In the case of fixed-angle deposition, the velocity vector is pre-determined and uniform for all deposition events. The trajectory is then gradually advanced until a surface particle (either a water molecule or a substrate atom) is encountered. The progression of the trajectory is then halted, and the position of the incoming particle is optimized according to local pairwise interactions with surface particles.

The positions of surface particles, following deposition or diffusion, are determined by calculating local minima via a Lennard-Jones potential, i.e.

image file: c7cp05966c-t1.tif(1)
where r is the distance between interacting particles, ε is the depth of the potential well, and σ is the optimal distance between particles. In practice, the above expression is adjusted, such that the potential gives a zero-value at a distance of 2.5σ or beyond. This is achieved by adding a small modifier to the expression in square brackets, then applying a small scaling factor to ensure that the potential minimum is fixed at ε. The standard expression would return a value around 0.8% ε at the truncation distance. The truncation of the potential simplifies the calculations, and is acceptable given that long-range interactions are of limited interest for the amorphous ices that are simulated here.

All chemical species within the model are treated as isotropic spheres with radius initially set to 1.6 Å, based on the van der Waals radius of a water molecule, following Garrod (2013);19 in eqn (1), σ is therefore set to a value of 3.2 Å. This value is further tuned in Section 3.2, to ensure appropriate ice densities.

The kinetic rates, such as diffusion and desorption, and how they are determined are briefly described in this paper, but a full description is provided by Garrod (2013)19 and Garrod (in prep.). The accretion/deposition rate is treated as a free parameter and is changed based on the desired conditions (101 to 106 molecules cm−2 s−1 for the ISM and protoplanetary disks, respectively and 1013 to 1016 molecules cm−2 s−1 for laboratory conditions).

The general temperature-dependent rate for thermal desorption or diffusion of a surface particle i is shown in eqn (2), where νi is the characteristic vibrational frequency of the molecule, Ei,n is the energy or barrier for process n, and T is the surface temperature.

Ri,n = νi[thin space (1/6-em)]exp(−Ei,n/T)(2)

The characteristic frequency for each individual molecule is determined using a one-dimensional harmonic oscillator relation37 based on the desorption energy of the molecule. The desorption energy (K) is simply the sum of all pairwise interactions with nearby particles (r ≤ 2.5σ). This is a departure from the method used by Garrod (2013),19 in which only contributions from contiguous particles were considered, which resulted in binding energies that were exactly proportional to the number of nearest neighbours. Here, weaker, somewhat longer ranged interactions are also considered, which produces a less discrete distribution of binding energies across all particles on the surface.

In the past version of the model,19 a successful desorption event led to the immediate removal of the desorbing particle from the surface. Here, a desorbing particle is assigned an outward trajectory vector, which is calculated simply as the vector sum of pairwise bond directions (unit vectors) for all immediately contiguous binding partners. pairwise binding strengths are not considered in the calculation of this direction vector. This vector sum therefore represents the mean geometric normal to the local surface. The particle trajectory is then traced outward until it exits the upper boundary of the simulation, or it collides with another surface particle, at which point it is re-adsorbed onto the surface, following the same procedure as regular deposition. This “re-accretion” process is important in cases where a particle desorbs from inside a pore.

Calculation of the barriers against diffusion is more complicated than the determination of the desorption energy, as these barriers are dependent on the spatial arrangement of the surrounding particles. Due to the amorphous nature of the ice that is produced in these simulations, not all diffusion pathways have the same barrier. To determine the number and directions of diffusion pathways, it is assumed that diffusion always occurs through local saddle points, which may be found by allowing the diffusing particle to rotate around an axis that joins a pair of surface particles that are contiguous with the diffuser (although not necessarily with each other). Simple geometric checks are also carried out to ensure that such pathways are physically allowed, and not immediately blocked by another particle. As explained by Garrod (2013),19 a plane may be defined that intersects with the initial position of the diffusing particle and the positions of the pair of particles that define the saddle-point. The positions of all other particles immediately contiguous with the diffuser are then assessed in relation to this plane; if all these surrounding particles fall on only one side of the plane or the other, then diffusion around the specific “saddle-point pair” is allowed. If one or more contiguous particles is on the other side of the plane from the others, then diffusion around the pair in question is deemed to be obstructed, and is not allowed. Importantly, no calculation of a post-diffusion position is carried out, until after such a diffusion event has been selected as the next process to occur in the simulation, which saves significant calculation time. For each diffusion pathway available, the energy barrier is determined by summing the energies that would be required either to break or stretch each pairwise surface bond to the furthest distance it could reach as a result of rotation around the axis, ignoring those bonds that would be made stronger by diffusion in the given direction. This method diverges from that described by Garrod (2013),19 which required the outright breaking of all bonds except those of the saddle-point pair. As a result, the surface diffusion barriers take on a smoother range of values, as well as generally being lower.

With the diffusion barriers for each direction determined, if diffusion by some particle or other is chosen as the next event in the sequence, the diffusing particle is rotated around the relevant saddle-point pair until it collides with one or more particles on the other side. Its position is then optimized according to its new location, in the same way as with deposition.

Only two pairwise potentials need to be defined in the models presented here. We initially assume ε(H2O–surface) = 550 K, with ε(H2O–H2O) taking a range of values from 700–1200 K, intended to take some account of the strong pairwise interaction produced by hydrogen bonding. This ε(H2O–H2O) value is varied in Section 3.2 (along with the inter-molecular distance, σ) to obtain the best fit to an experimentally-obtained ice density. The strength of the pairwise potential between water and a substrate atom is chosen so that surface binding is strong enough to ensure that ices can actually be formed at the 130 K substrate temperatures at which some experimental ices are formed in the data set of Brown et al. (1996).38 Test models suggest that as long as water molecules are able to stick to the surface at the given temperature and deposition rate, then the choice of surface–water pairwise potential has little effect on the results for the relatively thick ices modeled here.

The use of isotropic potentials, and thus the use of spherical particles to represent water molecules in particular, is simplistic when compared with more careful treatments of static ice structure. Low-density amorphous ice is typically found to exhibit broadly tetrahedral, albeit distorted, structures first reported by Narten et al. (1976) using X-ray diffraction to calculate radial distribution functions,23 and demonstrated by many others.24,25

While these structures may be approximated by the use of a pre-defined grid,28 the advantage of the off-lattice approach taken here is that it may be extended to various mixed molecular ices, for which a rigid grid may be a poor representation. An off-lattice method using directional potentials could in theory be employed, but this would introduce a significantly greater burden to the search for surface potential minima and the determination of diffusion pathways. The importance of surface diffusion kinetics in determining the larger-scale structure in the models presented here means that the above-mentioned simplifications are essential to achieve computational solutions on manageable timescales. (The CPU time required to run the models is discussed further in Section 3.1). Alternative non-directional potentials (versus Lennard Jones) for water molecules in particular are, however, a possible avenue for application in the model in the more immediate future.

The pairwise interaction strengths of 1200 K or less used in the models presented here are significantly lower than the dissociation energy found using more rigorous potentials and directional hydrogen bonding.39,40 However, the use of higher values in the present model produces ices that are far too porous and show little temperature dependence. This is to be expected, as the typical coordination number for individual bulk-ice water molecules in the model is around 8–10, rather than the 4 of tetrahedral water. The pairwise potentials used in this model are in fact around a factor of two to three smaller than typical evaluations in the literature.39,40

The lower pairwise interactions employed in this model are therefore not expected to be independently consistent with more careful treatments of water-ice structure. Indeed, an efficient chemical kinetics model that can address the requirements outlined above, while also matching the accuracy of the structures produced by the much more intensive molecular dynamics (MD) method, is out of reach for the time being. Instead, the model parameters (σ and ε) are tuned to produce appropriate densities for low-density amorphous water. The determination of mean molecular desorption energies from the model provides another measure of the fidelity of the model to macroscopic properties.

Future applications of the model will use the results, parameters, and porous structures determined here to model systems in which ice-surface reaction kinetics are involved, and to reproduce TPD experiments with mixed ices in which porosity is important; such applications will place even greater demands on the kinetics treatment used in the MIMICK model.

2.1 Deposition parameters

Experimental ice thicknesses are typically defined by the number of monolayers (ML). In these models, we define a monolayer as a single sheet of molecules (without pores), that will fit precisely on top of the bare substrate. Hence, simulations with substrate surface sizes of 80 × 80 and 160 × 160 have 160[thin space (1/6-em)]000 and 640[thin space (1/6-em)]000 molecules, respectively, for a 25 ML thick water ice. Typical dense cloud ices are ∼100 ML in thickness, on the order of tens of nanometers. Experimentally, ices have been studied in the range from tens of nm to μm thick. However, ices of thicknesses around 1–5 μm can be prone to spontaneous cracking, accompanied by the production of very large pockets that are phenomenologically different from the porosity that we aim to model here. Furthermore, MIMICK accepts only a uniform temperature value throughout the ice, and therefore does not take into account any variation in temperature associated with distance from the cooled surface, which may be important for thicker ices.34,38,41–43 Since results for thinner ices are also available in the literature, experiments closest to astronomical thicknesses will be used for comparison with the simulations.

Preliminary models indicate that surface sizes significantly lower than our smallest surface size of 80 × 80 substrate atoms (e.g. 10 × 10), produce surface and bulk-ice structures that are comparable in scale to the size of the periodically-repeated cell in which they form, raising the possibility that, in such a regime, the cell size itself may control the size of the simulated structures. Similar effects can also become apparent where ices are allowed to form that are significant thicker (taller) than the substrate is wide. To avoid such effects, we use at least an 80 × 80 ice, using the 160 × 160 ice where model run-times allow. The choice of using 25 ML as our standard thickness avoids this problem in the vertical direction. Surface effects and computational time are discussed further in Section 3.1.

It should be noted that the experimental data sets that we compare our results with most frequently in this work, which come from Brown et al. (1996)38 and Berland (1995),42 were obtained for ices ranging from ∼700 ML to 5 ML, respectively, with densities measured in each dataset that are reasonably consistent with each other, showing no significant divergence. However, the Berland et al. data, while closer in ice thickness to the 25 ML ices that we simulate, cover only a few temperature values, while the Brown data are well-sampled from 15–125 K. The Berland et al. ices may also be thin enough that the ice structures formed local to the substrate surface could influence the density values, possibly diverging from the values obtained in the bulk water ice. For these reasons, we mostly constrain our comparisons to the Brown et al. (1996)38 data.

In the models presented here, ices are deposited onto the surface at temperatures ranging from 10 to 130 K. Deposition rates are varied from 101 to 1016 molecules cm−2 s−1 (equivalent to the number of molecules delivered to the surface at 10−20 to 10−5 Torr). Additionally, a deposition rate of 106 molecules cm−2 s−1 was used to replicate conditions pertaining to protoplanetary disks.

Two different deposition methods are used to reproduce experimental treatments: (1) deposition angles are randomly generated individually for all molecules, to replicate experimental background deposition, and (2) incoming particles are all deposited at the same angle with respect to the surface (ranging from 0° to 75°). A number of ice thicknesses were investigated (10 ML to 100 ML).

2.2 Calculating density

The method used to calculate the density of the ices produced by the models must take into account the effect of porosity; it must therefore be an average over the entire bulk ice, rather than be based on small sections, whose densities may vary widely at a local level. Since the number of water molecules is known, the calculation of density reduces to a calculation of the volume of the ice produced. To calculate this volume, we divide the ice into vertical columns, each of base size σ × σ, and measure the height of the highest water molecule within it, to determine a volume. We then sum these volumes, and divide the total mass of water molecules by this value. Densities are reported in g cm−3.

2.3 Relaxation process

The basic MIMICK treatment considers only thermal diffusion. Results presented later, however, indicate that a non-thermal mechanism is also necessary. The relaxation process we introduce here accounts for the energy of an incoming particle, which can potentially allow it to hop (or indeed react, if reactive species are present) on the surface. The energy gained is a combination of the thermal kinetic energy that the particle brings from the gas-phase, and the energy gained when the incoming particle accelerates into a surface potential well at short range; the maximum energy available from this latter process is the total interaction energy of the particle with the surface, although at least some of this will be lost the inelastic collision of the particle with the surface.

In typical experiments, the gas from which the deposited particles derive is held at room temperature; we assume that it will maintain the same energy distribution when it is introduced to the vacuum chamber. In that case, the average kinetic energy is [/]kT. Expressing this energy in units of Kelvin, as is typical for astrochemical models, gives an average energy of 450 K, for a gas temperature of 300 K. It is noteworthy that such an energy is likely on the order of just 10% of the binding energy of a typical water molecule, or likely no more than one third the energy required to overcome a surface diffusion barrier for water. The kinetic energy gained due to acceleration into the surface potential well should therefore dominate the non-thermal energy available for immediate diffusion on the surface.

While detailed calculations for the efficiency of diffusion of a newly-adsorbed molecule could be obtained through a molecular dynamics treatment, this would be a major endeavour, and such calculations could not be directly incorporated into a kinetics model, due to the much greater timescales over which the kinetics models operate. The average efficiencies derived from separate MD calculations could potentially be incorporated into the present models, but even the most up-to-date MD treatments of excited molecules on water surfaces deal only with crystalline water ice,44 rather than the amorphous structure of interest here. We instead seek to implement a more empirical treatment to gauge the effects of such a process, and therefore to determine the efficiency required to reproduce experimental results using our kinetics model.

We use a Maxwellian distribution to randomly select the kinetic energy of each incoming gas-phase water molecule, KEinit. Upon adsorption to the surface (into a potential well whose strength is determined according to local conditions), a kinetic energy equal to the surface potential is added to this value. An efficiency factor (between 0 and 1), α, which is a free parameter, is then applied to this total energy, which sets the initial energy available to the molecule that can be applied to hopping over a surface diffusion barrier, i.e.

Einit = α(KEinit + PEinit)(3)

If the newly-adsorbed particle has a diffusion pathway accessible that has an energy barrier, Edif, less than or equal to the available energy, Einit, a diffusive hop spontaneously occurs. If more than one diffusion pathway is accessible, each is assigned an equal probability of occurring. Whichever hop is chosen, an energy equal to the operative diffusion barrier (corresponding to the energy cost in broken/stretched bonds) is subtracted from the available energy when the hop occurs. However, more energy is then added to the total, corresponding to the energy gained as a result of settling into the new binding site, ΔE, so that the available energy after one non-thermal hop is given by:

E1 = (EinitEdif + ΔE)β(4)
where E1 is the energy available for further hopping after 1 non-thermal hop; β is a second efficiency factor (also between 0 and 1) that accounts for inelastic losses associated with the particle's motion into the new binding site. The energy gained by moving into the new binding site, after factoring in the loss associated with leaving the old one, is defined by:
ΔE = PE1 − (PEinitEdif)(5)
where is the total binding energy associated with the new potential well, PEinit is the binding energy of the first binding site, and Edif is the diffusion barrier. If enough energy is available after the first non-thermal hop, another may occur, assuming one or more of the new diffusion barriers are low enough to be overcome with this energy. In such cases, the same value of β is assumed. If, at any stage, the available non-thermal energy is insufficient to overcome any diffusion barrier currently open to the particle, the period of non-thermal diffusion is over. The particle can then take part in thermal diffusion normally.

By defining global values for α and β, we allow average efficiencies to be assumed for the retention of kinetic energy following collisions relating to the initial impact of adsorption and subsequent diffusive hops to new binding sites. This avoids making detailed calculations of the energy losses and gains associated with each separate process, while retaining non-uniform behavior for each event that is dependent on local binding and diffusion conditions.

3 Results and discussion

We first discuss the computing requirements of the model. We then present results from preliminary models that allow the tuning of the model parameters, using experimental data in the regime where there is little/no temperature dependence for measured ice densities. We then present the results of this model as applied to fixed-angle deposition, and then to background-deposited ices at various temperatures.

3.1 Computational considerations

The MIMICK model runs in a single thread for each model setup; the advantages of parallelization are limited for this type of model, due to the sequential nature of the calculations. The code can, however, be optimized by storing previously-calculated particle positions, in the case where uninterrupted chains of surface hops occur for individual molecules; diffusing particles often return to the same sites multiple times. This chain information is saved even if the chain is broken by multiple particles taking turns to move on the surface, so long as they do not physically interact. This approach is essential to allow the code to be run on manageable timescales, but necessarily demands significant RAM allocation (averaging around 1–2 kB per particle in the model). However, in all the models presented here, the required RAM allocation is easily handled, and the main limitation on the size of model parameters is the run-time.

The bulk of the calculations carried out in the code relate to thermal hopping between binding sites. Therefore, the main factors that influence the CPU time of the model are the number of molecules being deposited (dependent on surface size and coverage), the rate at which the ice is deposited, the surface temperature, and whether the relaxation treatment is implemented. Higher dust temperatures and lower deposition rates will raise the time required for each run, due to the larger number of hops occurring prior to model completion. The relaxation process adds extra hops for each deposition event. The time required to run a model ranges from as little as 10 minutes (for a model at 10 K with a deposition rate of 1013, an 80 × 80 surface with maximum ice coverage of 25 ML, and no relaxation calculations) to around a week for high surface temperatures combined with the largest surface or low deposition rates. Given the dependence of the computational time on the number of molecules, we use the smallest surface size possible (for a given target ice thickness), while maintaining a large enough surface to prevent surface inhomogeneities from dominating the density value, and to avoid pore-scale issues (see Section 2).

3.2 Tuning of pairwise Interactions

A major influence on the simulated ice density is the strength of binding (ε) between each water molecule and its neighbours, which controls both the total binding energy of each molecule and its diffusion barriers. The general approach in the MIMICK model is to use non-directional potentials for all atoms and molecules, which in interstellar applications could include on the order of 100 species. In the models presented here, we require only two pairwise potentials: the surface–H2O interaction and the H2O–H2O interaction. The stronger the interaction, the higher the temperature must be to overcome the resultant diffusion barriers. When the value is low, the molecules have weaker local diffusion barriers, allowing them to diffuse more easily to find a deeper potential well, which produces a generally smoother ice. Conversely, with high values (assuming the same ice temperature), the strength of the H2O–H2O interaction can significantly slow – or even prohibit – diffusion, resulting in more irregular structures on large and small scales.

In order to achieve the appropriate water density for comparison with experimental values over a broad range of temperatures, we first tune the model to reproduce certain key values. The experimentally-determined densities of Brown et al. (1996)36 are plotted (along with test-model data) in Fig. 2, over the 80–130 K range; the full dataset extends down to 20 K. It may be seen that at temperatures around 110 K and higher, the experimental density values reach a plateau at around 0.87–0.90 g cm−3. We therefore assume, for the purposes of our parameter fitting, that thermal diffusion is efficient enough in this range of temperatures to produce the maximum surface rearrangement possible, thus achieving the densest structure available for low-density amorphous water. At lower temperatures, the thermal diffusion of water molecules is less able to rearrange the surface to reach this maximum density.

image file: c7cp05966c-f2.tif
Fig. 2 Densities of ices produced using pairwise interactions of 700, 850, and 900 K. Re-scaled densities are also shown, using scaling factors for each model determined by best-fit calculations to the experimental data of Brown et al. (1996).38

To find the best set of ε and σ parameters for our models, we therefore seek firstly a water–water pairwise potential that produces ice densities that achieve plateau behaviour beginning at around 110 K (regardless of the absolute value of the density achieved in the model). Fig. 2 shows the results from several models in the ε = 700–900 K range that achieve similar plateau behavior in the correct temperature range. A range of models was in fact run between 700–1200 K, for temperatures ranging from 80–130 K, in 10 K increments, with water deposited at random angles at a rate of 1013 molecules cm−2 s−1. High pairwise potential models tend to produce plateau behavior only at very high temperatures; lower potentials reach a plateau at lower temperatures.

We take the density results from these models and scale them with a uniform scaling factor, determined uniquely for each model (and representing divergence from the initial σ value of 3.2 Å). This factor is in each case calculated by minimizing the chi-squared function below to achieve the best match with the experimental data for each model.

image file: c7cp05966c-t2.tif(6)

We use the formulation in eqn (6) (typically used in Pearson's Chi-squared test) as a means to measure the goodness of fit. Index i indicates the temperature of interest, Oi is the observed density at a particular temperature and Ei is the corresponding density measured from Brown et al. (1996).38 The fitting uses the data in the 80–130 K range, corresponding not only to the plateau region, but also to somewhat lower temperatures, at which it is assumed that the density is still determined by purely thermal diffusion of water molecules (as is borne out by later models). The final best-fit values of chi-squared (per eqn (6)) determine not only the best-fit scaling factor for each model, but also the best-fit model out of all of those tested, across all values of ε. The best-fit values are listed in Table 1, while the unscaled and scaled density values are plotted over the experimental values in Fig. 2. Although the ε = 1200 K model shows moderately good agreement, with no scaling factor required, no plateau-like behaviour is in fact observed in those models; densities continue to grow for increasing temperatures. The model with ε(H2O–H2O) = 900 K did not achieve the plateau at the highest temperature desired. The overall best match to the experimental data was initially found to correspond to a value ε(H2O–H2O) = 800 K. An additional model was run using ε = 850 K, using a scaling factor of 0.735, which provided an even better match. We adopt the ε = 850 K value for all other models presented in this paper. Based on the required density-scaling factor for this model to match observations, we tune the value of the (one-dimensional) σ parameter from the initial value of 3.2 Å, to a value σ = 3.546 Å, so that the model results match the experimental values without a density-scaling factor.

Table 1 Best-fit density-scaling factors and associated values of the goodness-of-fit parameter (chi-squared) for models using various values of the pairwise interaction for H2O–H2O. The asterisk denotes the overall best-fitting pairwise potential, which was chosen for use in all subsequent models
H2O–H2O (K) Scaling factor χ 2
700 0.65 0.0034
800 0.71 0.0023
850* 0.73 0.0016
900 0.80 0.0041
1000 0.79 0.016
1100 0.97 0.0068
1200 1.00 0.0033

The adopted pairwise potential of 850 K results in an average total binding energy for surface molecules of ∼6700 K, based on the calculated averages at the end of the model runs at low temperatures. The values increase up to ∼7300 K for surface temperatures of 130 K. These values fall toward the upper end of typical ranges of binding energies for water in astrochemical models. However, binding energies for individual water molecules in the present MIMICK models can range to as low as ∼3000 K, making such molecules more likely to desorb. Since desorption of an individual molecule is likely to be preceded by multiple diffusion processes, the effective desorption energy may be lower than the average values quoted above, due to molecules finding and then preferentially desorbing from weak binding sites. To determine this effective binding energy, a full TPD model would need to be run. The value of 4800 K determined by Penteado et al. (2017)48 indeed falls between our lower value of 3000 K and the mean value of 6700 K for low-temperature deposition. A full analysis of TPD desorption characteristics is beyond the scope of the present project, but is of interest for future work.

3.3 Uniform incident angle

In order to test the effects of ice deposition at fixed incident angles, a number of ice simulations were run, using various angles of deposition, over a range of temperatures and deposition rates. Fig. 3 shows a side view of each of the resultant ices produced for deposition rates 101 and 1013 molecules cm−2 s−1, at temperatures of 10 and 30 K, with angles varied from 0 to 75°. The images are produced using the freeware ray-tracing software, POV-ray. While water is treated in MIMICK as a sphere, using non-directional potentials, it is nevertheless depicted in these images as a molecule, following the treatment by Garrod (2013),19 who chose this representation to aid identification of different species. The porous structure of the ices is clearly visible, demonstrating the effect of shadowing; voids are present between the column-like structures of H2O, especially for large angles of incidence. Ices deposited normal to the surface show gaps still present, but they appear smaller and less tubular than for high angles. Similar structures can be seen in Kimmel et al. (2001)34 and Dohnolek et al. (2003)33 where a ballistic model was used. Their model, while demonstrating angular dependence, did not investigate the effect of deposition temperature. An arbitrary number of hops or diffusion events were allowed to occur for each molecule, but the structure does not reflect the surface temperature. This is a critical component when considering an ice that may warm up and experience structural changes.
image file: c7cp05966c-f3.tif
Fig. 3 Simulated ices, shown in side-view, produced for deposition angles ranging from 0° to 75°. Results for two temperatures are shown (10 K and 30 K), divided into high and low deposition rates (101 and 1013 molecules cm−2 s−1). Columns of ice with interstitial voids are formed, whose angle traces the deposition angle. Greater porosity is produced for ices formed at lower temperature, higher deposition rate, and greater angle. Images were generated using the ray-tracing software, POV-ray.

These general trends in structure that are visible in Fig. 3 are also reflected in the density behaviour. Fig. 4 plots the simulated ice density against temperature for several incident angles (shown in blue), for temperatures ranging from 10 to 115 K; only the high deposition rate results are shown. Background deposition (i.e. random angle deposition), discussed in more detail later, is shown in green, and exhibits behaviour comparable to that of ices grown at a 60° angle. Higher temperatures during deposition produce higher densities, although at high deposition rates, little variation is observed below temperatures of 35 K. Each of the deposition angles shown produce density variations of a factor of ∼2, over the 10–115 K temperature range.

image file: c7cp05966c-f4.tif
Fig. 4 Simulated ice densities as a function of temperature for ices deposited at a rate of 1013 molecules cm−2 s−1. Ices grown at incident angles ranging from 0° to 75° are marked in blue; the results for background deposition are also plotted. The ice densities for background deposition are comparable to those deposited at a fixed angle of ∼60°.

Increased temperatures allow more loosely-bound water molecules to diffuse into nearby potential minima that provide stronger binding. The stronger sites will tend to be those with pronounced inward surface curvature at a microscopic level, as these will allow a water molecule to interact with more H2O neighbours. This effect may be seen visually in Fig. 3, for the low deposition-rate ices; higher temperatures produce smoother ices, indicating that microscopic-scale roughness is being smoothed away, while some larger-scale structures still remain.

The same structural trend seen in Fig. 4 was found by Stevenson et al. (1999) who measure porosity by N2 adsorption to water ices formed with various deposition angles. They found that N2 adsorption was greatest at low temperatures and angles, with the exception of a normal incident angle, for which a consistently low degree of porosity, with no temperature dependence, was observed.45 This particular result therefore disagrees with our simulation results, which show a clear temperature–density dependence for all angles. Unfortunately, it is difficult to compare directly the experimental N2 adsorption results with the simulated densities, as there is no means of direct conversion immediately apparent. It is possible that the degree of adsorption may also incorporate some measure of the connectivity of pores within the ices. The more irregular shapes of the pores produced by the simulations at normal angles are arguably consistent with a picture in which the efficient uptake of gas-phase N2 would be limited only to those pores closest to the outer surface. In this scenario, the pores throughout the ice would become isolated from each other, either during their formation, or as a result of blockage of any narrow connective passages by the adsorption of the N2 itself. It is unclear without more detailed modeling whether the latter effect might be possible for the ices shown in Fig. 3. The direct simulation of experimental adsorption of N2 on porous ices is a topic for future study.

3.4 Background deposition

Direct determinations of the densities of amorphous ices produced through background deposition are available in the literature. It is therefore valuable to attempt to reproduce these experiments using MIMICK, as we may compare both the density trends and the absolute values, in order to constrain the models and to shed light on the processes occurring as the ices are formed.

Background deposition was simulated for a range of deposition rates (101, 106, 1013, and 1016 molecules cm−2 s−1) and temperatures (10–120 K). The images and results are shown for ices with a thickness of 25 ML. Fig. 5 shows side-views of a selection of the resulting simulated ices. Sticking was not achieved for models with a deposition rate of 101 molecules cm−2 s−1 at temperatures above 40 K. The angle of incidence for each incoming particle is randomly generated, mimicking the likely behaviour of both experimental and interstellar ice deposition.

image file: c7cp05966c-f5.tif
Fig. 5 Simulated ices formed through background deposition at various rates, for a range of surface temperatures with immediate relaxation. Molecules deposited at the lowest rate are unable to form an ice at high temperature. Additional images are shown for temperatures 15 and 25 K for this deposition rate. The ice grown at 35 K at this rate has a density of 0.786 g cm−3, similar to the ice grown under laboratory conditions at 85 K. Ice densities are greatest for higher temperatures and lower deposition rates.

The total fluence of deposited water molecules is held constant between simulations, so the number of molecules is approximately the same for each case and for each image shown in the figure. The ices may be seen to become significantly more compact as temperatures increase from 10–120 K. Large-scale irregular structures or protrusions are less frequent at higher temperatures, and less small-scale porosity is apparent; the porous structures are also more regular and somewhat wider, consistent with more effective thermal diffusion. Decreasing the deposition rate for simulations at sufficiently high temperature also produces more compact ices; however, at lower temperatures, diffusion is slow, and there is little difference in the compactness produced for the different deposition rates.

The densities calculated for all the simulated ices are plotted in the left panel of Fig. 6, shown as coloured lines. Experimental data from Brown et al. (1996; black triangles)38 Berland et al. (1995; black squares)42 are also shown. Dohnalek et al. (2003) shows the same trend from the Brown data, but is not shown here for clarity. A rate of 1013 molecules cm−2 s−1 in the model is the most comparable to the laboratory experiments of Brown et al. (1996)38 with estimated rates between 1013–1014 molecules cm−2 s−1.

image file: c7cp05966c-f6.tif
Fig. 6 Densities of amorphous water ices grown at various temperatures and deposition rates using randomly-generated angles for each molecule. The gradient shows an increasing α and β for a deposition rate of 1013 molecules cm−2 s−1 comparable to the laboratory data of Brown et al. (1996),38 shown in black.

The figure demonstrates more clearly the increase in density with surface temperature. The ices become smoother and thinner at higher temperatures, while increasing the deposition rate has the opposite effect. This effect has been reported by Berland et al. (1995),42 who showed using refractive index measurements that increasing the deposition rate causes the density of ASW ices to decrease. If the deposition rate is much faster than the rate of diffusion, the molecules do not have time to rearrange into lower energy sites before another molecule is deposited nearby, locking the structure more firmly in place. At higher temperatures, the model results are in good agreement with the experimental data, at the appropriate deposition rate of 1013 molecules cm−2 s−1; both the general trend and the absolute density values are well reproduced.

In the low-temperature models, diffusion essentially does not occur rapidly enough to compete with deposition, and all deposition rates of a magnitude achievable in the laboratory (i.e. the two higher values) produce ices with the same low density. This marks a clear divergence from the experimental results. In the lowest deposition-rate case, appropriate to interstellar conditions, densities more comparable with the high deposition-rate experimental values can be obtained, indicating that thermal diffusion has time to re-arrange the surface structure before it becomes locked in.

These simulations therefore suggest that where thermal diffusion is sufficiently active, the modeled ices are a good representation of experimental densities and/or structures, over a range of temperatures. The break-down in fidelity suggests that at low temperatures, a non-thermal process may take over that was not included in the basic version of MIMICK that has been employed up to this point.

3.5 Relaxation process

Here we include the relaxation process described in Section 2.3, whereby the sum of (i) the residual gas-phase kinetic energy of the incoming water molecule, and (ii) the energy it gains as it is accelerated into a surface potential well, is made accessible to the adsorbed molecule to overcome diffusion barriers. It may continue to hop over barriers to diffusion until insufficient excess energy remains, by which point it has relaxed, becoming fully thermalized with the surface. The models presented in Section 3.4 may be considered to assume an immediate thermalization of all deposited molecules.

The relaxation process has two parameters associated with it that can be thought of as efficiencies, α and β, corresponding respectively to the fraction of energy retained upon the initial deposition, and the fraction retained after settling into a new potential well as a result of diffusion. While all of the energies considered in these calculations are determined locally, according to the nature of the molecular constituents of the surface and their relative positions, the efficiencies α and β are held constant for all molecules, throughout each ice deposition simulation.

Values for both α and β ranging from 0.1 to 0.9 were tested in the models; the smaller 80 × 80 surface was used for consistency throughout the model grid, due to the longer run-times required for the most extreme α, β values, which scale roughly with surface size. The resulting densities for some of the simulations are compared with the laboratory data of Brown et al. (1996)38 in the right panel of Fig. 6. It may be seen from the figure that higher-end α-values (∼0.7–0.8) produce a very good match to the experimental results. Variations in α have a strong effect on the overall agreement; the most extreme model we tested (α = 0.9, β = 0.9) is dominated by the non-thermal diffusion produced by the gradual relaxation of the molecule on the surface, showing no significant temperature dependence at all. The corresponding results from Section 3.3 (α = 0, β = 0) are also shown on this plot.

Fig. 7 shows the temperature dependence for each model in the grid, with each panel corresponding to a different α-value. Datasets for the various β-values employed are indicated by the colour of the dots. Data from Brown et al. (1996)38 are marked in black as before. It is again apparent that the α = 0.7, 0.8, 0.9 values provide the closest match to the experimental data, although the α = 0.8 dataset shows the best match over all possible β-values.

image file: c7cp05966c-f7.tif
Fig. 7 Modeled ice densities as a function of temperature, for various values of the α and β parameters. Each box is a different α value from 0.1 to 0.9. The colours are β values from 0.1 to 0.9. The Brown et al. (1996)38 data is shown in black for comparison.

To determine the best match more precisely, a χ2 test was applied to the densities obtained with the various temperature runs for each model setup. Fig. 8 shows a map of the χ2 values with the corresponding α, β pairs.

image file: c7cp05966c-f8.tif
Fig. 8 Divergence of the simulation results from the experimental densities of Brown et al. (1996),38 as a function of the relaxation parameters α and β. The divergence from experimental results is quantified as a chi-squared value. Smaller divergence is indicated by lighter colours.

A number of models were very close in error. An α-value of 0.7 or 0.8 was the best fit; β-values of 0.7 were the best fit for α = 0.7. A β-value of 0.7 was the best for α = 0.8. Based on these outcomes, several extra models were run in intervals of 0.05 for α and β. The resulting best fit was α = 0.7, β = 0.75. The χ2 error was a quarter of the error associated with the α = 0.7, 0.8 and β = 0.7 models. For the remainder of the study, the [0.7, 0.75] values are selected as the best match.

Images of the ices employing our relaxation treatment are shown in Fig. 9; the ice produced by a model with immediate relaxation, i.e. α = 0 and β = 0, is shown on the left. Increasing α produces a significant smoothing on small scales, removing porosity, as well as an increase in uniformity on large scales, such that the ices appear flatter. The β-values shown are the most extreme; this parameter has a relatively strong effect on structure when α is adequately high. The reason for this is more apparent in Fig. 11, which shows the efficiency and distribution of hops for a range of α and β values at 30 K. For an α-value of 0.7, all β-values are mostly dominated by single-hop events, until extreme β values; at lower α values, essentially all particles hop once. Starting at intermediate β values, the number of hops switches from 1 hop to a distribution of hops. At extreme β values 4 to 5 hops are dominating the structure. A broadly similar outcome is found for larger α-values, except in the case of α = 0.9, with high β; here, for β = 0.6, the two-hop outcome is as probable as a single-hop event, while for β = 0.9 the most probable outcome is four consecutive hops. In this case, four or more hops will occur with a probability close to 50%.

image file: c7cp05966c-f9.tif
Fig. 9 Images of ices deposited using various α and β values at 30 K. The deposition rate was maintained at 1013 molecules cm−2 s−1. The density increases as both α and β are increased.

For individual α, β pairings, the probabilities for different numbers of hops per deposition event are quite stable over the entire temperature range tested. The results presented in Fig. 10 are therefore a good representation for all temperatures.

image file: c7cp05966c-f10.tif
Fig. 10 Stimulated ices, shown in side-view, with the best fit α and β values of 0.7 and 0.75, respectively. The temperatures are varied from 10 to 130 K with a constant deposition rate of 1013 molecules cm−2 s−1.

It should be noted that for α = 0.7 and β = 0.75 the probability of one hop is greater than ∼25%, and for two hops ∼51%. While it is likely that the particle hopping twice can move back to its initial position, this only occurs in about 17% of cases. When β is large, more hops for a single particle are able to occur, because less energy is lost every hop. The number of particles and hops is also affected by the surface temperature (not shown), because diffusion may also occur thermally, which allows the molecules to find and occupy stronger sites, smoothing the ice and reducing the fraction of strongly-binding sites available. This reduces the energy required to hop, making multiple non-thermal hops more likely.

Fig. 10 shows the structure of model ices produced at several temperatures, using the best-fit α and β values. The structure of the ice clearly changes between each temperature shown. At 10 K the ice has its highest level of porosity. As the temperature is increased, from left to right, the structure becomes increasingly smooth and the pores gradually disappear. At high temperature (130 K), the structure is completely smooth with no porosity.

3.6 Interstellar and protoplanetary disk conditions

The simulations presented here allow water ices to form through direct deposition of water molecules; however, in many regimes in the interstellar medium, including dark clouds, ices are expected to form through surface chemical processes. One regime in which direct deposition of water molecules onto dust grains is a plausible formation process is within protoplanetary disks, where water may be repetitively desorbed and re-deposited, due to variable dynamical and physical conditions.46,47 The simulations so far have included model runs using deposition rates of 101 molecules cm−2 s−1 (similar to that suggested by Cuppen and Herbst, 200713), which are appropriate for a typical dark cloud with a density of 104 cm−3.13 Densities in protoplanetary disks, however, may be significantly higher. To test the effects on a water ice formed through direct deposition in this regime, models were run with a deposition rate of 106 molecules cm−2 s−1, which corresponds to a gas density of ∼109 cm−3.

The results of these simulations for various temperatures are plotted in Fig. 12, with data from Brown et al. 199638 for comparison. We use the best-fit values, α = 0.7 and β = 0.75 for the relaxation parameters. The dark cloud and the protoplanetary disk rates are shown in green and blue, respectively. When considering the relaxation process, the deposition rate relating to a dark cloud yields ice densities higher than the experimental values, indicating lower porosities. The porosity could be considerably lower, however, if the ice is formed through surface reactions. There is reasonable agreement between the protoplanetary disk conditions and the experimental data in a temperature range of ∼30–40 K, so laboratory data may be directly applicable to regions in disks that are within that temperature range. For other temperature values, and for the dark-cloud regime in general, porosities produced in experimental studies are unlikely to be a close match to interstellar ice structures, if formed through direct deposition. Closer correspondence may be obtained by producing laboratory ices at somewhat higher temperatures than interstellar conditions would suggest.

To this end, we provide in Table 2 a list of density values for interstellar and protoplanetary conditions, calculated at specified model temperatures, accompanied by the temperatures at which those same densities could be achieved under laboratory conditions. The temperatures were found by using a linear fit for the Brown et al. data for the temperature range between 20 to 110 K, shown in Fig. 11. The temperatures with laboratory conditions are significantly higher than those from the models using dark cloud and protoplanetary disk conditions to acquire the same density.

Table 2 Ice densities grown with deposition rates representative of dark cloud and protoplanetary disk regions are shown with the corresponding temperature at which ices grown using laboratory conditions have the equivalent density
Model temp. (K) Dark cloud ρ (g cm−3) Temp. equivalent (K) Protoplanetary disk ρ (g cm−3) Temp. equivalent (K)
10 0.618 32 0.609 28
20 0.648 43 0.624 34
30 0.676 53 0.647 43
40 0.759 84 0.679 54
50 0.721 70
60 0.765 86
70 0.838 113

image file: c7cp05966c-f11.tif
Fig. 11 Probabilities of varying numbers of consecutive hops, obtained from the simulations for different α and β values at 35 K. The bar colours correspond to number of hops carried out in succession.

image file: c7cp05966c-f12.tif
Fig. 12 Densities of deposition rates corresponding to a protoplanetary disk (blue) and a dark cloud (green). Laboratory data from Brown et al. (1996)38 are shown in closed squares with corresponding error.

4 Conclusions

The simulations presented here demonstrate a good match with experimentally determined amorphous water-ice densities. While the basic models produce a poor fit to lower temperature results (<80 K), they indicate that the production and prevalence of porous structures is significantly controlled by non-thermal diffusion. The modification of the model to incorporate a non-thermal diffusion mechanism based on the kinetic energy gained through the adsorption process allows the laboratory data to be well reproduced over all measured temperatures.

The degree of porosity under various conditions may be determined through ice-density measurements, which can also be calculated for the simulated ices. However, the models also indicate that the physical scales associated with the porosity can change, with higher temperatures and lower deposition-rates producing smoother ice surfaces with wider pores.

Simulations of ices produced with fixed angles of deposition are more difficult to compare directly with laboratory data, due to the lack of direct density measurements, but the trends with temperature, angle and deposition rate are appropriately reproduced, with the exception of the results for deposition normal to the surface. A more detailed computational reproduction of all laboratory conditions, including the adsorption of porosity tracers like N2 will allow such discrepancies to be explored in their proper context.

We note that, while the experiments and simulations considered here involve pure water ices, interstellar ices consist of mixed water, CO and CO2, as well as other components. Direct deposition of water ice is also unlikely to be the dominant interstellar ice formation mechanism, and experimental deposition rates are typically larger than those that are expected in the interstellar medium. The use of detailed kinetic models of laboratory ice deposition are therefore valuable in allowing experimental results to be understood in terms of microscopic kinetic processes, allowing extrapolation to interstellar conditions using the models.

We present several specific conclusions from this work below:

• Non-thermal diffusion in the formation of laboratory ices appears essential to explain low-temperature porosity measurements.

• We suggest that the energy gained by the adsorption of a water molecule to the surface is sufficient to produce enough diffusion to match experimental ice densities. A combination of one or two hops for ∼70% of newly-adsorbed water molecules is adequate to correct for the failure of the model at low temperatures. The majority of hops produce a net displacement of two surface sites from the original adsorption site. Around three quarters of the adsorption potential must be converted to kinetic energy of the adsorbed molecule to produce this result.

• Application of the model to protoplanetary disk conditions (in which direct deposition of water ice onto dust grains may be the dominant mechanism) suggests that laboratory ice porosities may be higher than those in protoplanetary disks at some temperatures. Under dark interstellar cloud conditions, laboratory ices were found to be too porous, under all temperature conditions tested.

Conflicts of interest

There are no conflicts to declare.


We would like to thank the anonymous reviewers for the helpful comments and suggestions. ARC would also like to thank Niklaus Dollhopf for useful discussions. RTG thanks the NASA Astrophysics Research and Analysis program for funding through the Grant NNX15AG07G. The authors would like to thank the FYRE program and UVa for funding and a research forum.


  1. D. C. B. Whittet, W. A. Schutte, A. G. G. M. Tielens, A. C. A. Boogert, T. de Graauw, P. Ehrenfreund, P. A. Gerakines, F. P. Helmich, T. Prusti and E. F. van Dishoeck, Astron. Astrophys., 1996, 315, L357–L360 CAS.
  2. E. L. Gibb, D. C. B. Whittet, A. C. A. Boogert and A. G. G. M. Tielens, Astrophys. J., Suppl. Ser., 2004, 151, 35–73 CrossRef CAS.
  3. K. I. Öberg, A. C. Adwin Boogert, K. M. Pontoppidan, S. Van Den Broek, E. F. Van Dishoeck, S. Bottinelli, G. A. Blake and N. J. Evans Ii, Astrophys. J., 2011, 740, 1–16 CrossRef.
  4. W. Hagen, A. G. G. M. Tielens and J. M. Greenberg, Chem. Phys., 1981, 56, 367–379 CrossRef CAS.
  5. A. G. G. M. Tielens and L. J. Allamandola, Composition, structure, and chemistry of interstellar dust, National Aeronautics and Space Administration Technical Report, 1987 Search PubMed.
  6. N. J. Mason, A. Dawes, P. D. Holtom, R. J. Mukerji, M. P. Davis, B. Sivaraman, R. I. Kaiser, S. V. Hoffmann and D. a. Shaw, Faraday Discuss., 2006, 133, 311–329 RSC ; discussion 347–374, 449–452.
  7. B. Sivaraman, V. Venkataraman, A. Kalyaan, S. Arora and S. Ganesh, Adv. Space Res., 2015, 56, 2428–2431 CrossRef.
  8. G. B. Hansen and T. B. McCord, J. Geophys. Res., 2004, 109, 1–19 CrossRef.
  9. A. C. A. Boogert, P. A. Gerakines and D. C. B. Whittet, Annu. Rev. Astron. Astrophys., 2015, 53, 1–51 CrossRef.
  10. E. Herbst and E. F. van Dishoeck, Annu. Rev. Astron. Astrophys., 2009, 47, 427–480 CrossRef CAS.
  11. B. T. Draine, Annu. Rev. Astron. Astrophys., 2003, 42, 241–289 CrossRef.
  12. J. S. Mathis, Annu. Rev. Astron. Astrophys., 1990, 28, 37–70 CrossRef CAS.
  13. H. M. Cuppen and E. Herbst, Astrophys. J., 2007, 668, 294–309 CrossRef CAS.
  14. M. P. Collings, M. A. Anderson, R. Chen, J. W. Dever, S. Viti, D. A. Williams and M. R. S. McCoustra, Mon. Not. R. Astron. Soc., 2004, 354, 1133–1140 CrossRef.
  15. A. Fresneau, G. Danger, A. Rimola, P. Theule, F. Duvernay and T. Chiavassa, Mon. Not. R. Astron. Soc., 2014, 443, 2991–3000 CrossRef CAS.
  16. P. Ayotte, R. S. Smith, K. P. Stevenson, Z. Dohnálek, G. A. Kimmel and B. D. Kay, J. Geophys. Res., 2001, 106, 33387–33392 CrossRef CAS.
  17. H. Kobayashi, H. Hidaka, T. Lamberts, T. Hama, H. Kawakita, J. Kästner and N. Watanabe, Astrophys. J., 2017, 837, 15 CrossRef.
  18. H. Hidaka, N. Miyauchi, A. Kouchi and N. Watanabe, Chem. Phys. Lett., 2008, 456, 36–40 CrossRef CAS.
  19. R. T. Garrod, Astrophys. J., 2013, 778, 158 CrossRef.
  20. T. Hama and N. Watanabe, Chem. Rev., 2013, 113, 8783–8839 CrossRef CAS PubMed.
  21. Y. Oba, N. Miyauchi, H. Hidaka, T. Chigai, N. Watanabe and A. Kouchi, Astrophys. J., 2009, 701, 464–470 CrossRef CAS.
  22. E. C. Fayolle, J. Balfe, R. Loomis, J. Bergner, D. Graninger, M. Rajappan and K. I. Oberg, Astrophys. J., Lett., 2016, 816, L28 CrossRef.
  23. A. Narten, C. Venkatesh and S. Rice, J. Chem. Phys., 1976, 64, 1106–1121 CrossRef CAS.
  24. I. M. Svishchev, P. G. Kusalik and N. Scotia, J. Phys. Chem., 1993, 99, 3049–3058 CrossRef CAS.
  25. L. Fu, A. Bienenstock, S. Brennan, L. Fu, A. Bienenstock and S. Brennan, J. Chem. Phys., 2009, 131, 1–10 Search PubMed.
  26. U. Raut, M. Famá, B. D. Teolis and R. A. Baragiola, J. Chem. Phys., 2007, 127, 1–6 CrossRef PubMed.
  27. K. Isokoski, J.-B. Bossa, T. Triemstra and H. Linnartz, Phys. Chem. Chem. Phys., 2014, 16, 3456–3465 RSC.
  28. S. Cazaux, J.-B. Bossa, H. Linnartz and A. G. G. M. Tielens, Astron. Astrophys., 2015, 573, 1–7 CrossRef.
  29. J. Keane, A. Boogert, A. Tielens, P. Ehrenfreund and W. Schutte, Astron. Astrophys., 2001, 375, L43–L46 CrossRef CAS.
  30. T. Lamberts, X. D. Vries and H. M. Cuppen, Faraday Discuss., 2014, 168, 327–347 RSC.
  31. T. Lamberts, H. M. Cuppen, S. Ioppolow and H. Linnartz, Phys. Chem. Chem. Phys., 2013, 15, 8287–8302 RSC.
  32. G. A. Kimmel, Z. Dohnálek, K. P. Stevenson, R. S. Smith, B. D. Kay, K. P. Stevenson, R. S. Smith, B. D. Kay, G. A. Kimmel and Z. Dohna, J. Chem. Phys., 2001, 114, 5295–5303 CrossRef CAS.
  33. Z. Dohnálek, G. A. Kimmel, P. Ayotte, R. S. Smith, B. D. Kay, G. A. Kimmel, P. Ayotte, R. S. Smith, B. D. Kay and Z. Dohna, J. Chem. Phys., 2003, 118, 364–372 CrossRef.
  34. G. A. Kimmel, K. P. Stevenson, Z. Dohnálek, R. S. Smith, B. D. Kay, R. S. Smith and B. D. Kay, J. Chem. Phys., 2001, 114, 5284–5294 CrossRef CAS.
  35. Y. G. Yang, X. W. Zhou, R. A. Johnson and H. N. G. Wadley, Acta Mater., 2001, 49, 3321–3332 CrossRef CAS.
  36. X. W. Zhou, R. A. Johnson and H. N. G. Wadley, Acta Mater., 1997, 45, 1513–1524 CrossRef CAS.
  37. T. I. Hasegawa, E. Herbst and C. Ming Leung, Astrophys. J., Suppl. Ser., 1992, 82, 167–195 CrossRef CAS.
  38. D. E. Brown, S. M. George, C. Huang, E. K. L. Wong, K. B. Rider, R. S. Smith and B. D. Kay, J. Phys. Chem., 1996, 100, 4988–4995 CrossRef CAS.
  39. W. Klopper, J. G. C. M. V. D.-v. D. Rijdt, F. B. V. Duijneveldt, T. C. Group, P. O. Box and T. B. Nl, Phys. Chem. Chem. Phys., 2000, 2, 2227–2234 RSC.
  40. E. M. Mas, K. Szalewicz and E. M. Mas, J. Chem. Phys., 1996, 104, 7606–7614 CrossRef CAS.
  41. M. S. Westley, G. a. Baratta and R. a. Baragiola, J. Chem. Phys., 1998, 108, 3321–3326 CrossRef CAS.
  42. B. S. Berland, D. E. Brown, M. A. Tolbert and S. M. George, Geophys. Res. Lett., 1995, 22, 3493–3496 CrossRef.
  43. Z. Dohnálek, G. A. Kimmel, P. Ayotte, R. S. Smith and B. D. Kay, J. Chem. Phys., 2003, 118, 364–372 CrossRef.
  44. A. Fredon, T. Lamberts and H. M. Cuppen, Astrophys. J., 2017, 849, 125 CrossRef.
  45. K. P. Stevenson, G. A. Kimmel, Z. Dohna, R. S. Smith and B. D. Kay, Science, 1999, 283, 1505–1508 CrossRef CAS PubMed.
  46. K. Furuya, M. N. Drozdovskaya, R. Visser, E. F. Van Dishoeck, C. Walsh, D. Harsono, U. Hincelin and V. Taquet, Astron. Astrophys., 2017, 599, 18 CrossRef.
  47. T. Albertsson and D. Semenov, Astrophys. J., 2014, 784, 11 CrossRef.
  48. E. M. Penteado, C. Walsh and H. M. Cuppen, Astrophys. J., 2017, 844, 71 CrossRef.

This journal is © the Owner Societies 2018