THREE-DIMENSIONAL, OFF-LATTICE MONTE CARLO KINETICS SIMULATIONS OF INTERSTELLAR GRAIN CHEMISTRY AND ICE STRUCTURE

The model adopts the numerical MC approach suggested and proved by Gillespie (1976) and outlined below; the general method has been successfully used in the past for astrochemical simulations by several authors (e.g., Charnley 2001; Chang et al. 2005; Cuppen & Herbst 2007; Vasyunin et al. 2009). In common with the work of Cuppen & Herbst, the new model takes a microscopic view of the grain-surface chemistry, but otherwise uses a new methodology for the treatment of the chemistry and the physical positions of atoms and molecules on the grain.

Gillespie (1976) treats the time evolution of a chemical kinetic system as a series of consecutive events. The two key quantities to be determined by the probabilistic approach at any time, t, are (1) which event occurs next and (2) how much time elapses between this and the previous event. To determine this information, each possible chemical process, i, is assigned a rate, R_i (in units of s⁻¹). In the new model, such processes include an accretion event, a desorption event, or an individual thermal hop from one binding site to another; in the case of desorption and diffusion, specific rates are calculated for each individual surface particle. Following Gillespie, each process is then assigned a probability equal to its own rate divided by the sum of all rates, R_tot. A computer-generated random number, 0 ⩽ N_{ran, 1} < 1, then chooses from these probabilities which process will occur next. Another random number, N_{ran, 2}, determines the time elapsed since the last event, equal to Δt = −ln (N_{ran, 2})/R_tot. The populations, positions, rates, or any other relevant information for each affected chemical species is then updated and the above steps are repeated until the user-defined end time of the simulation is reached. In the model, this method is used to determine the overall sequence of thermal hopping, reaction, desorption, or accretion of every surface atom/molecule, according to their individual rates.

The new model explicitly takes into account the positions of particles on the grain surface, without using a pre-determined position lattice. However, it does require the grain-surface structure to be defined initially. The model will allow chemical simulations on any arbitrary grain size and structure.

For this purpose, a simple code has been constructed to produce coordinates for every atom in an approximately spherical dust grain of radius 5 atoms; see Figure 1. The grain may be simplistically considered to be composed of carbon atoms, although the purpose of the generated grain is to provide a simple surface structure upon which the new model may be initially tested. Atoms are positioned in a simple cubic structure and assigned the representative radius of 1.6 Å (which varies only slightly from the canonical van der Waals radius for carbon atoms of 1.70 Å; Bondi 1964). The grain has a diameter of approximately 32 Å. Future applications of the generating code will be used to produce alternative/arbitrary grain sizes and morphologies, as well as more specific surface structures, compositions, and atomic sizes.

The total rate of accretion of gas-phase material onto the grain is dependent on the cross section presented by the grain and its surface ice mantle. In simpler models, this rate is typically approximated by assuming a spherical grain, using a mean radius corresponding to either a bare grain (Hasegawa et al. 1992) or the combined grain and mantle (Acharyya et al. 2011).

However, unlike most past models, the new method explicitly tracks the position of each surface species and allows the non-uniform build-up of surface structures. Thus, the treatment of accretion adopted in this model must account for the trajectory of the incoming particle as well as the rate.

In the new model, accretion from the gas phase is handled in several stages; first, a basic accretion rate is constructed for each gas-phase chemical species, based on the cross section of the smallest sphere that can contain the entire grain and ice mantle. Following Hasegawa et al. (1992),

$$\begin{equation} R_{\mathrm{acc}}(i) = \sigma _{\mathrm{sphere}} \, \langle v(i)\rangle \, n(i) \end{equation} \tag{ 1 }$$

for the sphere bounding a single dust grain, where v(i) and n(i) are the thermal velocity and concentration of the accreting gas-phase species, i, and σ_sphere = π r², where r is the radius of the sphere. (Because the bounding sphere is by definition larger than the combined mantle and dust grain, this basic rate is larger than the effective rate ultimately produced by the model, as explained below). An independent accretion rate for each gas-phase species is constructed in this way.

At such time that the MC routine selects, according to the rate given by Equation (1), the accretion of a particular species to be the next process in the sequence, a randomized entry point into the bounding sphere is then generated in spherical coordinates:

$$\begin{eqnarray} &&\theta = 2 \, \pi \, N_{\mathrm{ran,3}}, \end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray} &&\phi = \cos ^{-1} \, (2 \, N_{\mathrm{ran,4}} -1), \end{eqnarray} \tag{ 3 }$$

where N_{ran, 3} and N_{ran, 4} are random numbers between 0 and 1. Equations (2) and (3) result in a uniform distribution of points over the surface of the bounding sphere.

