Drew A. Christianson
Robin T. Garrod- 1Department of Chemistry, University of Virginia, Charlottesville, VA, United States
- 2Department of Astronomy, University of Virginia, Charlottesville, VA, United States
The degree of porosity in interstellar dust-grain material is poorly defined, although recent work has suggested that the grains could be highly porous. Aside from influencing the optical properties of the dust, porosity has the potential to affect the chemistry occurring on dust-grain surfaces, via increased surface area, enhanced local binding energies, and the possibility of trapping of molecules within the pores as ice mantles build up on the grains. Through computational kinetics simulations, we investigate how interstellar grain-surface chemistry and ice composition are affected by the porosity of the underlying dust-grain material. Using a simple routine, idealized three-dimensional dust-grains are constructed, atom by atom, with varying degrees of porosity. Diffusive chemistry is then simulated on these surfaces using the off-lattice microscopic Monte Carlo chemical kinetics model, MIMICK, assuming physical conditions appropriate to dark interstellar clouds. On the porous grain surface, the build-up of ice mantles, mostly composed of water, leads to the covering over of the pores, leaving empty pockets. Once the pores are completely covered, the chemical and structural behavior is similar to non-porous grains of the same size. The most prominent chemical effect of the presence of grain porosity is the trapping of molecular hydrogen, formed on the grain surfaces, within the ices and voids inside the grain pores. Trapping of H2 in this way may indicate that other volatiles, such as inert gases not included in these models, could be trapped within dust-grain porous structures when ices begin to form.
1 Introduction
Dust grains are known to play a major role in the chemical evolution of interstellar clouds and star-forming regions. Of critical importance is the production of molecular hydrogen on grain surfaces, through the reactions of hydrogen atoms adsorbed from the gas phase (e.g., Gould and Salpeter, 1963; Hollenbach and Salpeter, 1971); much of the H2 produced in this way desorbs again, to become the dominant component of the gas in dense, dark clouds. In clouds of sufficient visual extinction, the dust grains may build up ice mantles composed mostly of water, but with substantial amounts of CO, CO2, NH3 and several other simple species (Boogert et al., 2015). Much of the chemistry that produces these grain-surface molecules involves the surface diffusion of atomic hydrogen, which is understood to be mobile even under the cold (∼10 K) conditions of dark clouds (Senevirathne et al., 2017). By repetitive thermal hopping between surface potential minima (binding sites), H atoms may meet other reactive species, such as atomic O or N, and react. The adsorption of atoms, as well as molecules such as CO, from the gas phase and their subsequent hydrogenation leads to the build up of ice mantles many monolayers in thickness (Cuppen et al., 2017).
Interstellar dust grains are typically considered to be composed of either carbonaceous or silicaceous material, both of which populations may be identified in the interstellar medium. By fitting observational extinction values to models that consider the optical properties of dust grains, Mathis et al. (1977) provided a grain-size distribution for the Milky Way of
In chemical kinetic models of interstellar grain-surface chemistry, in which the grain-surface populations of various atomic and molecular species are traced over astronomical timescales, a canonical, spherical grain radius of 0.1 μm is typically chosen to represent all grains in the model. The grain surface is assumed to comprise around 1 million binding sites, based on a site density of ∼1015 cm−2 (Hasegawa et al., 1992). Some chemical models of interstellar clouds have considered a distribution of grain sizes, each having their own independent chemistry (Acharyya et al., 2011; Pauly and Garrod, 2016, Pauly and Garrod, 2018). However, without the existence of a strong temperature gradient across the size distribution, the influence of the grain size on the surface chemistry in those models appears to be fairly limited.
Most such chemical kinetic models use the so-called “rate-equation” approach (see Garrod, 2008), in which the grain-surface population of each chemical species is traced over time; each chemical species is assigned a fixed binding energy and barrier against diffusion. While efficient, this technique is unable to treat the microscopic structure of either the ice mantle or the underlying dust grain itself, as the positions of individual atoms or molecules on the grain are not considered. The dust grains in these models are thus implicitly smooth and spherical, as are the molecular ices that form on their surfaces. Most models therefore say nothing about the possible influence of grain morphology or porosity on interstellar molecular abundances, either on the grain surfaces or in the gas phase.
Experimental investigations of the direct deposition of interstellar ice analogs suggest the possibility for grain-surface ices to achieve substantial degrees of porosity (e.g., Raut et al., 2007), but no observational evidence has yet been found to indicate that true interstellar ices, which are formed primarily through surface chemistry rather than direct deposition, should be substantially porous; the dangling bond IR feature associated with porous ASW has not so far been detected in the ISM (Keane et al., 2001). Laboratory work indicates that the formation of ices through chemical reactions, rather than direct deposition from the gas phase, appears to be crucial to the degree of porosity. In experiments by Oba et al. (2009), hydrogen atoms and O2 molecules were co-deposited onto Al and D2O surfaces, leading to the build up of water ice films through chemical reactions. Those authors found no dangling-bond features that would indicate porosity, and they suggested that the reaction energy could be responsible for the compactness of the ice. Similar work by Accolla et al. (2013) involving D and O2 deposition likewise found no spectral identifiers of porosity.
Astrochemical models of interstellar ice formation by Garrod (2013a) used a microscopic Monte Carlo technique to trace the formation of water and several other chemical species through atomic addition reactions. Those authors found that the ability of oxygen atoms to diffuse on the grain/ice surface over long periods of time allowed them to find the strongest binding sites, prior to reaction with atomic H, so that the chemically-produced water molecules would automatically be formed in strong binding sites that would otherwise correspond to nascent surface inhomogeneities, thus allowing the ice to form smoothly at microscopic scales. Whether the production of compact water ice is due to the thermal diffusion of atomic oxygen, or to non-thermal, reaction-induced diffusion of the water molecule itself, the chemical production of water on dust grains appears to lead to little or no porous structure, in agreement with the existing observational evidence (although it has also been suggested that the filling of water-ice pores with molecules such as CO could remove the dangling-bond feature; e.g., He et al., 2019).
Similarly, until recently there has been little direct evidence to indicate either the degree or the precise form of porosity, if any, that interstellar dust-grains themselves might assume; however, some indirect evidence suggests that interstellar dust grains could be highly porous. While models of scattering by dust typically assume optical properties corresponding to spherical or spheroidal dust grains of uniform composition, Mathis (1996) tested models using composite dust grains, whose volume could also be composed of some fraction of empty space,
Dust grains detected in the comae of comets may also be porous, based on analysis of the range of grain sizes and densities. In the case of comet 67P, these grains appear to be loose aggregations of low density (0.1–1 g cm−3) and mean size around 10 μm that are made up of smaller sub-units (Agarwal et al., 2007). The sub-units themselves are of sizes more appropriate to interstellar grains, on average around 0.1 μm, and are thought to be closer in density to the native material from which they are composed (Mannel et al., 2019), indicating a low degree of porosity on those size scales. It is unclear, however, whether those sub-units could themselves be formed from smaller structures.
Dust-grain coagulation models appropriate to star- and planet-forming environments (Ossenkopf, 1993; Ormel et al., 2009) indicate that large, micron-scale grains can become porous, but their smallest sub-units are of similar size to the typical interstellar dust grains, and therefore do not indicate whether any porosity may exist in those smaller grains. Other dust models (Jones et al., 2014) have considered the accretion of carbon atoms onto existing grains in the outer regions of molecular clouds, as a mechanism to replenish the dust destroyed through energetic processes such as interstellar shocks, in order to explain the apparent imbalance between dust formation and destruction in the ISM. This accretion process, which would necessarily occur at low dust temperatures (no greater than around 30 K under typical interstellar conditions), could plausibly introduce irregularities to the grain surfaces, or, with a sufficient degree of accretion, outright porosity in the outer layers of the grain material.
Most recently, Potapov et al. (2020) carried out experiments in which interstellar dust-grain analog particles composed of amorphous carbon were formed and deposited onto a substrate. Ices of water, CO, and other molecular species were then deposited onto the grains, and their temperature-programmed desorption profiles analyzed. The authors reported behavior indicating that the ices were spread over very large surface areas, indicating an extreme degree of porosity. They proposed that interstellar grains may have surface areas enhanced by factors of several hundred, producing maximal ice thicknesses on the order of 1 monolayer in interstellar clouds. However, we note that the degree of porosity in interstellar grains themselves would be dependent on the precise mechanism and environment of their formation.
Here, we investigate the influence of interstellar dust porosity on the chemistry and composition of ices that form on the dust-grain surfaces. Dust porosity has the potential to affect the chemistry occurring on dust-grain surfaces, via increased surface area, enhanced local binding energies, and the possibility of trapping of molecules within the pores as ice mantles build up on the grains.
Garrod (2013a), mentioned above, introduced a microscopic Monte Carlo chemical kinetic model for interstellar grain-surface chemistry that was capable of simulating diffusive chemistry on fully three-dimensional surfaces. The model uses an off-lattice treatment, allowing surface particles to assume positions that are determined entirely by their interactions with nearby binding partners. In that study, the model was used to simulate the production of water and several associated molecules over astronomical timescales; an important effect noted by Garrod was the formation of porosity in the ices, which was dependent on both the temperature of the dust grain and the rate of accretion of gas-phase atoms onto the grain, as controlled by the gas density; under the dark-cloud conditions used in those models, the ices produced were found to be non-porous. Here, we use an updated version of the model with a larger chemical network, such that all of the main interstellar ice components may be modeled. The structure of the underlying dust grain may be directly controlled by the placement of the (carbon) atoms from which it is constructed, meaning that any structure can in principle be used; this includes the possible variation of both the large-scale morphology of a grain and its microscopic surface roughness. We run chemical models on both porous and non-porous grains of a selection of sizes, using broadly spherical large-scale structures. These models represent the first such attempt to model explicitly the effects of dust-grain porosity on interstellar grain-surface chemistry.
2 Methods
The model results presented here make use of the off-lattice microscopic Monte Carlo chemical kinetics code MIMICK (Model For Interstellar Monte Carlo Ice Chemical Kinetics), developed by Garrod (2013a). This model uses an exact description of a dust grain in which the positions of all of its constituent atoms are explicitly described. The model traces the chemistry occurring on the dust-grain surface by similarly tracking the positions and movements of the atoms and molecules that are adsorbed onto the surface from the gas phase. The positions of those species on the surface correspond to local potential minima (binding sites), as determined by the strength of the pairwise potentials experienced by each surface particle. Positions may change when mobile species make thermal hops between adjacent binding sites.
All of the various physical and chemical processes occurring on the grain, namely adsorption/accretion, thermal desorption, thermal hopping and reaction, are treated as a sequence of consecutive events. The Monte Carlo treatment used here follows that of Gillespie (1976), in which random numbers are used to determine 1) which process occurs next and 2) the time increment between each such event, based on the calculated rates of every possible process at that moment in time. In the case of thermal diffusion, each individual hop from one binding site to another is treated as a distinct event. Standard rate formulations are used for each of the various processes for each individual atom/molecule on the surface; thus, those rates are dependent not only on the type of atom/molecule in question but on the number, type, and spatial arrangement of its binding partners.
The model employed here uses a somewhat more advanced treatment for the calculation of surface binding energies and diffusion barriers than that of Garrod (2013a). Clements et al. (2018) provide an outline of this treatment, which is described in detail by Garrod and Deselm (in prep). In brief, each pairwise interaction between atoms, molecules or surface particles is treated as a Lennard-Jones (6–12) van der Waals potential, with a defined optimal separation, σ, and a potential well depth, ϵ. When an atom or molecule arrives on the surface through adsorption from the gas phase, its position is optimized by minimization of the sum of all pairwise potentials that it experiences (while the positions of the individual binding partners are assumed to be fixed). The 6–12 potentials are tuned to have a cutoff separation of
The process of diffusion may occur in multiple possible directions for a given surface atom/molecule, with each direction having its own associated energy barrier. Firstly, the model determines which directions are possible. To do this, all viable diffusion pathways are assumed to pass through saddle-points in the local surface potential; due to the spherical symmetry of the pairwise potentials, each saddle-point falls approximately between adjacent pairs of binding partners that are immediately contiguous to the diffuser (see Garrod, 2013a, for a diagram of such an arrangement). Diffusion through that saddle-point may be imagined as a rotation of the diffuser about the axis that joins two of those contiguous binding partners to each other; those binding partners are termed the “saddle-point pair.” Each saddle-point pair can be identified geometrically, and has a unique diffusion pathway associated with it. (Garrod, 2013a, also provides a more detailed description of the geometric considerations involved in identifying the saddle-point pairs). As an example, one may imagine an atom sited atop a triangular arrangement of three contiguous binding partners. Any given pair of those binding partners would constitute a saddle-point pair, with rotation around each such pair providing a unique direction of diffusion through a potential saddle point.
For every atom or molecule on the surface at a given moment, the model holds a list of saddle-point pairs (i.e., diffusion directions), for each of which a diffusion barrier is also defined. A diffusion barrier consists of the energy required to break or stretch each surface bond to the furthest distance the diffuser could reach as a result of rotation around the axis defined by the saddle-point pair. Some existing bonds would be entirely broken by such a procedure, while others would merely be stretched, while still remaining within the cut-off distance of the van der Waals interaction. The calculation excludes any bonds that would be shortened in this procedure. The bonds between the diffuser and the binding partners that comprise the saddle-point pair also do not contribute to the diffusion barrier, as the rotation around the saddle-point pair occurs at constant radial distance from each. The diffusion barriers are therefore always lower than the total binding energy, and they vary depending on the direction of diffusion; diffusion events that require the breaking or stretching of stronger pairwise potentials have higher barriers and are therefore less likely to occur. (The key difference between the present treatment and that of Garrod (2013a), is that the older model only considered contiguous binding partners as contributors to the diffusion barriers and binding energies). If and when a particular diffusion event is selected by the Monte Carlo routine according to its barrier-defined rate, the diffusing particle is allowed to rotate around the saddle-point pair until it encounters another particle on the other side; its position is then re-optimized based on all of the pairwise potentials experienced in its new location.
The full diffusion procedure described above, while being approximate in the determination of diffusion barriers and prescriptive in the choice of available directions, obviates the need to evaluate directly the potential energy surface involved in every possible diffusion event of every surface species; such a simplification is essential to allow this kinetic model to run on manageable timescales.
All surface diffusion considered in this model is based on thermal rates (as defined by the barriers described above); no tunneling contributions to diffusion rates are considered. At the 10 K temperature used in the models presented here, tunneling is not expected to contribute significantly to the diffusion rates of hydrogen atoms; quantum calculations by Senevirathne et al. (2017) indicate a crossing temperature (i.e., the temperature at which tunneling begins to dominate) of TC (H) = 9.7 K for H atoms on amorphous water ice.
In the present version of the MIMICK model, all surface particles, whether atoms or molecules (or atoms comprising the dust grain itself) are considered to be spherical, with a radius of 1.6 Å, such that the optimal separation used in the pairwise potentials is σ = 3.2 Å. This size is chosen to ensure that the physical size of the ice formed on the grains is representative of water ice, which is the dominant constituent. Variable values of σ will be introduced in future work.
Importantly, as in the models of Clements et al. (2018), the thermal desorption of a surface particle follows a specific trajectory, which is defined by the vector sum of the pairwise potentials at the starting position. The desorbing particle incrementally traces out this trajectory until it reaches the edge of the calculation-space used in the model. If the trajectory intersects with other particles then the desorbing particle may be re-accreted. This eventuality is especially important in cases where enclosed surfaces exist within the porous structures of the grain. Note that the present model does not use the treatment presented by Clements et al. for the gradual relaxation of a newly-adsorbed particle. Atoms and molecules accreted from the gas phase are assumed to thermalize immediately on the surface.
Only surface particles that are considered “free” are allowed to diffuse or to desorb. Free particles are defined as those that have available at least one surface diffusion pathway, as determined by their local geometry. Particles that are completely enclosed by contiguous binding partners (as would be the case for molecules within the bulk ice) have no such pathways. Thus, any particle that is capable of diffusion also has some probability of thermal desorption. Counting the number of free particles also provides a determination of the number of particles present in the surface layer of the ice.
Reactions between surface species are allowed to occur when the reactants become contiguous, whether that be through diffusion, adsorption, or as the result of an earlier reaction. Once contiguous, that reaction between specific particles is considered as a process that competes with all others to be the next event in the sequence. If the reaction has no activation energy barrier, then it is assigned a rate equal to the sum of the characteristic vibrational frequencies of the reactants. The rates of reactions that have activation energy barriers are further moderated by an efficiency factor representing the faster of either the thermal reaction probability (Boltzmann factor) or a probability based on tunneling through a rectangular barrier, as is typically employed in the usual rate-equation based models (e.g., Hasegawa et al., 1992).
Many typical (non-MC) kinetic models of grain-surface chemistry include a chemical desorption (CD) mechanism, whereby the energy released in the formation of a molecule via a surface reaction is sufficient to allow the molecule–surface bond to be broken and the reaction product to desorb into the gas phase (Garrod et al., 2007). This mechanism is usually assumed to occur with an efficiency of approximately (or sometimes precisely) 1%; in the remainder of cases, the reaction energy would be lost through vibrational interactions with the surface and the product would remain on the grain. Experimental studies (Minissale et al., 2016) have confirmed that this process can occur, typically with greater probability on bare-grain analogs than on the amorphous water ice that is expected to build up on the grain surfaces. For most reactions, those experiments obtained only upper limits to the efficiency, on the order of 10%. Recent ab initio molecular dynamics studies of the H + CO → HCO reaction system on crystalline water ice indicate that energy loss to the surface occurs too rapidly to allow product desorption (Pantaleone et al., 2020). In most kinetic models that include the CD mechanism, while its inclusion may be important to the abundances of certain gas-phase species (such as methanol) that do not apparently have a viable gas-phase formation mechanism, the influence on the abundances of surface species is relatively minor. We therefore leave the inclusion of this desorption mechanism in the MIMICK model to future work.
2.1 Construction of the Dust Grains
Three sizes of non-porous grain are constructed for use in the chemical model, with radii of approximately 50, 75, and 100 Å. In each case, this begins with the generation of a spherical grain composed of atoms in a simple cubic arrangement, placed at their optimal separations, σ. Because the chemical model MIMICK requires only the surface of the grain to be defined, a sphere of material is removed from the center, leaving a shell that is substantially thicker than the cutoff distance of the potentials used by MIMICK. The process up to this point is identical to that used by Garrod (2013a), albeit with much larger grains in the present case.
To produce an amorphous surface structure, the simple-cubic grain is first compressed in all three dimensions by a factor 0.75. The position of each atom is then offset in a randomized direction, by a random distance ranging from 0–σ. These two procedures randomize the positions of the atoms sufficiently to remove all trace of the original simple cubic micro-structure, while retaining the spherical shape. Each atom in the grain is then allowed to relax, according to a randomized sequence, using Lennard-Jones type potentials of arbitrary, uniform strength. This ensures that all of the atoms in the grain are at appropriate distances from each other. The upper row of Figure 1 shows the non-porous grains produced in this way; the surface of each is amorphous, while retaining a well defined spherical morphology at large scales.
FIGURE 1. Ray-traced images of the five grain types used in the chemical models. All grains are composed of a collection of spheres representing carbon atoms. Top row: the three non-porous grains; middle row: the two porous grains, of comparable size to the two largest non-porous grains; bottom row: the cross sections of the two porous grains. From left to right, grains have a radius approximately 50, 75, and 100 Å. Note that the inner portions of the grains are removed to simplify the chemical kinetics calculations; the inner void is completely inaccessible to surface atoms and molecules.
To produce a porous grain, a smaller grain core is initially created using the method described above. The largest porous grain uses the smallest non-porous grain as its core, while the medium-sized porous grain uses a core of radius 25 Å. The chemical model is then employed to deposit more (carbon) atoms from the gas phase onto this core surface. The trajectories of the incoming atoms are entirely randomized. Upon landing, no chemistry or surface diffusion is allowed, but the positions of the accreting atoms are optimized, again using Lennard-Jones potentials. This hit-and-stick type approach results in open, porous structures. The accretion of atoms is halted when the maximum radial position of any atom (centered on the zero-position of the grain core) is equal to the maximum radial position of atoms in either the intermediate- or the large-sized non-porous grain; the resulting porous portions of each grain have a maximum thickness ∼50 Å, but the average is closer to 40 Å. The porosity of the porous portions of both grains is measured to be ∼64%. Taking into account the non-porous core of each grain, this produces overall porosities of 61 and 53% for the medium and large grains, respectively. The latter is not far from the value used by Mathis (1996).
The five different grains created for the chemical simulations range in size from 50–100 Å in radius. This is smaller than the typical value of 0.1 μm (1,000 Å) used in chemical models, although smaller grains such as these would be more numerous in the distribution and should provide a larger fraction of the available surface area for accretion, assuming non-porous, spherical grains. Smaller grains are also more computationally favorable using the microscopic Monte Carlo kinetics approach, as larger grains generate larger ices composed of greater numbers of atoms and molecules, each of which must be handled individually.
During the main chemical simulations, the atoms comprising the dust grains are fixed in place and do not participate in any chemical reactions, either directly or as a catalyst. Their main function is to define the physical surface of the grain.
2.2 Physical and Chemical Conditions
To represent a dark interstellar cloud (or, more precisely, a dark cloud core), the gas density is set to nH = 2 × 104 cm−3 and both the grain and gas temperatures are set to 10 K (e.g., Bergin and Tafalla., 2007); similar conditions have been used in past models of dark cloud chemistry (e.g., Hasegawa et al., 1992; Garrod et al., 2007). Although no gas-phase chemistry is included in these models, the gas density and temperature determine the rates of accretion of gas-phase species onto the grain.
The gas-to-dust ratio by number is tailored for the specific grain size used in each model run, to maintain a consistent mass ratio of 100:1 appropriate to the Milky Way. This influences the total number of atoms initially available in the gas phase in our models. The chosen grain size used in each model run is assumed to be representative of the entire dust population, assuming a mass density of 3 g cm−3 for grain material. As grain size increases, the number of representative dust grains must decrease, thus increasing the amount of gas per dust grain. For the large grain, a gas-to-dust ratio by number of 1010 is used, which corresponds specifically to the number of hydrogen atoms present in either H or H2 at the start of the models. The medium and small grains use lower values, as scaled by volume. For comparability between the chemical models, the gas-to-dust ratios used for the large- and medium-sized porous grains are the same as those of their non-porous equivalents.
The total budget (i.e., number of atoms) of each element used in the model is generated based on the gas-to-dust ratio and the abundance of each, relative to total hydrogen. The assumed gas-phase composition at the beginning of each model is shown in Table 1, based on the initial gas-phase elemental abundances used by Garrod (2013b). The grains in our models begin with no ice mantles. To better represent the composition of ices in dark clouds, a proportion of the carbon and oxygen budgets is incorporated into CO. Most of the hydrogen budget is assumed to be in the form of molecular hydrogen, with the atomic H abundance also chosen to represent typical dark cloud values. Gas-phase H2 is not allowed to adsorb onto the grains, following Garrod (2013a); this simplification is necessary to allow the models to run on manageable timescales.
TABLE 1. Initial gas-phase abundances as a fraction of total hydrogen.
As the models progress, atoms and molecules in the gas accrete onto the grain surface, depleting the gas-phase budget of each and thus reducing the amount still available for the further growth of the ice mantles. However, an exception is made in the treatment of the atomic H budget. To ensure that its rate of adsorption is always appropriate to dark cloud values, its abundance in the gas phase is held steady at its initial value. In full gas-grain models of dark cloud chemistry, gas-phase processes such as cosmic-ray induced dissociation allow the gradual release of some hydrogen locked up in H2 back into atomic form, meaning that the fractional abundance of H is maintained at some base value, even though other species may become strongly depleted from the gas phase over time. This simple fix to the code avoids having explicitly to model the gas-phase chemistry.
The model assumes that the sticking coefficient for all species is unity, following typical assumptions used in other models. Experimental work by He et al. (2016b) indeed found a sticking coefficent of 1 for CO (the only adsorbing molecule in our model), although their low-temperature value for H2 was in the 0.6–0.7 range. The sticking coefficient for atomic H on an amorphous water surface was calculated to be 0.980 by Dupuy et al. (2016) at dust/gas temperatures of 10 K. Sticking coefficients for H on amorphous carbon surfaces are not well defined, although Pirronello et al. (1997) suggest it could be appreciably less than one. Experiments by Chaabouni et al. (2012) find a value of 1 for silicate surfaces, however.
The model includes 39 chemical species and 70 different reactions, allowing C, O, N and CO to be hydrogenated all the way to CH4, H2O, NH3 and CH3OH. Molecules such as H2, C2, O2, N2 and CO may also be formed on the grain surfaces, along with various others. Table 2 details each of the reactions included in the model, with the assumed activation energy barriers and barrier widths indicated; this network is a simplified version of that presented by Garrod (2013b). Most of the reactions occurring on the grains are initiated by the diffusion of atomic H to meet its reaction partners, although reactions may occur in any situation in which two reactants come into contact. Table 3 provides a selection of important pair-wise potentials used in the model, which determine both the binding energies and diffusion barriers experienced locally by surface species. The potentials are somewhat lower than those used by Garrod (2013a), which reflects the fact that contributions to the binding of all local species within the cutoff distance are now considered, not only those species that are contiguous to the particle in question, which was the previous method. As described by Garrod (2013a), the pairwise potentials are intended to reproduce the overall binding energies of the species in question on surfaces of, for example, pure amorphous water.
TABLE 2. The reaction network used in MIMICK. Activation barriers, where present, are in units of Kelvin, with barrier widths in Angstrom.
TABLE 3. Pairwise potentials, ϵ (K), between selected surface species.
The barriers against diffusion in any direction are always lower than the corresponding binding energies for a particular particle in a particular binding site, although different directions of diffusion typically have different values. The models produce as an output the distribution of the ratios of diffusion barrier vs. binding energy for each chemical species in the network (in increments of 0.05), which are calculated at every output time in the model. In around 99% of atomic-H surface diffusion events, the barrier of the successful diffusion direction lies in the range 0.25–0.7 times the local binding energy, with around 18–20% of those barriers falling in the 0.4–0.45 range, which is the most populous bin. The average binding energy of an H atom on the final (outer) surface of the ice mantles is approximately 460 K. This compares well with the value used in many past gas-grain models (e.g., Garrod, 2013b) of 450 K; DFT calculations by Wakelam et al. (2017b) indicate values of 400 K and 680 K for H on water, dependent on the method used. Using our model’s average binding energy of 460 K, the diffusion barriers in the above fractional range produced by the model provide absolute values in the range 115–322 K, with the most common falling in the 184–207 K bin. Senevirathne et al. (2017), based on their calculations for H on amorphous water ice, recommended a representative ratio of diffusion barrier to binding energy for H of 0.37; combined with their representative binding energy of 661 K (57 meV), this suggests a diffusion barrier of 243 K. This latter value in particular seems fairly consistent with our range of values, in spite of the higher overall binding energy. Indeed, looking at the actual distribution of diffusion barriers produced by those authors, their most common diffusion barriers lie in the ∼175–230 K range, which coincides well with our most common bin.
For H2, the model produces average binding energies of around 380–395 K; Garrod (2013b) use a value of 430 K. The calculations of Wakelam et al. (2017b) for the H2–H2O dimer suggest a binding energy of 800 K, but full cluster calculations by Ferrero et al. (2020) for amorphous solid water give a value in the range 226–431 K, in keeping with our model values. The experimental measurements of He et al. (2016a) indicate binding energies in a range between 322 and 505 K, with the former corresponding to high H2 coverage and the latter to low coverage. Our models produce binding energies for CO and O2 of around 1015 K and 1240 K, respectively; these are somewhat higher and somewhat lower, respectively, than the calculations by Ferrero et al. (2020) although our values fall comfortably within the range of experimental values quoted by those authors (870–1600 K for CO and 914–1520 K for O2).
For H2O, the model produces an average binding energy around 4,750 K. The outer ice surface contains a selection of molecules besides water, so the effective binding energy of a water molecule is somewhat lower than might be expected on a pure water ice; a common value used in rate-equation models is 5,700 K (Garrod, 2013b). The calculations of Wakelam et al. (2017a) indicate values of 4,800 or 5,300 K based on water dimers or clusters, while Ferrero et al. (2020) obtain a range of 3,605–6,111 K. Experiment also provides a range of values, although Sandford and Allamandola (1990) obtain 4815 K for unannealed water ice that is very close to our value. Again, we note that our value is produced with a water-dominated ice mixture that includes all of the typically observed interstellar solid-phase species, rather than pure water. In practice, at the low temperatures used in the present models, the rates of diffusion and desorption of water molecules are negligible.
The model does not consider the influence of the condensation energies of adsorbed species, the energy released by chemical reactions, nor the effects of high-energy photons on the overall temperature of the grain/ice. Although the associated temperature fluctuations could be substantial for very small grains (e.g., Cuppen et al., 2006), they are beyond the scope of the present model. The recent modeling simulation of water-ice deposition by Clements et al. (2018) included a treatment for the non-thermal diffusion of newly-adsorbed molecules, caused by the excess energy gained as they are accelerated into the surface potential. That treatment was tuned specifically for water molecules, and we do not extend it here to other species. In any case, the ability of the adsorbed atoms in the present model to exhibit some degree of thermal diffusion following adsorption makes this simplification less important under the conditions tested here.
The models are run for simulated times on the order of tens of thousands of years, and by the end of each model run around half of the heavy atoms have been depleted from the gas phase, forming ices on the grain surfaces. By the end of its run-time, each model is thus on the cusp of experiencing strong gas-phase depletion which would slow the further growth of the ice mantle. The time taken to build up an ice mantle of a particular size is dependent on the cross-sectional area presented by the combined grain and ice mantle, averaged over all solid angles (the details of the accretion process are described in detail by Garrod, 2013a). The different sizes and morphologies of the grains thus affect the cross-sectional area presented to gas phase species, while those cross sections also rise as the ice mantles grow over time. For purposes of comparison between models, it is therefore helpful to define output times in relative terms, corresponding to the amount of depletion from the gas phase. This measure is more instructive than, for example, a determination of the number of ice molecules on the grains or the thickness of the ice, as it makes no assumptions about the chemical or physical form that the ice mantles take. Each model is set up to produce outputs for 1,000 different time values, which here we will refer to as “intervals.” Those intervals are evenly spaced in the amount of O, C, N and CO that have been accreted from the gas phase as a fraction of the initial amount (which varies depending on the size of the grain, via the gas-to-dust ratio). Each interval corresponds to a roughly equivalent amount of ice build-up, but as time goes on and depletion from the gas phase increases, the timescale required to complete one interval lengthens. The final populations of the grain-surface molecules achieved in each grain-size model scale roughly with the volume ratio of the underlying grains.
The run-time (wall time) of the models ranges from a few days to two weeks, depending on the size of the dust grain.
For each grain size and morphology used in the models, a total of three physically and chemically identical runs are carried out, with the only difference being the initial random number seed used in the Monte Carlo algorithm. This ensures that the model results are reproducible, and to some extent demonstrates the degree of random variation in the chemical abundances.
3 Results
Figure 2 shows ray-traced images of the cross-sections through the large porous and non-porous grains, at selected points in the model runs. From left to right, the images depict the growth of the ice mantles at comparable moments (based on the degree of gas-phase depletion of heavy species). For the non-porous grain, the material from the gas phase builds up relatively uniformly over the entire surface. There is also no indication of empty pockets. Based on the results of the hydrogen/oxygen system modeled by Garrod (2013a), the ice produced on dust grains under dark-cloud physical conditions should be non-porous, as heavy atoms accreted from the gas phase would have ample time to diffuse on the surface, finding and filling any nascent surface inhomogeneities, which tend to manifest initially as stronger binding sites due to their local curvature.
FIGURE 2. Ray-traced, cross-sectional images of the large (100 Å radius) non-porous (top) and large porous (bottom) grain types at different stages of ice build-up; left to right: interval 0, interval 100, interval 500, interval 1,000. Gray spheres represent the carbon atoms composing the dust grain. Grain-surface atoms are represented by spheres of white (H), red (O), black (C) and blue (N). Molecules are similarly represented as collections of such atoms in appropriate structures. H2 molecules are colored yellow for ease of identification.
For the porous dust grains, after a sufficient amount of ice mantle has been built up (interval 100), the surface appears to be relatively uniform as well. However, looking below the surface, there are some differences. Firstly, rather than filling up the pores in the grain material, the ice fills them only partially, with some porosity remaining as the ice builds up over the top, leaving empty pockets that are not present in the non-porous grains. Secondly, there appears to be a collection of molecular hydrogen gathered in those remaining porous structures.
Figure 3 shows cross-sectional close-up images of the large porous grain at early times in the ice-mantle evolution. Porous structures are gradually filled, sometimes leaving voids in the ice. While much of the volume of the pore voids present in the original grain is filled, some of the bare grain surface remains and there are internal ice surfaces formed where the pores are not completely filled. H2 molecules, marked in yellow, can also be seen more clearly in Figure 3; some are apparently bound directly to the grain material itself, while others are embedded in the ice. Notably, the H2 molecules present at interval 25 (left-most images) seem mostly to retain their positions at later intervals, indicating that they are well fixed into those structures. It may be noted that the ices that form once the pores in the grain material are covered over are themselves non-porous; as found by Garrod (2013a), under dark-cloud conditions, oxygen and other atoms have long enough pre-reaction lifetimes before they react with the more mobile H atoms (producing e.g., water), that they may diffuse enough to find strong surface binding sites, leading to a smoothing of the ice surface. The porosity that is present in the ice by the end of the model runs is due to the initial presence of porosity in the grain material. The final grain/ice does not retain that same degree of porosity, as the ice partially fills the original voids, and then forms non-porous ice once those are fully enclosed.
FIGURE 3. Ray-traced, cross-sectional images of the large (100 Å radius) porous grain type at the early stages of ice build-up, corresponding to the grain shown in the bottom row of Figure 2; here, the cross-section at a different angle is shown. Left to right: interval 25, interval 50, interval 100. The three upper panels show the full cross-section, while the lower panels show close-up images of the same small region of the grain at each of the three intervals. Gray spheres represent the carbon atoms composing the dust grain. Grain-surface atoms are represented by spheres of white (H), red (O), black (C) and blue (N). Molecules are similarly represented as collections of such atoms in appropriate structures. H2 molecules are colored yellow for ease of identification.
Table 4 shows the absolute number of atoms and molecules resident on the grain surfaces at the end of each model. Values shown are the mean final populations of the three random number-seed runs for each setup. The final time in each model is also indicated. Table 5 shows the same population data given as a fraction of the total water on the grains, which allows an easier comparison between models. The porous grains take a little longer to reach their final state (i.e., the same amount of ice) as, compared with their non-porous equivalents, the porous grains present slightly smaller accretion cross sections.
TABLE 4. The total abundance of some notable ice species averaged over three runs. The final time of each model is also indicated, which corresponds to interval 1,000 in each case.
TABLE 5. The fractional abundance of some notable ice species relative to the total water ice, averaged over three runs.
Comparing between models, firstly we may consider the effect of grain size, by contrasting the relative abundances achieved in each of the three non-porous grain models. The abundance data suggest that variation in size has little effect on the relative abundances of the main ice constituents, such as CO, CO2, and other stable molecules. CO reaches a population a little over 4% of total water, while CO2 is at around 49%. In comparison to observational values toward background stars, the latter is a little high, while the former is a little low; observational lower and upper values are 9–67% and 14–43%, respectively (Boogert et al., 2015). This may be explained by our use of static physical conditions throughout the model run, combined with the lack of gas-phase chemistry. The initial chemical abundance of CO is representative of a somewhat evolved gas-phase chemistry, while that of atomic oxygen is more appropriate to early times; this combination therefore leads to a somewhat higher degree of conversion of CO to CO2 than observations would suggest.
The fractional abundances of volatile species on the grains do not scale so well with grain size; C, O, N and OH are seen to fall with increasing grain size. The absolute abundances of these species increase with greater grain radii, but not at the r3 ratio seen for stable species. However, these values are likely reflective of the lifetime of such species against reaction with a newly-accreted H atom on the grain surface, which is more dependent on grain surface area than on total volume. The differences between grain-size models do not otherwise seem to affect the resultant chemical compositions of stable molecules.
The two porous grains may also be compared. Fractional abundances between the medium- and large-sized porous grains are again similar for stable species, while the atoms and radicals show somewhat less variation; indeed the fractional abundance of atomic H is identical between the two. We note also that the quantities of atomic H shown in the Table correspond largely to H trapped within the bulk ice; the total, instantaneous population of free (surface) H atoms on the grains in any of the models hovers around one.
Comparison of the abundances of each chemical species between the porous and non-porous models of the same grain size demonstrates an interesting outcome. Beginning with the medium grains, the final fractional abundances are fairly consistent with one another for most species, especially stable molecules. The one molecular abundance that stands out in comparison between porous and non-porous grains is that of molecular hydrogen, which is about two times greater in the porous cases than it is in the non-porous cases, by the end of the model runs. The absolute numbers of H2 on the grains are on the order of hundreds, so the differences cannot be due to random small-number fluctuations. Consideration of the large porous and non-porous grains shows similar behavior for all species, while the discrepancy in the absolute abundance of H2 is further exacerbated. Atomic H also shows a much stronger variation than the other atoms when comparing between porous and non-porous grains.
In order to pinpoint the cause of this, the time-dependence of the H2 population may be examined. Figure 4 shows the absolute populations of H2 for the large porous and non-porous grains, averaged over the three random number-seed runs for each case; each of the individual runs is close to its mean curve. For most intervals, the porous and non-porous grains show almost identical trends, except that the porous-grain model has a very rapid up-tick at early times–times at which the porous grain-surface is still exposed. Ignoring this offset, there is otherwise convergence in behavior after that initial jump. It is evident that the enhancement in H2 population at those early stages becomes a permanent feature of the ice, as the offset is retained right until the end of the models. The growth in the H2 populations toward the end of both sets of models is caused by the falling rates of accretion of heavy species from the gas phase, as their limited budgets run out, while H atoms remain plentiful throughout, increasing the probability of H2 production (the same effect is seen in Figure 5).
FIGURE 4. The amount of molecular hydrogen present on the grain vs. relative time intervals, for the large porous and non-porous grain models. The solid lines show a moving average (of period 10) over the three runs for either grain type, with the data from the individual runs plotted as individual points. Interval 1,000 corresponds to the end time of the models shown, which are shown individually for each model in Table 4.
FIGURE 5. The total count of hydrogen molecules formed in the model up to a given moment, vs. relative time intervals. All H2 formed is shown, regardless of whether it is trapped in the ice, desorbed into the gas phase, or destroyed by chemical reactions. Data from the three runs for the large porous and non-porous grain models are shown. Interval 1,000 corresponds to the end time of the models shown, which are shown individually for each model in Table 4.
There are several possible explanations for the larger population of H2 on the porous grains. We propose the following as the most likely: 1) the porous surface enhances the overall conversion of H into H2, with mobile H2 formed in the pores diffusing away to be incorporated into the general bulk ices; 2) the porous surface enhances production rates locally within the pores, leading to the build-up of H2 where it is formed; 3) H2 that is formed anywhere on the grain surface may diffuse into the pores, where it becomes trapped.
Figure 5 shows the total number of H2 molecules formed throughout each model; the plot shows all of the three runs for both the large porous and large non-porous grains. The overall production of H2 is extremely similar between models and grain sizes, including at the early times when the effects of the porosity manifest. Any difference may be accounted for by slight variations in the sizes of the porous and non-porous dust grains. We can therefore rule out a difference in the production efficiency of H2 on the porous grain as the cause of its greater H2 retention.
It is instructive also to consider the amount of “free” H2 throughout each model, shown in Figure 6 (large grains). This is a measure of the number of H2 molecules that are capable of undergoing diffusion or desorption, at whatever rate, and which may therefore be classified specifically as surface species. Molecules that are not “free” are necessarily trapped under other species. The behavior of the free H2 is qualitatively similar to that of the total H2 population, and it also grows strongly at early times in the porous grain models. The absolute number of free H2 molecules in the porous model reaches a peak of around 45 until late times, which is a factor ∼15 lower than the total H2 on the grain at that time. Over the same period in the non-porous case, there is only a few free H2 molecules at any one time, corresponding to roughly one in fifty of the total H2 on that grain.
FIGURE 6. The amount of free molecular hydrogen present on the grain plotted vs. relative time interval, for the large porous and non-porous grain models. Free species are able to undergo surface diffusion and desorption, but may reside on enclosed surfaces within pores. The solid lines show a moving average (of period 10) over the three runs for either grain type, with the data from the individual runs plotted as individual points. Interval 1,000 corresponds to the end time of the models shown.
Clearly, the much greater absolute numbers of surface H2 on the porous grain indicates that those particular molecules at least are stored on the surfaces of the open porous structures that become closed off to outside chemical influence by around interval 100. With such a build-up in surface H2 abundances in the pores, it is natural that some H2 should become incorporated into the ice mantles that form specifically within the pores (i.e., no longer free). Figure 7 shows an image of the H2 present on the large porous grain at the end of one of the runs, with no other atoms, molecules, or grain material depicted. The image includes all H2 molecules present in three dimensions, not just a cross-section. The distribution of H2 clearly favors the pore structures, although some H2 is still present in the outer mantles.
FIGURE 7. Ray-traced image of the large porous grain type at the final stage of ice build-up showing only the H2 molecules and excluding the grain and all other ice species.
We again consider the feature seen in Figures 4–6 in which the amount of H2 (whether the amount on the grain, the number of free molecules, or the total formed) increases over the last few hundred intervals. This effect is the result of the run-down in the gas-phase budget of species other than hydrogen, i.e. the depletion of those species. While, for example, atomic O becomes less and less abundant in the gas phase over time, due to its adsorption and incorporation into the ice mantle, the amount of gas-phase H stays the same; this means that for any one adsorbed H atom, its surface reaction partner is increasingly likely to be another H atom, which drives up the H2 production rate. Meanwhile, the time required for a fixed number of gas-phase atoms (or CO) to be adsorbed from the gas phase also rises, due to their smaller remaining budgets. The amount of time associated with each interval therefore lengthens (see Section 2.2), making the increase in H2 molecules appear more dramatic.
Figure 8 shows the amount of free H2 on the grain, just as in Figure 6, but is shown here as a percentage of the total population of free species. Note that the number of free species in this model corresponds approximately to the number of species (atoms or molecules) present in the surface of the ice or on the surface of the grain (during the early stages before a mantle has built up). This count of surface species includes not only the outer surface but any enclosed surfaces that remain within the pores. As a fraction of the total, we see initially that the difference between the porous and non-porous grains is quite stark, as the pores in the grain material retain a substantial amount of H2. The difference gently declines for later intervals for two reasons: 1) much of the porous structure within the grain material becomes filled up and smoothed out, so that some of the surface H2 trapped there is locked into the bulk ice that resides in the grain pores, and 2) the outer surface of the ice mantle grows, so that the enclosed surfaces that remain deep down in the grain/ice account for a smaller fraction of the total surface available. By the end of the two runs, the percentage of surface coverage by H2 is similar (within a factor of 2); the difference that remains nevertheless corresponds to the H2 that is trapped within porous structures, residing on their enclosed surfaces and unable to escape. Plotted as a percentage of total coverage, the apparently dramatic rise in free H2 (seen in Figure 6) is far more subdued in both models.
FIGURE 8. The percentage of free molecular hydrogen present on the grain as a function of all free species, plotted vs. relative time interval for the large porous and non-porous grain models. Free species are able to undergo surface diffusion and desorption, but may reside on enclosed surfaces within pores. The solid lines show a moving average (of period 10) over the three runs for either grain type, with the data from the individual runs plotted as individual points. Interval 1,000 corresponds to the end time of the models shown.
4 Discussion
Aside from some small variations in fractional abundances, neither the grain size nor porosity appear to have a major effect on the grain-surface ice composition, except in the case of H2. Atomic H is also retained to a somewhat higher degree on the porous grains, but the absolute numbers are relatively small compared with H2.
Porous grains retain molecular hydrogen uniquely, either in the ice mantles that form within the pores, or on the surfaces of the pore mantles, in cases where voids exist within the icy pores (having been partially filled and then closed over by ice). The porosity does not induce any greater production of H2 (or other stable molecules) on the grains; rather it acts simply to trap those molecules. However, even if the total H2 produced is essentially unchanged, it is a different question to ask whether the presence of pores affects where on the grain that fixed amount of H2 is formed. Trapping of atomic hydrogen on pore surfaces could lead to local production of H2 that may likewise become trapped; alternatively, H2 produced elsewhere on the grain could also become trapped, following surface diffusion into the pores. Both are plausible explanations for the excess of trapped H2 in the pores and pore ices.
As described by Garrod (2013a) in the case of porous ices produced on smooth grains, the reason molecular hydrogen becomes trapped in pores at all has to do with the effect of surface structure on binding energies. The surface of a smooth grain is flatter; thus, fewer grain-surface binding partners maintain close proximity to a given hydrogen molecule at any time. However, the interior of the pores is more irregular and concave on very small scales, allowing more grain atoms or bulk-ice molecules to contribute to the binding energies and diffusion barriers experienced by a particle within that structure. The resulting stronger binding site allows for the molecular hydrogen to be held within the pore for long enough for other ice material to cover it up, or at least to seal over the pore before it is able to escape. The larger proportion of “free” H2 molecules in the porous, vs. non-porous, grain models indicates that the sealing of porous voids is a non-negligible effect, although the complete incorporation of excess H2 into the pore ices is the primary cause of retention.
Atomic hydrogen is mobile on the grain surfaces, but upon entering a pore it may encounter a strong binding site and be retained there for a longer period than if it were in a weak binding site. The arrival of another H atom would convert it to H2, allowing that molecule to be trapped in that site. It may be seen, though, that other molecules–whose progenitor-atoms such as O, C and N are much less mobile on the grains–can also fill up the pores; the pores are not filled solely or even predominantly with H2. Also, the overall ice structure outside of the grain-pores does not become porous itself once the underlying pores are covered over. This is because the heavier atoms in the model are still sufficiently mobile to find strong binding sites for themselves, including those located within the pores, prior to being hydrogenated by atomic H (which leads to a stronger surface bond and much lower mobility, at which point that species becomes fixed).
Why, then, should H2 be uniquely affected by the presence of many strong binding sites within the pores? H2 is unique also in being a mobile surface species that is largely unreactive. It may react indeed with OH, for example, but all of its few reactions have an activation energy barrier. On the other hand, atoms like H, O, C, and N, even though reasonably mobile, will react immediately upon meeting some other reactive species, ending their diffusion. Their diffusion path length on the surface is therefore limited by meetings with potential reaction partners, as well as by their somewhat slower diffusion rates (vs. their main reaction partner, atomic H). An H2 molecule, however, has the chance to sample a much larger fraction of the grain surface, seeking out strong binding sites, with its main loss route ultimately being desorption into the gas phase rather than reaction. As well as having a greater chance to find those strong sites, the lifetime against desorption of an H2 molecule will be substantially increased if it does find one.
Given a large selection of strong binding sites on the grain, due to their high mobility H2 molecules are most likely to fill them. The availability of a larger number of strong sites in the porous case therefore invites the retention of H2. When strongly bound species are formed on top of the H2, or in close quarters, that H2 becomes trapped. The disappearance of strong binding sites as the ices build up and become more uniform in structure removes this effect.
While the numbers of H2 being trapped in the present models are on the order of hundreds, larger grains would likely store proportionately more. Also, the model specifically leaves out the accretion of H2 directly from the gas phase; all of the trapped H2 in the models must first be formed on the grains through chemical reactions. In reality, grains would be under a constant flood of H2 adsorption, even if the individual lifetimes of those molecules were short. We should expect that any strong binding sites in the grain pores would be rapidly and comprehensively filled by H2 in that case. The potential for trapping could therefore be much higher in a more realistic simulation. Unfortunately, as the H2 abundance in the gas phase is around 1,000 times higher than the other species collectively, at current computational speeds those simulations would take years rather than weeks.
Similarly to molecular hydrogen, the trapping of other unreactive, volatile species such as helium could also be important. Interstellar He abundances are around 10% with respect to total H (by number), and could therefore similarly flood the surface. Neon, as the next inert element, might also be sufficiently mobile on the grain surfaces to find strong binding sites and therefore experience trapping. If these or other inert species could be trapped in the deepest layers of the ices, close to the grain surfaces themselves, then they could be retained on the grains long into the star formation process, even under warm conditions where much of the ice material might have already desorbed. Such a mechanism may be contrasted with the ice annealing process explored by, e.g., Bar-Nun et al. (1987), whereby highly-porous water ice was used to trap Argon gas, with the pores being closed over through a diffusive creeping mechanism of the ice at the pore openings.
It should be noted that it is specifically the microstructure of the porous grains that is most important to the retention of H2, or other volatiles similarly affected, in the models presented here. The coagulation of large, smooth grains would not likely produce such an effect. However, the presence of an icy coating on the grain pore surfaces may be just as effective as the bare grain itself in trapping H2, so long as the build-up of the ice within those pores allows the underlying microporous and/or irregular structure to be retained.
The release of the chemical energy of surface reactions has been suggested as a mechanism for both the general heating of a dust grain and the non-thermal desorption of H2 in particular (see e.g., the review of Wakelam et al., 2017a). Calculations of the heating effect of H2 production on silicate grains through the addition of H atoms by Navarro-Ruiz et al. (2014) suggested a strong influence on the temperatures of very small grains. Scaling their calculations to the size of our “large”, 100 Å radius grains, the expected change in temperature would be less than 0.1 K per H2 production event, with a time constant for the radiation of the excess energy of less than 1 μs, i.e., much faster than the rate of H2 production. The whole-grain effect of chemical heating (at least via H2 formation) is therefore unlikely to have a significant influence on our models, although the somewhat smaller local structures produced in our porous grains could be more strongly affected. Chemical heating could affect the local structure of the ice also; however, Oba et al. (2009) suggest that in the case of amorphous ice formation through surface chemical reactions, the structure is not influenced by this mechanism, in line with calculations by Kouchi et al. (1994).
The chemical energy of reactions could allow H2 formed on the grains to desorb through direct chemical desorption, although in our models the desorption process follows explicit trajectories, so that an H2 molecule desorbed in this way could certainly be re-adsorbed. However, if the surface of the grain is indeed flooded with H2 molecules (see above), the influence explicitly of the H2 formation process should be far less important than the diffusion of H2 into the pores. Indeed, the models presented here show that the H2 does not need to be formed in situ within the pore structures to become trapped there.
The kinetic models used here sample only a small region of parameter space in terms of grain size, porosity and small-scale structure/roughness. The simulated dust-grains are constructed crudely to produce high degrees of porosity. A more detailed determination of the precise microstructures of interstellar dust grains is therefore necessary to constrain the proposed effect using chemical kinetic models. It is also the case that the carbon atoms that constitute the grain material in our treatment are not allowed to participate in any way in the surface reactions; we do not, for example, consider chemisorption in the production of H2 or other species. As discussed by Wakelam et al. (2017a), reactions between physisorbed H atoms are likely to be the dominant H2 formation mechanism under translucent and dark cloud conditions, while substantial water or other ice mantles are not expected to build up before visual extinctions greater than around 3 magnitudes (i.e., translucent) are achieved (Whittet et al., 2001).
Our model also uses only a simple, static treatment of the physical conditions of dark clouds, ignoring the lower density and visual extinction that would pertain at early times, and which could lead to a somewhat higher dust temperature. The initial chemical conditions are also quite simple, with CO and H2 assumed already to have formed in substantial quantities at the time the models begin. The lack of a gas-phase chemistry could also affect the results; it is possible, for example, that the gas-phase production of water would lead to some component of the surface water ice originating through direct deposition from the gas phase; such conditions would be a little closer to those frequently employed in the laboratory for the study of porous amorphous water ices, and might therefore lead to a higher degree of porosity in the ices themselves that would not be associated directly with the grain material’s porosity. However, past dark-cloud gas-grain chemical models by our group (e.g., Garrod and Pauly, 2011) have not found any substantial contribution from gas-phase water in the build-up of the ice mantles, even under conditions of low extinction and elevated dust-grain temperatures.
The possibility of thermal desorption and re-deposition of water ice onto the dust grains in protoplanetary disks, however, might provide a regime in which ices of substantial porosity could be formed, through direct deposition rather than chemical formation on the grain surfaces. H2 or other volatiles might be expected to be trapped in such structures also.
5 Conclusion
Our chemical kinetic models of interstellar grain-surface chemistry indicate that, given a highly porous dust grain, molecular hydrogen may be trapped in strong binding sites caused by the irregular surface structure of the pore surfaces. The coating of those pore surfaces with water and other common dust-grain ice mantle constituents locks the H2 into the ice, although the irregular surface structures of the ice itself continue to allow trapping until the pores are completely filled and the grain/ice takes on a more regular, spherical morphology. The closing over of the pores, creating voids, is only a small part of the trapping effect, but a fraction of the trapped H2 is found to be free on the surfaces of those enclosed voids. A full-scale treatment of H2 adsorption onto the grains using this microscopic treatment might be expected to exhibit even greater degrees of trapping. The trapping effect could provide a means for other volatile and relatively inert species to be retained on the grains, enclosed in the deepest ice layers. Such an effect would allow chemical species trapped under interstellar cloud conditions to be preserved in the ices into much later stages in the star- and planet-formation process.
Data Availability Statement
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.
Author Contributions
DC ran computational models, analyzed the data, and co-wrote the manuscript. RG wrote the MIMICK computational code, helped analyze the data and co-wrote the manuscript.
Funding
This work was funded by the NASA Astrophysics Research and Analysis program, through grant number NNX15AG07G.
Conflict of Interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Acknowledgments
We thank the referees for comments and suggestions that improved this paper.
References
Accolla, M., Congiu, E., Manicò, G., Dulieu, F., Chaabouni, H., Lemaire, J. L., et al. (2013). Morphology of the solid water synthesized through the pathway D + O2 studied by the sensitive TPD technique. Month. Not. Royal Astro. Soc. 429, 3200–3206. doi:10.1093/mnras/sts578
Acharyya, K., Hassel, G. E., and Herbst, E. (2011). The effects of grain size and grain growth on the chemical evolution of cold dense clouds. Astrophys. J. 732, 15. doi:10.1088/0004-637X/732/2/73
Agarwal, J., Müller, M., and Grün, E. (2007). Dust environment modelling of comet 67p/churyumov-gerasimenko. Spa. Sci. Rev. 128, 79–131. doi:10.1007/s11214-006-9139-1
Bar-Nun, A., Dror, J., Kochavi, E., and Laufer, D. (1987). Amorphous water ice and its ability to trap gases, Phys. Rev. B. 35, 2427–2435. doi:10.1103/PhysRevB.35.2427
Bergin, E. A., and Tafalla, M. (2007). Cold dark clouds: the initial conditions for star formation, Annu. Rev. Astron. Astrophys. 45, 339–396. doi:10.1146/annurev.astro.45.071206.100404
Boogert, A. C. A., Gerakines, P. A., and Whittet, D. C. B. (2015). Observations of the icy universe. Annu. Rev. Astron. Astrophys. 53, 541–581. doi:10.1146/annurev-astro-082214-122348
Chaabouni, H., Bergeron, H., Baouche, S., Dulieu, F., Matar, E., Congiu, E., et al. (2012). Sticking coefficient of hydrogen and deuterium on silicates under interstellar conditions. Astron. Astrophys. 538, A128. doi:10.1051/0004-6361/201117409
Clements, A. R., Berk, B., Cooke, I. R., and Garrod, R. T. (2018). Kinetic monte carlo simulations of water ice porosity: extrapolations of deposition parameters from the laboratory to interstellar space. Phys. Chem. Chem. Phys. 20, 5553–5568. doi:10.1039/C7CP05966C
Cuppen, H. M., Ioppolo, S., Walsh, C., Lamberts, T., Semenov, D., Garrod, R. T., et al. (2017). Grain surface models and data for astrochemistry. Spa. Sci. Rev. 212, 1–58. doi:10.1007/s11214-016-0319-3
Cuppen, H. M., Morata, O., and Herbst, E. (2006). Monte Carlo simulations of H2 formation on stochastically heated grains. Mon. Noti. Roy. Astronom. Soc. 367, 1757–1765. doi:10.1111/j.1365-2966.2006.10079.x
Dupuy, J. L., Lewis, S. P., and Stancil, P. C. (2016). A comprehensive study of hydrogen adsorbing to amorphous water ice: defining adsorption in. Classi. Mol. Dyn. 831, 54. doi:10.3847/0004-637X/831/1/54
Ferrero, S., Zamirri, L., Ceccarelli, C., Witzel, A., Rimola, A., and Ugliengo, P. (2020). Binding energies of interstellar molecules on crystalline and amorphous models of water ice by ab initio calculations, Astrophys. J. 904, 11. doi:10.3847/1538-4357/abb953
Garrod, R. T. (2008). A new modified-rate approach for gas-grain chemical simulations. Astron. Astrophys. 491, 239–251. doi:10.1051/0004-6361:200810518
Garrod, R. T. (2013b). A three-phase chemical model of hot cores: the formation of glycine. Astrophys. J. 765, 29. doi:10.1088/0004-637X/765/1/60
Garrod, R. T., and Pauly, T. (2011). On the formation of Co2 and other interstellar ices. Astrophys. J. 735, 15. doi:10.1088/0004-637X/735/1/15
Garrod, R. T. (2013a). Three-dimensional, off-lattice monte carlo kinetics simulations of interstellar grain chemistry and ice structure. Astrophys. J. 778, 14. doi:10.1088/0004-637X/778/2/158
Garrod, R. T., Wakelam, V., and Herbst, E. (2007). Non-thermal desorption from interstellar dust grains via exothermic surface reactions, Astron. Astrophys. 467, 1103–1115. doi:10.1051/0004-6361:20066704
Gillespie, D. T. (1976). A general method for numerically simulating the stochastic time evolution of coupled chemical reactions. J. Comput. Phys. 22, 403–434. doi:10.1016/0021-9991(76)90041-3
Gould, R. J., and Salpeter, E. E. (1963). The interstellar abundance of the hydrogen molecule. i. basic processes. Astrophys. J. 138, 393. doi:10.1086/147654
Hasegawa, T. I., Herbst, E., and Leung, C. M. (1992). Models of gas-grain chemistry in dense interstellar clouds with complex organic molecules. Astrophys. J. 82, 167. doi:10.1086/191713
He, J., Acharyya, K., and Vidali, G. (2016a). Binding energy of molecules on water ice: laboratory measurements and modeling. Lab. Measuremen. Model. 825, 89. doi:10.3847/0004-637X/825/2/89
He, J., Acharyya, K., and Vidali, G. (2016b). Sticking of molecules on nonporous amorphous water ice. Astrophys. J. 823, 56. doi:10.3847/0004-637X/823/1/56
He, J., Clements, A. R., Emtiaz, S., Toriello, F., Garrod, R. T., and Vidali, G. (2019). The effective surface area of amorphous solid water measured by the infrared absorption of carbon monoxide, Astrophys. J. 878, 94. doi:10.3847/15384357/ab1f6a
Hollenbach, D., and Salpeter, E. E. (1971). Surface recombination of hydrogen molecules. Surf. Recomb. Hyd. Mol. 163, 155. doi:10.1086/150754
Jones, A. P., Ysard, N., Köhler, M., Fanciullo, L., Bocchio, M., Micelotta, E., et al. (2014). The cycling of carbon into and out of dust. Farad. Discuss. 168, 313. doi:10.1039/C3FD00128H
Keane, J. V., Tielens, A. G., Boogert, A. C., Schutte, W. A., and Whittet, D. C. B. (2001). Ice absorption features in the 5–8 μm region toward embedded protostars. Astron. Astrophys. 376, 254–270. doi:10.1051/0004-6361:20010936
Kouchi, A., Yamamoto, T., Kozasa, T., Kuroda, T., and Greenberg, J. M. (1994). Conditions for condensation and preservation of amorphous ice and crystallinity of astrophysical ices. Astron. Astrophys. 290, 1009–1018.
Mannel, T., Bentley, M. S., Boakes, P. D., Jeszenszky, H., Ehrenfreund, P., Engrand, C., et al. (2019). Dust of comet 67p/churyumov-gerasimenko collected by rosetta/midas: classification and extension to the nanometer scale. Astron. Astrophys. 630, 14. doi:10.1051/0004-6361/201834851
Mathis, J. S. (1996). Dust models with tight abundance constraints. Astrophys. J. 472, 643. doi:10.1086/178094
Mathis, J. S., Rumpl, W., and Nordsieck, K. H. (1977). The size distribution of interstellar grains. Astrophys. J. 217, 425–433. doi:10.1086/155591
Minissale, M., Dulieu, F., Cazaux, S., and Hocuk, S. (2016). Dust as interstellar catalyst. I. Quanti. Chem. desorp. proce. 585, A24. doi:10.1051/0004-6361/201525981
Navarro-Ruiz, J., Sodupe, M., Ugliengo, P., and Rimola, A. (2014). Interstellar h adsorption and h2 formation on the crystalline (010) forsterite surface: a b3lyp-d2* periodic study. Phys. Chem. Chem. Phys. 16, 17447–17457. doi:10.1039/C4CP00819G
Oba, Y., Miyauchi, N., Hidaka, H., Chigai, T., Watanabe, N., and Kouchi, A. (2009). Formation of compact amorphous H2O ice by codeposition of hydrogen atoms with oxygen molecules on grain surfaces, Astrophys. J. 701, 464–470. doi:10.1088/0004-637X/701/1/464
Ormel, C. W., Paszun, D., Dominik, C., and Tielens, A. G. (2009). Dust coagulation and fragmentation in molecular clouds. i. how collisions between dust aggregates alter the dust size distribution. Astron. Astrophys. 502, 845–869. doi:10.1051/0004-6361/200811158
Ossenkopf, V. (1993). Dust coagulation in dense molecular clouds : the formation of fluffy aggregates. Astron. Astrophys. 280, 617–646.
Pantaleone, S., Enrique-Romero, J., Ceccarelli, C., Ugliengo, P., Balucani, N., and Rimola, A. (2020). Chemical desorption versus energy dissipation: insights from ab initio. Mol. Dyn. HCO.⋅Form. 897, 56. doi:10.3847/1538-4357/ab8a4b
Pauly, T., and Garrod, R. T. (2018). Modeling co, co2, and h2o ice abundances in the envelopes of young stellar objects in the magellanic clouds. Astrophys. J. 854, 1613. doi:10.3847/1538-4357/aaa96a
Pauly, T., and Garrod, R. T. (2016). The effects of grain size and temperature distributions on the formation of interstellar ice mantles. Astrophys. J. 817, 16. doi:10.3847/0004-637X/817/2/146
Pirronello, V., Liu, C., Shen, L., and Vidali, G. (1997). Laboratory synthesis of molecular hydrogen on surfaces of astrophysical interest. Lab. Synth. Mol. Hydrogen Surf. Astrophys. Int. 475, L69–L72. doi:10.1086/310464
Potapov, A., Jäger, C., and Henning, T. (2020). Ice coverage of dust grains in cold astrophysical environments. Phys. Rev. LettersPhys. Rev. Lett. 124, 221103. doi:10.1103/PhysRevLett.124.221103
Raut, U., Famá, M., Teolis, B. D., and Baragiola, R. A. (2007). Characterization of porosity in vapor-deposited amorphous solid water by methane adsorption. J. Chem. Phys. 127, 204713. doi:10.1063/1.2796166
Sandford, S. A., and Allamandola, L. J. (1990). The volume- and surface-binding energies of ice systems containing CO, CO2 and H2O. Icarus 87, 188–192. doi:10.1016/0019-1035(90)90028-8
Senevirathne, B., Andersson, S., Dulieu, F., and Nyman, G. (2017). Hydrogen atom mobility, kinetic isotope effects and tunneling on interstellar ices (Ih and ASW). Mol. Astrophys. 6, 59–69. doi:10.1016/j.molap.2017.01.005
Wakelam, V., Bron, E., Cazaux, S., Dulieu, F., Gry, C., Guillard, P., et al. (2017a). H2 formation on interstellar dust grains: the viewpoints of theory, experiments, models and observations. Mole. Astrophys. 9, 1–36. doi:10.1016/j.molap.2017.11.001
Wakelam, V., Loison, J. C., Mereau, R., and Ruaud, M. (2017b). Binding energies: new values and impact on the efficiency of chemical desorption. Mol. Astrophys. 6, 22–35. doi:10.1016/j.molap.2017.01.002
Weingartner, J. C., and Draine, B. T. (2001). Dust grain-size distributions and extinction in the milky way, large magellanic cloud, and small magellanic cloud. Astrophys. J. 548, 296–309. doi:10.1086/318651
Keywords: astrochemistry, chemical kinetic simulation, dark interstellar clouds, dust porosity, interstellar dust, interstellar ice, molecular hydrogen, surface chemistry
Citation: Christianson DA and Garrod RT (2021) Chemical Kinetics Simulations of Ice Chemistry on Porous Versus Non-Porous Dust Grains. Front. Astron. Space Sci. 8:643297. doi: 10.3389/fspas.2021.643297
Received: 17 December 2020; Accepted: 02 February 2021;
Published: 15 April 2021.
Edited by:
Ryan C. Fortenberry, University of Mississippi, United StatesReviewed by:
Gianfranco Vidali, Syracuse University, United StatesAlbert Rimola, Autonomous University of Barcelona, Spain
Copyright © 2021 Christianson and Garrod. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Drew A. Christianson, dac5tt@virginia.edu