Expressions of the same form as Equations (2) and (3), using new random numbers, are then used to determine a trajectory for the incoming particle that is randomized over solid angle. However, only inward trajectories are considered (outward trajectories are reversed to point inward), thus avoiding the double counting of cases where a particle enters from an alternative point on the sphere.

The randomly generated trajectory is traced through the interior of the sphere until a grain/ice particle is met by the incoming particle. At this stage, the accretion is judged to be successful and another routine finds the nearest potential well to the contact point, assuming Lennard-Jones potentials for all interactions between the accreting particle and the grain/ice particles (see Section 2.3.2).

If the trajectory of the incoming particle does not intersect with any atoms/molecules on the grain or surrounding ice mantle, then the accretion was unsuccessful. In such a case, the clock is nevertheless advanced and the MC selection process continues as usual. The entry of a particle into the bounding sphere is still a valid process, even if no accretion ultimately results. The choice of radius of the bounding sphere could be made arbitrarily large (producing a very low accretion efficiency for particles entering the bounding sphere), with no loss of accuracy to the model results. The purpose of the adoption of the smallest bounding sphere possible is simply to avoid unnecessary and time-consuming calculations.

In the results presented, only atomic H and O are allowed to accrete from the gas phase, with an assumed sticking coefficient of unity. Fractional abundances of 2 × 10⁻⁴ are adopted for each species. The overall gas density is varied in the models, but the fractional abundances are assumed to remain constant over time and no gas-phase chemistry is treated in the current model. Each of these restrictions may be lifted in future work.

MC treatments of gas-phase chemistry need consider only the evolution of, and rates associated with, the total population of each gas-phase chemical species. However, a microscopically exact treatment of ice structure requires the positions and behavior of individual grain-surface atoms and molecules to be simulated explicitly.

In this model, rates for diffusion and desorption are calculated for each individual grain-surface atom or molecule. Diffusion and desorption are treated as thermal processes, with individual rates of the form:

$$\begin{equation} R(j) = \nu _{j} \, \exp (-E_{j}/T), \end{equation} \tag{ 4 }$$

where j is a specific atom or molecule, ν_j is the characteristic vibrational frequency of the species in question, T is the dust temperature, and E_j is either the binding energy or diffusion barrier, depending on which process is being considered. In this paper, the terms "thermal desorption" and "evaporation" are considered synonymous.

In practice, the MC procedure initially selects a particle to undergo a thermal process, using the rate of the fastest thermally activated process allowed for each surface particle in order to make the selection (although these rates also compete with the accretion rates described above). From here, a random selection is made between desorption and each and every diffusion process that the chosen particle may undergo, weighted according to each individual rate. This two-stage approach is necessary because, for example, the probability that a particle attains sufficient energy to diffuse also contains the small probability that it attains sufficient energy to desorb entirely; Equation (4) represents the rate at which a particle reaches an energy of E_j or higher. These processes are therefore in competition with one another and must be treated as such. In the same way, for the case of a particle that has, for example, four available diffusion paths with equal energy barriers, the diffusion rate is not increased four-fold, but split four ways.

The atoms that constitute the dust grain surface itself are not allowed to diffuse or desorb, but do contribute to the binding energies of surface-bound species.

The treatments applied to desorption and diffusion processes are described in more detail below.

2.4.1. Thermal Desorption

When an atom has just accreted onto the grain surface, it initially resides in a potential well whose position and binding energy are calculated at the time of accretion. (It should be noted that the position assigned to each particle represents only its mean position; in reality, all particles would be vibrating in their potential wells). In the model, a particle is considered to be bound to a surface if it is bonded to three or more other atoms/molecules—these could be atoms of the dust grain itself or other particles that are themselves bound to the grain surface. Two particles, A and B, are considered to be bonded to each other if their center-to-center separation, s_AB, is approximately equal to σ = 3.2 Å (specifically, in the range 0.9σ < s_AB < 1.1σ; closer bonding is explicitly prohibited in the model, in consideration of van der Waals repulsion).

For each individual atom/molecule, the binding energy (i.e., the energy required to desorb/evaporate into the gas phase), E_des, is equal to the sum of the interactions, , with the other particles to which it is bound (note that Equation (5), below, is not used for these calculations). This value is used to evaluate the evaporation rate given by Equation (4) and in the calculation of ν_j, following Hasegawa et al. (1992). Any particles that are spatially "boxed-in," under the hard-sphere approximation, are not allowed to evaporate or diffuse. Thus, molecules within bulk ice mantles are fixed in place unless the outer-surface molecules are removed.

Interaction potentials are defined for the pairing of every type of chemical species in the model, as shown in Table 1. Potentials between any chemical species and a particle in the dust grain itself are chosen to reproduce approximately the binding energies used in previous, rate-based gas-grain chemical models, under the assumption that binding to the bare grain involves bonding with four dust grain atoms. For example, the binding energy of atomic hydrogen assumed by Garrod & Herbst (2006) of 450 K may be compared with a value in this model of 400 K, which would be obtained for bonding to a flat face of the dust grain used in the model (see Figures 1 and 2). Bonding to only three grain atoms, or to as many as five, is possible in some positions on the bare grain surface, which would produce a binding energy of 300 K or 500 K, respectively. For atomic oxygen, the 200 K interaction potential with a dust grain atom adopted in this model reproduces the commonly used binding energy of 800 K (Tielens & Allamandola 1987), for most positions on the grain.

Zoom In Zoom Out Reset image size

Table 1. Pair-wise Interaction Potentials, in Units of K

	Grain	H	H2	O	O2	OH	H2O	H2O2
H	100	100	50	100	100	100	100	100
H2	50	50	50	50	50	50	50	50
O	200	100	50	200	200	200	200	200
O2	300	100	50	200	300	300	300	300
OH	400	100	50	200	300	400	500	600
H2O	500	100	50	200	300	500	1000	1000
H2O2	600	100	50	200	300	600	1000	1200

Download table as:  ASCIITypeset image

Interaction potentials between particles that are not a part of the grain itself are nevertheless guided by the grain-binding values. In general, the interaction potential of a pairing A–B is assigned the lower of the values for A–gr and B–gr (where "gr" represents a dust grain atom). An exception is permitted for species expected to demonstrate strong hydrogen-bonding properties. Thus, the interaction potentials involving H₂O, H₂O₂, and OH (when bound to one of those two molecules) are augmented by several hundred Kelvin. For H₂O–H₂O and H₂O₂–H₂O₂ bonding, the A–gr values are doubled. This produces a four-partner binding energy for water of 4000 K, although in practice the model may yield water ice that is sufficiently amorphous to produce many surface water molecules that are bound to five or six other waters. The generic value adopted by Garrod & Herbst (2006) was 5700 K.

The binding treatment assumed throughout the model may be characterized as a "hard-sphere, nearest-neighbor" approach. The interaction potentials assumed in the model are isotropic; polar molecules may in reality exhibit more directionality, but such considerations are currently beyond the scope of this model.

2.4.2. Diffusion

Because the model uses an "off-lattice" approach, the positions of the adjacent potential wells into which a particle could diffuse are not initially specified and must be calculated—as needed—from the current physical state of the system, as defined by the positions of all nearby particles. The method set out for this model allows the energy barrier for each possible diffusion pathway of a given surface particle to be calculated at the same time as the total binding energy, prior to any calculation of a final post-diffusion position associated with it. The putative final position resulting from a thermal hop need therefore only be calculated if that particular hop has actually been selected as the next process.

To determine diffusion barriers, the model requires information about the possible pathways that a diffusing particle may take from the current binding site to another. For the illustrative case of an H atom bound to a flat face of the grain (see Figure 2), the atom is bound to four surface-binding partners. The atom may also be seen to have four possible diffusion pathways, corresponding to diffusion along the surface between pairs of binding partners. While in theory there are an infinite number of paths that a diffusing atom may take, these four constitute those that pass through a saddle point in the surface potential: they are, by far, the most probable paths. For each saddle-point pair, a diffusion barrier is defined to be equal to the total binding energy, minus the interaction potentials of the diffusing particle with each of the pair of atoms through which the saddle point would pass. In other words, diffusion requires the breaking of two of the four bonds, while the other two remain intact, even after the hop is finished (Figure 2). In a different case where, for example, the atom were bound to only three binding partners, only one bond would need to be broken. This approach is consistent with typical values of the ratio of diffusion barriers to binding energies, estimates of which can range from ∼0.3 to 0.8 in current models (see, e.g., Garrod & Pauly 2011). Such values may be achieved through surface binding to different numbers of binding partners; the minimum of three partners (of equal interaction potential) would produce a value E_diff: E_des = 1/3, while values closer to unity may be attained for binding to highly irregular (rough) surfaces.

The calculated diffusion barriers are stored and used to determine the selection of a particle and a thermal process as described in Sections 2 and 2.3. If diffusion is selected, the next step is to determine the final position of the particle after diffusion through the selected saddle point.

In order to determine this position, the diffusion process is treated as a rotation around the saddle-point pair (with which the bonds remain intact following diffusion). The particle is rotated at a fixed radius around the axis joining the saddle-point pair, until it is close enough to be considered "bonded" to the atoms on the other side of the saddle point.

At this stage, the position of the particle is optimized iteratively to find the deepest part of the combined interaction potentials of the particles to which it is now bound. Each of these interaction potentials is treated as a Lennard-Jones 6–12 potential:

$$\begin{equation} V_{{\rm LJ}}=\, \epsilon \, \left[\left(\frac{\sigma }{s}\right)^{12} - 2 \, \left(\frac{\sigma }{s}\right)^{6}\right], \end{equation} \tag{ 5 }$$

where s is the current center-to-center separation, is the depth of the potential, and σ is the separation at which the potential minimum is reached, which is set to 3.2 Å for the present model. (This expression is used wherever distance-dependent, pair-wise potentials are required in the model; for example, the binding positions of newly accreted particles are also calculated iteratively using this expression.) The derivatives of the potentials are used to calculate the net force vector resulting from the combined field at the current position. The position of the particle is advanced incrementally according to the force vector, to yield the position of the potential minimum to a tolerance of 0.005 Å.

After the position is optimized in this way, the new binding energy is, again, taken simply as the sum of the strengths of the interaction potentials, , with all the species to which the particle is bound (i.e., 0.9σ < s < 1.1σ); new diffusion barriers for all viable saddle-point pairs are also calculated at this time.

The viability of saddle-point pairs is defined in the model as those pairs for which a rotation would be unhindered by other particles to which the diffusing species is already bound. For example, in the case shown in Figure 2, while there are four pairs that would allow rotation (and thus diffusion), a pairing of diagonally arranged particles from this group (Figure 2(c)) is disallowed, because rotation/diffusion around the axis between the diagonal particles would be obstructed in either direction by another binding partner.

All rotation pairings of particles to which the species in question is bound (i.e., N(N − 1)/2 pairings, for N bonds) are initially allowed. Each pairing is then assessed according to the following test: a plane is defined by the initial positions of the diffusing particle and each of the particles of the rotation pairing. If the other particles to which the surface species is bound are either all above or all below the plane, then diffusion is allowed (in a direction away from those particles). Otherwise, diffusion is disallowed, as a collision would result (as for the diagonal pairing shown in Figure 2(c)). Typically, only pairings between adjacent (although not necessarily bonded) particles are allowed by this test, although irregular surface morphologies may allow other cases to be permitted.

In this way, the possible diffusion paths on an arbitrary surface may be determined. Species for which all rotation pairs are forbidden are considered "boxed in" and are thence also prohibited from desorption, until such time as this condition changes due to diffusion or evaporation of other particles.

It should be noted that, following accretion, desorption, or diffusion, the diffusion and desorption properties of all affected species—such as binding partners both prior to and following a hop—are re-calculated and the rates and diffusion viabilities are re-assessed.

The rotation itself should not be considered an accurate description of the motion of the particle, but rather a convenient method for finding the potential minimum to which the particle hops, via a potential saddle point. Actual diffusion paths for arbitrary binding structures would require intensive molecular dynamics or quantum chemical calculations.

In the case given in the example (Figure 2), the diffusion barrier associated with each of the four possible hops is the same. However, if the particle were bound to different chemical species, differences in the magnitudes of the interaction potentials would lead to differences in the probabilities of diffusion along each path. On an irregular surface, the particle could be bound to three or more other particles; the number of binding partners as well as their chemical species is therefore important to the energy barrier to diffusion (i.e., the energy required to break all but two of the bonds).

The model is also optimized such that previous hops of the currently active particle are stored in memory, so long as no other particles diffuse, desorb, or accrete in the meantime. This decreases run times by several orders of magnitude; the calculations are often dominated by the diffusion of a single weakly bound atom (until it finds a stronger binding site).

2.4.3. Reaction

Many models include a number of other grain-surface processes not included here, such as photo desorption and photo dissociation. The latter is often suggested to be an important mechanism for the formation of complex organic molecules in the interstellar medium (e.g., Garrod & Herbst 2006). All such processes are expected to be included in future versions of the model, but are omitted here so that the basic approach may be tested.

The inclusion in the model of photo induced desorption based on estimated or measured rates (Öberg et al. 2009) would be technically trivial. Photo dissociation would be somewhat more complex, as it would require the positioning of two surface products that originated from a single molecule.

Consideration of the formation of complex organic molecules using the new model is anticipated, but will depend on the adoption of a significantly larger reaction network that includes CO and its related products. Furthermore, only surface processes are considered in the current model; photo processing of the sub-surface ice mantle material would involve a considerable increase in technical complexity.

In order to investigate the formation of structure in the ices, simulations have been run to produce ice mantles several times thicker than the diameter of the underlying grain. Four sets of models have been produced, corresponding to gas densities in the range 2 × 10⁴–2 × 10⁷ cm⁻³. For each density, the model has been run three times, adopting a different initial random number seed. This is done to demonstrate that the behavior of the models is a general feature and not dependent on the specific random numbers used in each run. All the simulations are run until 200,000 water molecules have been formed on the grain. As before, dust and gas temperatures are held at 10 K and the gas-phase fractional abundances of H and O atoms are fixed at 2 × 10⁻⁴.

In addition, a single model has been run in which water molecules are accreted directly from the gas phase, with no accretion of H and O atoms. An arbitrary gas density of 2 × 10⁷ cm⁻³ is used; however, for the consideration of ice structure, the precise rate of accretion is unimportant, due to the lack of surface chemistry and the negligible rate of water diffusion caused by the low temperature.

3.1.1. Structure and Porosity

Figure 4 shows the final state of each density model (run 1 is shown in each case), including pure accretion (panel a). Hydrogen molecules (H₂) are highlighted in blue to aid the eye.

A number of features are immediately apparent; first, the large- and small-scale structure of the ice is strongly dependent on the gas density. The most irregular structure results from the higher-density simulations, culminating in the accretion-only model—in this case, water molecules remain wherever they land as the result of accretion. In the simulations with active chemistry, higher densities result in more rapid hydrogenation of surface oxygen atoms, before they have the opportunity to find alternative, stronger binding sites.

The higher-density ice mantles exhibit a highly porous or "creviced" structure. Panels (b) and (c) in particular show a "cauliflower-like" structure, wherein the crevices act to isolate large-scale nodules of ice. In the lower density simulations, the nodules become broader and the crevices become less deep. The lowest density model indeed shows a fairly smooth ice structure with no deep crevices or voids.

In the accretion-only simulation, the ice appears irregular on all scales. The small-scale porous structure is largely the result of the lack of surface diffusion, such that new species are unable to move even a short distance into the stronger binding sites that are provided by areas of extreme small-scale curvature of the ice surface. However, the large-scale irregularity of the ice is caused by the combination of fully randomized trajectories with a three-dimensional grain; the growth of randomly produced surface irregularities, or "bulges," is amplified by their protrusion into the accretion field, picking up more material from a greater solid angle, while blocking out those trajectories for other regions of the grain/ice surface. Such effects are still prevalent in the active-chemistry models—panel (c) shows an ice mantle whose underlying dust grain is centered in the image; a large protrusion may be seen on the right. While a moderate degree of surface diffusion may act to close over the small-scale pores or crevices in the ice, oxygen atoms are required to diffuse a significant distance over the ice surface to counteract the large-scale irregularities. Consequently, while smooth on the small scale, the lowest density simulation (panel e) nevertheless shows an overall irregular shape.

Figure 5 shows cross sections of the same final ice mantles. Each cross section passes through the dust grain center at an arbitrary angle, showing molecules that sit within a "slice" 3σ deep—this allows space for up to two molecules along a line of sight. Animations of panels (a), (b), and (e) are available in the online version of the journal, showing cross sections over 360 deg.

The detail of the porosity within the ice mantles is clearly visible in Figure 5. Panel (a) shows significant and fairly uniform porosity, with pores or crevices that pass deep within the mantle almost down to the dust grain itself. The highest density simulation with active chemistry is not as extreme, but also shows very deep pores/crevices in places and exhibits some closed pores. The degree of porosity falls away at lower gas densities, until, for n_H = 2 × 10⁴ cm⁻³, there are essentially no open pores and the ice is extremely compact.

One of the most striking features of Figure 5 is the arrangement of hydrogen molecules (H₂). Shown in blue in Figures 4 and 5 for ease of identification, H₂ is almost completely segregated from the water molecules, especially in the higher-density simulations. Furthermore, the H₂ molecules are arranged in "veins" that are similar in appearance to the empty porous structures. As may be seen three-dimensionally in Figure 4, the H₂ indeed fills in crevices between larger structures. The veins of H₂ are typically no more than around five molecules across (commensurate microporous structure), while many of the unfilled pores are significantly larger.

Following its formation on the grains, an H₂ molecule may diffuse around the surface. If it finds a weak binding site, it may evaporate entirely, but if it finds a strong binding site—such as may be found in a pore, where multiple binding potentials converge—then it may remain bound in place. The addition of further hydrogen molecules could then contain it and allow a pore to fill up with hydrogen. The presence of larger, unfilled pores in the higher-density simulations suggests that these pores have an insufficient degree of inward curvature to produce binding potentials large enough to retain H₂ for long enough before the deposition of new material can lock it into place.

Figure 6 shows a bar chart of the number of surface bonds of each molecule (of any species) within the ice mantle, for each of the simulations shown in Figures 4 and 5. The accretion-only model shows a very different distribution from the active-chemistry models. Approximately 35% of all molecules in this simulation experience only 3–5 bonds, making them (typically) surface molecules. Conversely, only around 6% of molecules in even the highest density active-chemistry simulation are so weakly bound and the number of species with only three bonds is negligible. However, this does not imply that potential minima that would allow only three bonds do not exist—only that they are unoccupied. Pores, and the associated H₂ veins, may begin to form in such positions, with ice structure building up around them.

Zoom In Zoom Out Reset image size

The porosity visible in Figures 4 and 5 is reflected in the bonding distributions of Figure 6, with the most compact ice (produced by the lowest density model) showing the greatest bias toward 10–12 bonds, which are necessarily bulk ice molecules. Such distributions may therefore be useful as a measure of porosity in these simulations.

Figure 7 shows the fraction of bonds shared with an H₂ molecule, averaged over all H₂ molecules. The highest density model shows the greatest degree of clustering of H₂ molecules. Based on a purely statistical comparison of abundances, the expected value would be around 13%–14% for this density (based on values in Table 3). For the lowest density model, the statistical expectation would be ∼8%. In each case, the actual fractional bonding of H₂ to H₂ is close to 3 times higher than the statistical expectation. It should also be borne in mind that the binding potentials used for H₂ (Table 2) are the same for all binding partners. The H₂ clustering is thus the result of the structural arrangement of the ice mantle to produce pores or concentrations of strong binding sites in which H₂ molecules may become trapped, rather than a preference for H₂ molecules to bind to each other, per se.

Zoom In Zoom Out Reset image size

Table 2. Surface Reactions

#	Reaction
1	H + H → H2
2	O + O → O2
3	H + O → OH
4	H + OH → H2O
5	OH + OH → H2O2

Note. All reactions are presumed to proceed without an activation barrier.

Download table as:  ASCIITypeset image

3.1.2. Chemistry

Figure 8 shows time-dependent, grain-surface chemical abundances for each of the stable molecules included in the model, for all four gas density values. Results for all three runs at each density are plotted in the same panel. The end time of each model corresponds to the time at which 200,000 water molecules have been formed; the precise value varies according to the gas density and, to some degree, according to random variation between runs at the same density. For O₂ and H₂O₂, which have generally low abundances, there is some random fluctuation apparent at early times and there are small systematic differences between runs of the same density. However, as a proportion of the total quantities, these variations fall away as the ice mantles grow.

Zoom In Zoom Out Reset image size

Table 3 shows the final values for each model run, as well as the total amount of H₂ formed, much of which evaporates back into the gas phase. H₂O₂ shows the largest variation between same-density runs, of around 7%–8%, but variations for the other species are typically 1%–2%. The elevated value for H₂O₂ may be due to its more intricate formation mechanism in this model (see Section 3). For reference, Table 3 also indicates the end times of the models, t_f, and the greatest radial distance of any ice-mantle particle from the centroid of the grain, r_max, achieved by the end of each run.

Table 3. Final Time, Maximum Radius, and Final Number of Grain-surface Particles, by Species, Bound to the Grain (and, in the Case of H₂, the Total Formed on the Grain), for Each Model Run

nH	Run	tf	rmax	H	O	OH	O2	H2O2	H2	H2
(cm−3)		(yr)	(Å)						(on grain)	(total)
2 × 104	1	7.047 × 104	147.1	3	1	3	1,238	240	18,284	207,313
	2	6.916 × 104	152.6	6	3	3	1,241	261	18,399	206,096
	3	6.994 × 104	149.1	1	4	0	1,260	257	18,407	208,248
2 × 105	1	8.471 × 103	153.9	6	7	8	1,599	304	23,011	209,973
	2	8.214 × 103	148.3	9	4	5	1,547	301	22,998	206,612
	3	8.175 × 103	160.0	7	6	8	1,523	325	23,494	208,044
2 × 106	1	8.216 × 102	164.9	37	30	19	3,073	620	28,872	217,227
	2	8.292 × 102	176.0	27	23	19	2,974	618	29,216	217,463
	3	8.298 × 102	178.3	27	20	18	3,031	596	29,620	216,183
2 × 107	1	6.687 × 101	191.9	72	42	52	5,308	1,321	32,824	225,258
	2	6.695 × 101	177.3	60	53	45	5,385	1,233	31,484	228,142
	3	6.978 × 101	183.6	81	42	42	5,381	1,281	33,319	227,244

Note. All models were stopped after the formation of 200,000 water molecules.

Download table as:  ASCIITypeset image

Abundances of species other than H₂O are seen to be somewhat lower for the lower density models, including H₂. However, the total amount of H₂ formed (see Table 3) is much less affected than the quantity that is stored on the grain at each gas density. Accreted atomic hydrogen is typically able to remain on the grain for long enough for another reactive species to accrete and react with it, for all density models. In the case where the product is H₂, it is more likely to evaporate for the lower density models, which have a more regular surface that provides fewer strong-binding positions or pores. In each of the simulations, the amount of H₂ formed is only slightly larger than the total amount of water formed.

The running of the new model over various density conditions allows an assessment of the computational efficiency to be made. Each simulation was run with an Intel Xeon X5482 cpu, using a single thread.

The run time for all of the three higher density simulations (2 × 10⁵, 2 × 10⁶, and 2 × 10⁷ cm⁻³) is approximately 90–95 hr, depending on which random number seed run is considered. Only the lowest density simulations show an apparent density dependence, taking ∼180 hr to run.

The majority of cpu time is used up on the calculation of diffusion pathways and on the assignment of position and pairing information associated with each surface hop. Because the model keeps diffusion paths in memory, the most intensive calculations are typically done only once per particle, providing the positions of surface potential minima for an individual, mobile particle for large numbers of subsequent, identical hops, prior to reaction, desorption, or trapping in a strong binding site.

It appears that the larger drain on cpu time in the low-density case is caused by the relative smoothness of the ice surface. In this case, a diffusing particle (typically atomic hydrogen) on the grain has a smaller probability of becoming trapped in a strong binding site (thus ceasing its diffusion), while the probability of another particle being accreted onto the grain and subsequently reacting with the first particle is smaller (due to lower gas density). This results in a longer sequence of diffusion before a reaction occurs.

It is to be expected that the use of larger grains would ameliorate this effect, as the current model dust grain is extremely small, providing low overall rates of accretion. Conversely, higher temperatures would most likely make it worse, due to the resulting faster diffusion, although desorption rates would also increase.

The inclusion of a broader network of surface species is expected to make the surface more heterogeneous, such that the smoothness of the surface potential would be reduced, minimizing the possibility that surface species may spend long periods of cpu time without being trapped in a strong binding site.

This model allows porous structures within the ice mantle to form naturally, as the result of chemical and physical conditions. However, the consideration of porosity in simulations of interstellar dust grain chemistry is not new; Perets & Biham (2006) used a rate-based model to study the effects of porosity within the dust grain structure itself. The models indicated that porosity could increase the efficiency of H₂ formation over a wider range of temperatures than would otherwise be possible. Taquet et al. (2012) extended the treatment of Perets & Biham (2006) to a full gas-grain chemical model that included the formation of a porous ice mantle, using pre-defined values for the number and sizes of pores. They found only moderate increases in the production of certain key ice-mantle species.

Cuppen & Herbst (2007) used an MC technique to treat the formation of ice on a surface with periodic boundary conditions, using a fixed lattice for the positions of particles on the grain or within the ice mantle. Depending on the physical conditions, ice structures of varying compactness and porosity were produced, including some with tower-like structures of apparently single-molecule thickness. The models that show the most extreme porosity, however, were obtained assuming temperatures rather higher than the 10 K used here.

In the present model, only a limited reaction set is considered, which does not include any activation barrier-mediated reactions. Nevertheless, the basic reactions that allow the formation of the most important species, such as H₂ and H₂O, are indeed present. The use of an off-lattice approach in this model requires that the binding of atoms or molecules to the surface involve at least three other binding partners. Thus, the "sky-scraper" structures produced by the models of Cuppen & Herbst are not found here, both because binding to a single partner cannot occur and because the direction of the binding is not pre-determined, but depends on the arrangement of the binding partners.

Perhaps most importantly, the off-lattice approach of the present model allows the effects of the curvature of the surface to be automatically taken into account. Inward curvature, such as may be found within the pores, produces stronger binding conditions resulting from the ability to bind with a greater number of partners; likewise, the outward curvature found on the bare grain and on the "nodules" or "bulges" that form on the ice, tends to minimize the number of binding partners. These curvature effects drive mobile particles into the pores as the pores are forming, so that they become filled with H₂, producing veins. Larger pores—which are formed in the higher density models and which experience less extreme curvature—tend not to fill up in this way.

The kind of micropores investigated by Taquet et al. (2012), which appear to be comparable in size to the H₂ veins that form here, may not therefore be long-lived enough to affect the chemistry significantly. The present model indicates that such pores would be filled up, although others may develop in their place. In general, however, the low-density model (which perhaps provides the most appropriate comparison) tends to predict a compact, albeit segregated, ice structure.

The gas-phase chemical abundances used in the model presented here represent dense-cloud conditions, under which most hydrogen takes the form of H₂, resulting in similar H and O gas-phase abundances. However, higher values for H are entirely plausible, which would result in greater H fluxes onto the grains and therefore a shorter time period before surface oxygen atoms could find alternative binding sites. It is therefore likely that a cloud with large quantities of H not yet converted to H₂ should produce more porous dust grain ices. It is also probable that, other model quantities being equal, the smallest of these pores would rapidly fill up.

The consideration here of a model that involves the direct accretion of water onto the grain allows a basic comparison with laboratory ASW ices, which are also formed by the deposition of water molecules directly onto a cold surface. It is well known that the angle of deposition has a strong influence on the degree of porosity in the ice (e.g., Kimmel et al. 2001; Raut et al. 2007). The present model does not investigate the effect of deposition angle directly, but uses a completely randomized field of accreting particles. The result is an extremely porous ice, which is qualitatively different from the forms produced by active grain-surface chemistry, because of the minimal diffusion of the water molecules following accretion. In the active-chemistry models, the slow formation of the ice from its constituent atoms allows voids and strong binding sites to be filled more effectively, so that the lowest density simulations show an absence of pores and a relatively small quantity of H₂ veins as compared with higher gas-density (i.e., faster-accretion) models.

It therefore seems appropriate to suggest that interstellar ice analogs formed in the laboratory may be significantly more porous than actual interstellar ices. This would imply that effects observed in laboratory ices, relating to the absorption of other molecules into pores and their subsequent release at higher temperatures (e.g., Collings et al. 2004), may be of less importance in interstellar ice mantles. However, the present model does not include any mechanism for super-thermal surface diffusion of the accreted particle or indeed of newly formed surface molecules. The acceleration of the accreting particle into a surface potential well, or the release of chemical energy from the surface formation of a molecule, could provide sufficient energy for diffusion. This could allow an otherwise immobile particle to find a nearby binding site with stronger binding properties. However, in view of the experimental evidence for significant porosities in laboratory ASW ice (e.g., Westley et al. 1998 and references therein), the importance of this effect may be small.

As shown in Table 3, for decreasing density the H₂:H₂O production ratio approaches 1:1. This may be explained through a simple analysis: taking the 1:1 gas-phase abundance ratio for H:O used in the present model, the ratio of accretion rates is 4:1, due to the 16 times higher mass of oxygen, as entered into Equation (1). Under conditions where atomic hydrogen mobility dominates all surface processes and accretion dominates evaporation, this rate should indeed produce a 1:1 formation rate for H₂:H₂O. (For every 1 oxygen that accretes, 4 H atoms will accrete; two are used up, as they rapidly find the O/OH on the surface, to form water. The other two will, on average, meet each other and form H₂.) This analysis ignores O₂ and H₂O₂ formation, each of which would otherwise remove four and two more hydrogen atoms, respectively, if their oxygen atoms had instead formed surface water. This difference accounts for around 40% of the H₂ excess over H₂O in the low-density models.

It may be seen, therefore, that for the lowest density simulations in particular, the production of H₂ is very close to what may be called a "smooth-grain" limit, in which hydrogen atom mobility is not significantly affected by the presence of strong binding sites.

With increasing density, there are more strong binding sites (due to the more irregular ice structure), while the waiting period for the accretion of a new particle is also shorter. These effects increase the probability that an accreted hydrogen atom will wait for another H to react with it, rather than diffuse to find a free oxygen atom elsewhere on the dust/ice surface, thus producing a greater proportion of molecular hydrogen.

For gas densities less than ∼2 × 10⁴ cm⁻³, one should expect that the ice surface would be similarly smooth and perhaps yet more spherical, while the waiting period for accretion would be greater, all of which should produce H₂:H₂O ratios even closer to unity. Under such circumstances, the use of rate-based models that employ generalized binding and diffusion characteristics should not produce strongly divergent results for H₂ production from those obtained here, provided that appropriate characteristic values are chosen. Such an approach would also require that the stochastic behavior of the chemistry be treated adequately (e.g., the use of the modified rate equations of Garrod 2008 and Garrod et al. 2009), as the scenario described clearly falls into the case where standard rate equations should fail. It is unclear, however, whether such an approach could accurately reproduce the quantity of the resultant H₂ that is retained on the grain surface, which is still determined by heterogeneous structural considerations.

The choice of physical conditions investigated in the present study has been limited to the gas density, which controls the rate of accretion onto the grain. The gas density demonstrates a clear effect on the resultant ice structure and therefore on the retention of H₂ molecules on the grain. The relative formation rates of O₂ and H₂O₂ are also seen to be affected by the density variation, due to the competition between chemical reactions whose rates are dependent on surface diffusion.

Each of the effects discussed above concerns the interplay between accretion, diffusion, and desorption. The gas density is thus not the only parameter that may affect the chemistry and structure of the ice formed on the grain. Future work will address the importance of dust temperature; however, one may predict that temperatures higher than the 10 K used here would produce more diffusion of oxygen atoms, meaning more compact ices. Likewise, lower temperatures would result in more porosity and a greater retention of H₂ within those structures.

The present, preliminary model uses a very small grain, of radius 16 Å, as compared with a canonical value of 0.1 μm (=1000 Å). However, the final ice-mantle radii are significantly larger, at around 150–190 Å. Future models will investigate chemistry on larger grains. The application of the model to a porous dust grain, of explicitly defined structure, would also be quite possible.