Dynamics of CO in Amorphous Water Ice Environments
L.J. Karssemeijer, S. Ioppolo, M.C. van Hemert, A. van der Avoird, M.A. Allodi, G.A. Blake, H.M. Cuppen
aa r X i v : . [ a s t r o - ph . GA ] J a n Draft version July 25, 2018
Preprint typeset using L A TEX style emulateapj v. 04/17/13
DYNAMICS OF CO IN AMORPHOUS WATER ICE ENVIRONMENTS
L.J. Karssemeijer , S. Ioppolo , M.C. van Hemert , A. van der Avoird , M.A. Allodi , G.A. Blake and H.M.Cuppen Draft version July 25, 2018
ABSTRACTThe long-timescale behavior of adsorbed carbon monoxide on the surface of amorphous water ice isstudied under dense cloud conditions by means of off-lattice, on-the-fly, kinetic Monte Carlo simula-tions. It is found that the CO mobility is strongly influenced by the morphology of the ice substrate.Nanopores on the surface provide strong binding sites which can effectively immobilize the adsorbatesat low coverage. As the coverage increases, these strong binding sites are gradually occupied leav-ing a number of admolecules with the ability to diffuse over the surface. Binding energies, and theenergy barrier for diffusion are extracted for various coverages. Additionally, the mobility of CO isdetermined from isothermal desorption experiments. Reasonable agreement on the diffusivity of COis found with the simulations. Analysis of the 2152 cm − , polar CO band supports the computationalfindings that the pores in the water ice provide the strongest binding sites and dominate diffusion atlow temperatures. Keywords:
Astrochemistry, Diffusion, Methods: laboratory: molecular, Methods: numerical, Molecu-lar processes, ISM: clouds INTRODUCTIONMolecules are important constituents of the interstel-lar medium (ISM). In molecular clouds, they functionas coolants and can be used to trace physical condi-tions like temperature, hydrogen density, and the life-time of the cloud (Herbst & van Dishoeck 2009). Fur-thermore, simple interstellar molecules are necessary pre-cursors to more complex biomolecules which are es-sential for the formation of life (Charnley et al. 1992;van Dishoeck & Blake 1998).In the ISM, a rich gas phase chemistry leads pre-dominantly to unsaturated (organic) molecules whereassaturated molecules are mostly formed on dust grainsurfaces. The diffusive Langmuir-Hinshelwood mech-anism is the primary formation process on ice man-tles. As such, the diffusion of reactants is of key im-portance (Hasegawa et al. 1992; Garrod et al. 2008) andshould be well understood. In many rate equation mod-els, the rate of a simple addition reaction A + B → C isimplemented asd [C]d t = k A+B ( k D,A + k D,B ) [A] [B] , (1)where the square brackets denote the surface concentra-tions of species. For reactions with a barrier, k A+B isthe rate for crossing this barrier and k D,A specifies therate with which species A scans the grain surface. Eq (1) Theoretical Chemistry, Institute for Molecules and Materi-als, Radboud University Nijmegen, Heyendaalseweg 135, 6525AJ Nijmegen, The Netherlands. Division of Geology and Planetary Science, California In-stitute of Technology, 1200 E. California Blvd., Pasadena, CA91125, USA Gorlaeus Laboratories, Leiden Institute of Chemistry, LeidenUniversity, P.O. Box 9502, 2300 RA Leiden, The Netherlands Division of Chemistry and Chemical Engineering, CaliforniaInstitute of Technology, 1200 E. California Blvd., Pasadena, CA91125, USA [email protected] demonstrates the importance of accurate diffusion rates.Experimentally however, measurements of diffusion ratesare tedious and often prone to errors. This has led to alack of accurate information, leaving diffusion as an un-certain process in many models. Diffusion barriers arenow frequently assumed to be a fixed fraction of the sur-face binding energy. Even though this fraction influencesoutcomes of the models significantly (Vasyunin & Herbst2013), a wide range of values is used, ranging from0.3 to 0.8 (Hasegawa et al. 1992; Ruffle & Herbst 2002;Cuppen et al. 2009; Chang & Herbst 2012). Theoreticalstudies into the dynamical behavior of specific adsorbateson model interstellar ices can thus greatly add to thevalue of rate equation models by specifically calculatinginput parameters like energy barriers for diffusion anddesorption.In previous work, we have studied the mobility of COon hexagonal water ice and found that, for this system,the diffusion barrier is equal to one third of the bindingenergy (Karssemeijer et al. 2012). In the present paperwe will apply the same methodology to determine the dif-fusion barriers of a more astrophysically relevant system:adsorbed CO on amorphous solid water (ASW). Further-more, we will present results from isothermal desorptionexperiments, from which we determined the diffusion co-efficient of CO in amorphous water ice environments.There are several reasons for taking H O CO as amodel system. Firstly, CO is the second most abundantmolecule in the ISM and is a key precursor for morecomplicated species like carbon dioxide (Ioppolo et al.2011b), methanol (Watanabe & Kouchi 2002), andformic acid (Ioppolo et al. 2011a). All these reac-tions occur on ice mantles of which H O is the maincomponent (Gibb et al. 2000) so the dynamics of ad-sorbed CO on water ices forms an integral part of themolecular cloud’s chemistry. Secondly, because COforms in the gas phase and freezes out on the grainmantles only under certain conditions (Pontoppidan
Karssemeijer et al. O CO system has been exten-sively studied both experimentally (Bar-Nun et al. 1985;Devlin 1992; Allouche et al. 1998; Manca et al. 2000;Collings et al. 2003b; Ayotte et al. 2001) and theoret-ically (Al-Halabi et al. 2004a,b; Manca et al. 2001) sothere is ample reference material.Finally, we can compare our results on amorphous sys-tems to our previous study on ice Ih to learn aboutthe importance of the surface structure of the ice sub-strate. This difference between crystalline and amor-phous substrates is interesting because this is not ex-amined in astrochemical models, which assume homo-geneous grain surfaces, nor in laboratory experiments,which only probe average properties. Also given theamount of discussion about the porosity of interstellarices (Bossa et al. 2012; Ayotte et al. 2001), we will em-phasize the effect of surface inhomogeneity on the diffu-sion coefficient and binding energy of adsorbed CO.The first experimental measurements of the diffusioncoefficient of CO in ASW were reported only very re-cently by Mispelaer et al. (2013), as part of a larger setof studied species. In our experiments, we use a similartechnique but we have focused only on CO and studieda broader temperature window to get a more accuratevalue for the diffusion barrier. Theoretically, no simula-tions reaching beyond the molecular dynamics timescale(roughly nanoseconds) have yet been performed on amor-phous substrates. We will present these simulations us-ing an off-lattice kinetic Monte Carlo (KMC) approachwhich can probe long timescales but which also has a highlevel of detail because it gives access to the positions ofall individual atoms throughout the simulation. We willshow that the surface structure, and especially its poros-ity, plays a critical role in the CO mobility making theamorphous system essentially different from crystallinesystems we studied before.We will start in Section 2 with a description of theKMC simulations and their results. In Section 3, theexperiments are presented and these are compared withthe outcome of the simulations in Section 4. Astrophys-ical implications, as well as the consequences for largerscale astrochemical models are discussed Section 5. Theconclusions are summarized in Section 6. SIMULATIONSThe dynamics of the water-CO systems are simulatedwith the Adaptive Kinetic Monte Carlo (AKMC) tech-nique (Henkelman & J´onsson 2001; Pedersen & J´onsson2010), which combines the atomistic level of detail frommolecular dynamics with the ability of probing longtimescales from KMC simulations. The technique is im-plemented in the EON code and is described in moredetail in Karssemeijer et al. (2012), which also demon-strates its applicability to molecular systems. In thissection we will only give a short summary of the AKMCmethod before proceeding to the simulations themselves,and their results. http://theory.cm.utexas.edu/eon/ Computational Details
The time evolution of a system in AKMC is describedin exactly the same way as in normal KMC simula-tions (Bortz et al. 1975; Gillespie 1976; Charnley 1998;Chang et al. 2005). A sequence of steps between discretestates with corresponding time increments is generatedbased on computer-generated, pseudorandom numbers.This procedure requires every state to have a unique ta-ble of events (TOE). This table contains all possible pro-cesses and their rates, which can take the system outof its current state, into the next. In traditional KMC,states are defined by a set of occupation numbers on alattice and the TOEs have to be specified before the startof the simulation. In AKMC however, this is not the case.States are defined as local minima on a potential energysurface (PES) which in turn, relies on atomic coordinatesthrough any kind of force field or higher level method.Hence, the particles are not confined to lattice positionsand atomistic details of the system are available (the po-sitions of all individual atoms are known). New statesare discovered on-the-fly, as the simulation progresses,by performing swarms of transition-state searches onthe PES which fill the TOEs by calculating transitionrates. In this work, the transition states are first-order saddle points (SPs) on the PES. The minimum-mode following method (Henkelman & J´onsson 1999;Malek & Mousseau 2000; Olsen et al. 2004) is used tolocate the SPs and rates are estimated using harmonictransition state theory (HTST) (Vineyard 1957).As explained in our previous paper (Karssemeijer et al.2012), an advantage of the KMC method is that, oncethe TOEs of a system are known for a specific tem-perature, one can easily adjust them for simulations atdifferent temperatures, without having to perform newtransition-state searches. We therefore typically performAKMC simulations of a system at a high temperaturefirst, where states are easily discovered. We then use theresulting TOEs for KMC simulations at lower tempera-tures. From the simulations, we obtain the trajectories ofadsorbed CO molecules on ASW substrates. The meansquared displacements of these trajectories are then usedto extract the diffusion coefficient of the adsorbates.Interactions in the system are described by means ofclassical pair potentials for both inter- and intramolec-ular forces. Since we have two molecular species inour systems, three interaction potentials are needed:H O H O, CO CO, and H O CO. The details of allthree potentials are given in Appendix A. Because theH O CO and CO CO potentials were fitted directlyto ab-initio calculations, corrections were made to ac-count for the zero point energy contribution for all bind-ing energy calculations in this work. This procedure isoutlined in Appendix B.In this study, three different amorphous ice substrateswere studied. Each of these has a unique surface mor-phology, due to their different initial structures. By usingthree different substrates, instead of one, we get a betterhandle on the spread in the results due to a particu-lar surface morphology. The substrates were preparedin the following manner using molecular dynamics sim-ulations. An initial sample is generated containing 480water molecules with a density of 0.94 g cm − , the ex-perimental density of low-temperature vapor-deposited YNAMICS OF CO IN AMORPHOUS WATER ICE ENVIRONMENTS O ice (Jenniskens & Blake 1994), at random (non-overlapping) positions with orientations such that theentire system has no net dipole moment. Periodic bound-ary conditions are applied along all three Cartesian axes.This system is equilibrated at 300 K using a Nos´e-Hooverthermostat for 100 ps after which the thermostat is in-stantaneously set to 10 K and the system is left for an-other 20 ps. The thermostat is then turned off and thesystem is left to equilibrate for another 100 ps after whichthe bottom 120 molecules are constrained to their in-stantaneous positions in this bulk amorphous structure.Next, the periodic boundary condition in the z -directionis removed to create a surface and the temperature isincreased to 100 K over a period of 100 ps. In the x -and y -directions, the periodic boundary conditions re-main applied in order to mimic an infinitely large sur-face. After this, the system is once again equilibrated for100 ps, then cooled back to 10 K in 100 ps, and left toequilibrate for another 100 ps. All heating, cooling, andequilibration periods in this procedure were chosen suf-ficiently long to stabilize the energy fluctuations in thesystem to their expected values. Finally the system isrelaxed to the nearest minimum on the PES. The threedifferent amorphous ice substrates generated in this waywill be referred to as S , S , and S .Even though the substrates each contain 360 watermolecules which are free to move, we will constrain theircoordinates in some of the simulations. In this case, anadditional superscript c will be added to the substratename. 2.2. Computational Results
The discussion of the simulation results starts with aninvestigation of characteristics of the amorphous ice sub-strates and the nature of the CO binding sites on theirsurfaces. We then discuss the dynamics of a single ad-sorbed CO molecule on each surface and evaluate theeffect of the CO surface coverage. These results will bedirectly relevant for surface chemistry in the ISM andcan be compared to the diffusion measurements. Fi-nally, we turn our attention to systems with multipleadsorbed CO molecules, corresponding to the late stagein cloud evolution, where CO-dominated ice layers startto form ( ¨Oberg et al. 2011). These systems allow forcomparison to TPD experiments reported in literature.2.2.1.
Substrate Morphology
The morphology of the ASW substrates is expected todetermine the behavior of any adsorbate. We thereforestart with an analysis of the bulk density and porosity ofthe substrates. The bulk density of the ice was obtainedby averaging the density in a set of spheres of radius4 ˚A, centered on a regular grid with a lattice spacingof about 2 ˚A. With increasing z -coordinate, the densityin the spheres drops because they approach the rough,nanoporous surface. If we leave these spheres out of thecalculations, to avoid effects of the surface, we find a bulkdensity value of 1 . ± .
03 g cm − . There is little varia-tion between the three different substrates. The densityis somewhat higher than the experimental value for va-por deposited ice of 0.94 g cm − (Jenniskens & Blake1994). This is probably explained by the presence ofbulk macropores in the laboratory ices, which are absent in our samples.The porosity of the surface itself is illustrated in Fig-ure 1. This figure shows the height of the surface for theamorphous sample S , as well as for a hexagonal ice sam-ple (containing 360 water molecules). The surface heightwas calculated from a regular grid in the x, y -plane witha spacing of about 2 ˚A. Each point on this grid was thenassigned a z -value, the surface height, corresponding tothe point 2.5 ˚A (roughly the H O CO distance) awayfrom the center of mass of the nearest substrate molecule.The amorphous sample clearly has a much larger surfaceroughness than the crystalline sample. Its z -coordinatecovers an almost 10 ˚A range, which is entirely explainedby the nanopores in the surface. A typical cross sectionof such a pore, containing a CO molecule, is shown inFigure 2. The snapshot was taken from the AKMC sim-ulation of CO on substrate S . The position of this crosssection is also indicated in Figure 1. The deeper well inthe upper left corner on the crystalline substrate is not aspecial site but is an artifact of the mapping method. Itcorresponds to one of the hexagons in the ice Ih crystalwhich has a grid point almost exactly above its center,making it appear deeper than the other hexagons. Allthe hexagons in ice Ih are too small for CO to diffuseinto. From the surface height map, the surface area ofthe substrates can also be calculated by performing apolygonal triangulation. Again, little variation betweenthe amorphous substrates was found. The average areais 807 ± , which is 1 . ± .
01 times the area ofthe base of the simulation box. For ice Ih, this factor is1 . ± . AKMC Simulations and Binding Sites
We performed AKMC simulations at T = 50 K onall three substrates, starting from a configuration witha single CO admolecule relaxed on the substrate, at arandom position. The temperature of 50 K was chosenbecause, when compared to lower temperatures, it re-duces the number of KMC steps necessary to explore allbinding sites on the surface. The TOEs are then readilyadjusted for simulations at lower temperatures. Whilethe simulation explores the PES and the system evolvesin time, not only does the admolecule diffuse, but thesubstrate itself also evolves. In contrast with our pre-vious calculations on crystalline ice (Karssemeijer et al.2012), we observe here that the water molecules in theamorphous substrate are rather mobile. Many states arefound where the CO molecule remains in roughly thesame position but where the water molecules have moved Karssemeijer et al.
Figure 1.
Surface height map of ice Ih sample 1 from Karssemeijer et al. (2012) (left) and amorphous substrate S (right). The dashedcontour on the amorphous substrate shows the cross section from Figure 2. The z -coordinate is measured from the lowest atomic coordinatein the system. Periodic boundary conditions are applied in the x - and y -directions. Figure 2.
Typical cross section of an amorphous ice substratewith a surface nanopore containing a CO molecule. The figurecorresponds to an actual state found from AKMC simulations onsubstrate S . The cross section is also indicated in Figure 1, wherethe surface pore is clearly visible (although it has been shifted forclarity). enough for the states to be considered distinct. This is amanifestation of the much rougher and more complicatedPES, which has many more shallow minima compared tocrystalline systems. The hydrogen bond network, how-ever, remains mostly unchanged during the simulation. This roughness has an important consequence for thesimulations. On crystalline substrates, the number ofstates entered by the simulation remained limited sincethe water molecules did not move significantly and theCO molecule could only occupy a finite number of bind-ing sites. For the amorphous systems this is not thecase; since the water molecules also move, the simula-tion keeps entering new, unexplored states and expensivesaddle point searches are almost continuously required.We therefore had to stop the AKMC sampling manually,once confident that a sufficiently large set of states wasexplored.AKMC simulations were also performed on the con-strained substrates ( S c , S c , and S c ), where the numberof states is limited due to the frozen substrate (now onlythe CO molecule is allowed to move). Here, we foundabout 80 states per substrate (see Figure 3). To be surethat we sampled long enough on the unconstrained sam-ples, we visually verified that the binding sites found onthe constrained substrates are also present on the corre-sponding unconstrained sample.For all states entered by the AKMC simulations, wecalculated the binding energy, E B , of the CO moleculeto the ice surface by taking the difference between theenergy of the system with the adsorbed CO and the en-ergy of the substrate, after relaxing it without the CO.The distribution of binding energies and the number ofstates on both the free and the constrained substratesare shown in the left panels of Figure 3. The bindingenergies are found to be broadly distributed between 60and 250 meV, much broader than the distribution forcrystalline ice (between 100 and 150 meV) we found inour previous study.The origin of the broad distribution of binding energiesis revealed by correlating, for each state, the binding en-ergy with the number of H O molecules neighboring theCO. We defined the number of neighbors as the number
YNAMICS OF CO IN AMORPHOUS WATER ICE ENVIRONMENTS N u m be r o f s t a t e s ( no r m a li z ed ) Binding energy (meV) a) S (0), N S =354S (0), N S =547S (0), N S =291
50 70 90 110 130 150 170 190 210 230 250 c) S c1 (3), N S =84S c2 (3), N S =96S c3 (3), N S =75
50 70 90 110 130 150 170 190 210 230 250 b) S c1 (0), N S =81S c2 (0), N S =83S c3 (0), N S =77
50 70 90 110 130 150 170 190 210 230 250 d) S c1 (6), N S =83S c2 (6), N S =86S c3 (6), N S =78
50 70 90 110 130 150 170 190 210 230 250
Figure 3.
Distribution of CO binding energies on the different amorphous ice substrates studied. Panel a) contains the free substrateswith one admolecule. Panels b), c), and d) show binding energies on constrained substrates with 0, 3, and 6 additional CO molecules,respectively. The classification of the substrates is explained in Section 2.1; N s is the number of states explored by the AKMC simulationson each substrate. of H O molecules which have their center of mass withina radius 4.5 ˚A of the center of mass of the CO molecule.The correlation with the binding energy is shown in Fig-ure 4 and has a P-value of 0.79. Physically, this meansthat the nanopores in the ASW substrate provide thestrongest binding sites, which is consistent with exper-imental observations of trapped CO in ASW (Devlin1992; Collings et al. 2003b). Previous theoretical investi-gations of the adsorption of CO on ASW report a similarcorrelation between the CO binding energy and its num-ber of neighbors, though the reported maximum bindingenergy of 155 meV is significantly lower than our high-est binding energy (Al-Halabi et al. 2004a). A possibleexplanation for this, besides the different interaction po-tential, is that in this work, the binding energy was foundfrom geometry optimizations starting from the final con-figurations of 15 ps molecular dynamics trajectories. Aswe shall see at the end of Section 2.2.3, this may not leavesufficient time for the adsorbate to find its way into oneof the pores, where the binding energy is highest.The analysis above shows that when classifying thebinding sites on the surface a good first criterion iswhether the site is a pore site or a surface site. Thepore sites have a higher binding energy, more neighbors,and, as we shall see in the next section, critically influ-ence the mobility of adsorbed CO. This criterion couldbe used as an improvement in astrochemical models, for N u m be r o f ne i ghbo r s (r c o m < . Å ) Binding energy (meV)P = 0.79
Figure 4.
Correlation between the binding energy and the num-ber of neighbors for a single adsorbed CO molecule on ASW. example by including two CO populations.2.2.3.
Single Admolecule Dynamics
The complicated structure of the substrate and thePES also influences the efficiency of the AKMC simula-tions to simulate long-timescale diffusive behavior. Oneaspect is that the roughness of the PES leads to manystates, separated by only low barriers. It is a known prob-
Karssemeijer et al. lem that the KMC algorithm tends to get stuck in thesestates, crossing and recrossing a low barrier many timesbefore eventually evolving over a higher, physically moreinteresting barrier. As we demonstrated before for ice Ihsubstrates, this effect can be countered by using a coarsegraining algorithm (Pedersen et al. 2012). We employedthis algorithm in all AKMC simulations presented here.A second aspect is that the admolecule resides muchlonger in the regions of the substrate with nanoporesthan in the regions with just surface sites, which becomesclear by analyzing the time spent by the admolecule onthe different parts of the substrate. Especially for sys-tems with an unconstrained water substrate, this givesrise to a problem. In this case, we observe in our simu-lations that the CO molecule is able to diffuse into oneof the pores quickly, but once it is there, many differ-ent states are found by the AKMC method where theCO molecule remains more or less in the same positionbut where all molecules, including the surrounding H Omolecules reorient themselves somewhat, in a concertedmotion. One should be mindful of this effect of thenanopores when looking at the histograms of binding en-ergies on the unconstrained substrates in Figure 3. Eventhough the number of physical pores on each substrateis limited (typically about three per substrate), each ofthem contains a large number of states. As for the KMCsimulations, these states inside a pore often have low bar-riers between them, making it difficult for the simulationto escape from a pore region in a reasonable number ofKMC steps. For this reason we were unable to extracta reliable diffusion coefficient on the unconstrained sub-strates. Only on S (0) were we able to make an estimateof 6 . ± . × − cm s − at T = 50 K.On the constrained substrates this effect from thenanopores is much less prohibiting since the substratemolecules cannot move to contribute to a large numberof shallow pore states. This is the reason why Figure 3shows fewer sites with high binding energies on the con-strained than on the unconstrained substrates. The COis able to enter the pores on the constrained samples, butthe water molecules cannot reorient themselves to ac-commodate for the CO molecule, whereby increasing thebinding energy and the number of states. We performedKMC simulations and extracted the diffusion coefficientson all three constrained substrates in the temperaturerange between 25 and 50 K. These results, as well thevalue for the unconstrained substrate S (0), are shownin Figure 5.The diffusion coefficient at T = 50 K for the con-strained substrate S c (0) is more than one order of mag-nitude higher than the value for the unconstrained sub-strate. Even though the latter is not as reliable dueto the bad statistics of the KMC simulations, this dis-crepancy is consistent with comparisons we performedbetween CO diffusion on constrained and unconstrainedice Ih substrates. Also, Batista & J´onsson (2001) reporta lower diffusion energy barrier for water surface diffu-sion on constrained ice Ih surfaces than on freely movingsubstrates. This is because substrate relaxations lowerthe potential energies of minima on the PES more thanat the saddle points. For our systems, the process bar-rier heights on the unconstrained substrates are on aver-age 10% higher than on their constrained counterparts. −24 −20 −16 −12 −8 −4 D ( c m s − ) S (0)S (3)S (6)S (0)S (3)S (6)S (0)S (3)S (6) S (0)Exp. (this work)Exp.(cid:9)(a)Ice 1h (b) Figure 5.
CO diffusion coefficient on frozen substrates with vary-ing coverage. For comparison, the green triangle shows the dif-fusion coefficient on the unconstrained substrate, S (0). Exper-imental data from the isothermal desorption experiments is alsoshown, as well as experimental results from Mispelaer et al. (2013)(a). The value for ice Ih is that of substrate sample 1 fromKarssemeijer et al. (2012) (b). Based on this rough estimate, the diffusion coefficient at50 K, assuming a barrier height of 100 meV, leads to aten times higher CO diffusivity on the constrained sam-ples. For this reason, the diffusion coefficients on theconstrained substrates given in this paper should be in-terpreted as upper limits for the true diffusivity of thecompletely free system.Figure 5 clearly shows that the diffusion coefficient, D ,of a single CO on the constrained substrates is describedby an Arrhenius expression D ( T ) = D exp (cid:18) − E D k B T (cid:19) . (2)By fitting the equation to our data we extracted theeffective diffusion activation energies, E D , and the pre-exponential factors, D , for all three substrates. Thesevalues are listed in Table 1.For reference purposes, we also show our earlier re-sults on the diffusivity of a single CO molecule on ice Ih.On the crystalline substrate, the diffusion coefficient isat least four orders of magnitude larger than on any ofthe constrained amorphous substrates. Regardless thatthe diffusivities on the constrained amorphous substratesshould be considered as an upper limit, it is clear thatthe single CO mobility on an amorphous ice substrateis very low. Our KMC simulations show that, even at T = 50 K on S c (the substrate with the highest mobility),the admolecule spent 98% of the simulation time trappedin one of the pores. So if, in a molecular cloud, a COmolecule lands on an H O dominated ice mantle, it willbe mobile for only a short time, until it reaches a strongbinding pore site. Averaged over all three unconstrainedsubstrates, the time it takes the CO to reach a pore is7 ns at T = 50 K. We determined this time by averag-ing KMC trajectories which we started from a randomweakly bound surface site ( E B <
100 meV) and stoppedonce a strong binding site was entered ( E B >
150 meV).
YNAMICS OF CO IN AMORPHOUS WATER ICE ENVIRONMENTS T = 50 K, a bindingenergy of 150 meV and a diffusion barrier of 100 meVleads to a timescale difference of five orders of magni-tude between diffusion and desorption. This differenceonly increases for lower temperatures. Consequently, theabsence of desorption as a process in our TOEs is not aproblem for our mobility analysis.2.2.4. Filling the Pores
The simulations of a single adsorbed CO on amorphousice have revealed that the porous nature of the substrateeffectively immobilizes the admolecule. However, bothin molecular clouds and in laboratory experiments, thereare multiple CO molecules. Given the previous section,some of these molecules will fill the nanopores which nat-urally raises questions regarding the mobility of the re-maining CO molecules, which are not trapped.We address this question by occupying the three, orsix, strongest binding sites on each substrate with COmolecules. This corresponds to filling the nanopores withadsorbates. These configurations were then relaxed andthis entire system is constrained from movement. Thenwe added one more CO admolecule, free to move, whichwe studied by means of AKMC. These systems with threeor six additional CO molecules are denoted respectivelyby S ci (3) and S ci (6) where i = 1 , , Higher Coverages
In dense cloud conditions, CO typically freezes out ongrains which are already covered with an H O dominatedice ( ¨Oberg et al. 2011). It is therefore interesting to con-sider layered ices which already have considerable CO
Table 1
Diffusion coefficients at T = 50 K, and fitted Arrheniusparameters for diffusion on all CO-amorphous icesystems studied. Experimental data and valuesextracted from previous calculations on hexagonal iceare also listed.System D (50 K) D E D (cm s − ) (cm s − ) (meV)S (0) · · · · · · · · · S c (0) 3 . × − . × − c (3) 2 . × − . × − c (6) 1 . × − . × − (0) · · · · · · · · · S c (0) 2 . × − . × − c (3) 1 . × − . × − c (6) 1 . × − . × − (0) 6 . × − · · · · · · S c (0) 4 . × − . × − c (3) 9 . × − . × − c (6) 5 . × − . × − . × − . × − ± a · · · · · · ± b . × − . × − a Experimental data from Mispelaer et al. (2013) b Theoretical value from Karssemeijer et al. (2012) buildup. Also, from an experimental point of view, sim-ulating the dynamics of just a few admolecules is hardlyrepresentative. In TPD experiments, the surface cover-age is typically much higher (ranging from 0.1 to manymonolayers) and the measured molecules are those whichare not trapped in the pores.For these reasons we have studied the adsorption en-ergy of CO on amorphous ice substrates when they arealready partially covered with CO. Starting from the barewater substrate S , we generate a grid at a distance of3 ˚A above the surface with a lattice spacing of roughly2 ˚A. A CO molecule is then placed at a randomly cho-sen grid point with a random orientation, the geometryis relaxed, and the binding energy of the CO is regis-tered. Next, a new grid is generated and a second COmolecule is deposited in the same way. The energy dif-ference between the relaxed configurations with one andtwo CO molecules is registered as the binding energy ofthe second CO molecule. In this way, 100 CO moleculeswere deposited and the whole procedure was repeated450 times. The average binding energy and its standarddeviation are shown in Figure 6. The surface coverage isalso shown in monolayers, defined as 10 molec cm − ,which is the commonly-used definition. Visual inspectionreveals however, that one monolayer for this system cor-responds more accurately to 0 . × molec cm − , or40 CO molecules on the surface. Because we have aboutthree nanopores on each substrate, about 10% of CO canbe trapped in pores at monolayer coverage.Because our AKMC simulations show no significantdiffusion of CO at low temperature and coverage, thismethod may well represent the experimental situationwhere a CO molecule is deposited from the gas phase ata random position and immediately sticks in the posi-tion where it landed on the ice. This is of course underthe assumption that the CO molecules carry no signifi-cant kinetic or internal energy which would allow them to Karssemeijer et al. diffuse over the surface until they thermalize. In princi-ple the method could be extended by including a certainperiod of equilibration, either by AKMC or moleculardynamics, after deposition.The calculated binding energies show a large variation.This is mainly because the CO molecules are depositedat random positions, and they therefore probe a set ofbinding sites which, as we know from Section 2.2.2, hasa broad distribution of binding energies. Furthermore,the minimization after deposition of the n ’th moleculemay also trigger restructuring of the previous n − n ’th molecule andcould lead to a systematic overestimation of the bindingenergy but we could not verify this. Despite the widedistribution, the mean binding energy follows a smoothtrend and shows the decrease in binding energy with thesurface coverage. This is because the probability of anew molecule finding a strong binding site on the sub-strate decreases with increasing coverage as the numberof the relatively weak CO CO interactions grows whilethe stronger H O CO interactions decrease in number.It is seen that the average binding energy drops fromabout 125 meV at zero coverage to about 75 meV whenthere is already one monolayer of CO present. This is ingood agreement with TPD experiments by Collings et al.(2003a) who report an activation energy of 100 meV forthe desorption of CO, at monolayer coverage, directlyfrom a non-porous ASW surface. From more recent,sub-monolayer, desorption experiments by Noble et al.(2012) the coverage dependence of the binding energycould be extracted from the TPD data. This dependenceis also well reproduced by our calculations. Data fromboth experiments are also shown in Figure 6.
20 40 60 80 100 120 140 160 180 0 10 20 30 40 50 60 70 80 90 1000.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 B i nd i ng ene r g y ( m e V ) Adsorbed CO molecules on the surfaceCoverage (10 molec cm -2 ) Figure 6.
Binding energy of CO on ASW as a function of sur-face coverage. The solid (blue) line shows the coverage dependentbinding energy fitted to sub-monolayer TPD results by Noble et al.(2012), the dashed (black) shows the value Collings et al. (2003a)extracted from TPD experiments. One monolayer is defined as10 molec cm − . EXPERIMENTS Experimentally, we determine the diffusion coefficientof CO in an amorphous water ice environment from therate with which CO desorbs, after it has traveled througha layer of ASW. This is achieved by depositing a slab ofmixed H O:CO ice onto a substrate, followed by a layerof pure H O. This system is then kept at constant tem-perature while the amount of CO in the ice is monitoredusing infrared (IR) spectrometry. The desorption rate ofCO from the ice depends on its diffusivity in the H Ooverlayer. The experimental procedure closely resemblesthe method used by Mispelaer et al. (2013). In our ex-periments however, we study a broader temperature win-dow (32 to 50 K) and use a H O:CO mixture covered byan ASW layer, instead of a pure CO ice to avoid thediffusion of CO upon deposition of the water overlayer.3.1.
Experimental Details
The experiments were performed with the Caltech as-trochemical ice spectroscopy setup which is described indetail by Allodi et al. (2013). It consists of a high vac-uum chamber (base pressure < − Torr) containing asilicon substrate which can be cooled down to 8 K usinga closed-cycle helium cryostat. Gas mixtures can be pre-pared in a separate metal deposition line to be depositedonto the substrate. The IR spectra of the samples arerecorded by means of a Fourier Transform-IR spectrom-eter in transmission mode at a resolution of 0.5 cm − .3.1.1. Experimental Procedure
The isothermal desorption experiments were per-formed at 32, 37, 40, and 50 K starting from an ice whichwas prepared as follows. First, a deposition of a mixtureof H O:CO is carried out at T = 10 K through a 1 / ′′ di-ameter stainless steel pipe that faces the substrate. Theend of the pipe is positioned about 1 ′′ away from the sub-strate and is capped with a metal mesh with a 38 micronhole size to ensure a uniform ice deposition. To removethe most loosely bound CO, this ice is then annealedtwice to 32 K at a rate of 5 K min − , cooling back down,with no additional waiting time, to 8 K in between. Thiscreates a CO-rich ASW film with an H O:CO mixing ra-tio of around 2:1. Finally, an additional H O layer isdeposited at 8 K to form a porous overlayer. The ice isthen heated to the desired temperature at 10 K min − ,and the IR spectra are recorded while the CO diffusesthrough the ice and desorbs from the surface. The mo-ment when the ice reaches the desired temperature istaken as t = 0 in the analysis of the results.The column densities of H O and CO in the ice aremonitored by integrating the characteristic bands of themolecules and dividing by the band strengths of thepeaks. For water we use the 1660 cm − bending modeand for CO the 2139 cm − stretching mode. The corre-sponding band strengths are 1 . × − and 1 . × − cm molec − , respectively (Gerakines et al. 1995). Thearea of the water band was determined by numerical in-tegration while the area of the CO band was calculatedby fitting two Gaussians to the spectrum to distinguishbetween CO in water-rich (polar) and water-poor (apo-lar) environments. The polar component has a distinctpeak around 2152 cm − while the larger, non-polar, com-ponent peaks around 2139 cm − (Sandford et al. 1988;Bouwman et al. 2007). YNAMICS OF CO IN AMORPHOUS WATER ICE ENVIRONMENTS .
94 g cm − forH O (Jenniskens & Blake 1994) and 0 .
81 g cm − forCO (Loeffler et al. 2005). Before the start of the isother-mal experiments, the ice consists of a 1 . ± . µ m thickmixed layer with an 0 . ± . µ m layer of pure H Oon top. This ratio was similar for every temperature andis in good agreement with the deposited amounts of gas(3 Torr of mixture followed by 1 Torr of H O) as mea-sured in the dosing line by a mass independent activecapacitance transmitter. We are therefore confident thatthe sample we start from is the same for all experiments.Furthermore, by performing the annealing cycles beforedepositing the H O overlayer, the most weakly boundCO molecules desorb. This amounts to about 4% of thetotal CO. The annealing procedure limits diffusion of COinto the H O cap during its deposition and makes the icesat the start of each isothermal experiment more similar.After the isothermal experiments, a TPD experimentis performed with a heating of 1 K min − , until the icehas fully desorbed. IR spectra are taken every minuteduring the heating.3.1.2. Determination of the Diffusion Coefficient
The concentration profile of CO, n ( z, t ), in the ice isdescribed using a solution to Fick’s second law of diffu-sion in one dimension: ∂n ( z, t ) ∂t = D ( T ) ∂ n ( z, t ) ∂z . (3)This approach was shown to give good results byMispelaer et al. (2013). The solutions of this equationdepend on the initial conditions. For our purpose we im-pose that n ( h, t ) = 0 to reflect immediate desorptionof molecules which reach the surface (at z = h ) and ∂n (0 , t ) /∂z = 0 because no CO can escape from the bot-tom of the film. Furthermore, the concentration profileat t = 0 is chosen to be either n s or n h : n h ( z,
0) = n , if 0 < z < h, (4) n s ( z,
0) = n n , if 0 < z ≤ d, , if d < z < h. (5)The first function, n h , is used when CO is homoge-neously distributed in the ice by the time it reaches thedesired temperature. The second expression, n s , de-scribes the situation where CO is initially confined toa slab of height d at the bottom of the ice film. For theseconstraints the solutions to Fick’s second law read n h ( z, t ) = ∞ X i =0 n ( − i µ i h cos ( µ i z ) exp (cid:0) − µ i Dt (cid:1) , (6) n s ( z, t ) = ∞ X i =0 n µ i h sin ( µ i d ) cos ( µ i z ) exp (cid:0) − µ i Dt (cid:1) , (7)for n h and n s respectively. Here µ i = (2 i + 1) π/ h and D is the diffusion coefficient. From these expressions,the column density of CO molecules in the ice is readilyfound by integrating over z from 0 to h . These can beconverted to band areas, A h and A s , which can then be fitted to the experimental data. The final expressions are A h ( t ) = s + ∞ X i =0 A − s ) µ i h exp (cid:0) − µ i Dt (cid:1) , (8) A s ( t ) = s + ∞ X i =0 A − s )( − i µ i hd sin ( µ i d ) exp (cid:0) − µ i Dt (cid:1) , (9)where A is the initial band area and s is an offset whichwe have to include to reproduce the experimental data.Physically this corresponds to CO which is completelytrapped in the water ice and cannot diffuse out.More complicated models, with separate diffusion co-efficients for CO in the upper and lower layers, were alsotried. With these models, better fits to the data could beobtained, but we did not find two clearly distinguishablediffusion coefficients. We therefore attribute the betterfit to the increased number of parameters in the modeland decided to stick with the simpler models of Eqs (8)and (9). 3.2. Experimental Results and Analysis
In the isothermal desorption experiments, the IR spec-tra are recorded once the deposited ice reaches the de-sired temperature. The spectrum of the CO peak,recorded after 60 minutes at T = 37 K, is shown in Fig-ure 7. At this time, CO is diffusing through the porousASW overlayer and desorbing from the surface. This isschematically shown in the inset of the figure. The twofitted Gaussians, which were used to estimate the COband area, are also shown. The time evolution of theband areas corresponding to the polar and apolar peaks,are shown in Figure 8. The polar component is alwaysjust a small fraction (6 ± t = 0 value. Due to the small contribution fromthe polar band, the total CO band area is almost identi-cal to the apolar component and it decreases in time, dueto desorption from the ice. The decay times are seen todecrease with increasing temperature which we attributeto faster diffusion of CO through the ASW overlayer. At T = 32 K, an increase is seen in the band area of the ap-olar peak during the first ∼
30 minutes of the experiment.We believe that this is due to changes in band strengtharising from structural changes in the ice. These includethe dilution due to diffusion of CO into the upper layerand possibly the local crystallization of CO. At 32 K wethink this process is slow enough, in combination withthe low CO desorption rate, to be observed in the IRspectra. At higher temperatures, the changes will alsooccur, but these then proceed too rapidly to be observed.The behavior of the polar component, correspondingto CO interactions with dangling OH bonds, differs fromthat of the apolar component. As seen from the lowerpanel of Figure 8, the polar component only decreasessignificantly at 50 K, the highest temperature. At thelower temperatures, the polar component remains almostconstant, while the total amount of CO decreases. At T = 32 K, the polar component even increases with time.This suggests that the ‘polar CO population’ is less mo-bile than the apolar CO and thus corresponds to COoccupying strong binding sites within the ASW. The in-crease at 32 K is then naturally explained because, as0 Karssemeijer et al.
Wavenumber (cm -1 )Wavelength ( µ m) Apolar COPolar COIR spectrum 21002110212021302140215021602170 4.764.744.724.704.684.664.644.62 H O COCOH O:CO
Figure 7.
Infrared spectrum of the CO peak taken at t = 1 hr.in the isothermal desorption experiment at T = 37 K. The insetschematically shows the structure of the ice. diffusion into the upper layer progresses, more of theseenergetically favorable sites become available for the CO,leading to an increase of the polar band. From the spec-tra taken during the TPD, following the 32 K isothermalexperiment, it seems that the polar component decreasesmost rapidly around 40 K, consistent with the data fromFigure 8. Unfortunately, the data is too noisy to drawdefinitive conclusions and the analysis of the IR spectraduring the TPD remains speculative.From the analysis above, we conclude that the apolarCO band is the best measure of the mobile moleculesand we thus use the data from this component to fit ourFickian model. The band area contribution from the po-lar peak is omitted from the fit because it is a measureof molecules which are generally bound too strongly todiffuse. The rapid initial decrease in band area at tem-peratures of 37 K and higher suggest that CO has alreadydiffused into the upper H O layer at t = 0. This meansthat CO in homogeneously distributed in the ice at thestart of the isothermal measurements, so we use Eq (8) tofit the data in these cases. At T = 32 K there is an incu-bation period of about one hour before the CO band areastarts decreasing. This behavior indicates that the CO isstill confined to the lower layer at t = 0 and we thereforeuse Eq (9) as a model for this temperature. Because themodel cannot describe the increasing band area duringthe first 30 minutes at this temperature, these data wereexcluded from the fit.To fit the data, the diffusion coefficient D , the offset s ,and the initial band area A were used as fitting parame-ters. The fit of the initial band area was needed becausethere is too much noise in the data to keep it fixed atthe measured value at t = 0. The thickness of the COcontaining slab d and of the total ice film h were fixedparameters determined from the spectra taken right af-ter the second annealing of the ice and after deposition ofthe H O overlayer, respectively. All parameters, includ-ing the diffusion coefficients, are listed in Table 2 and theArrhenius behavior of the diffusion is shown in Figure 5.From the latter we extract a diffusion energy barrier of26 ±
15 meV.From Figure 8, we see that the model is able to de-scribe the experimental data with reasonable accuracy.Especially the incubation time at T = 32 K is well de- B and a r ea ( c m − )
32 K 37 K 40 K 50 K0.40.60.81.01.2 0 5 10 15 20 25 30 35 40Time (hours)Polar CO component N o r m a li z ed band a r ea Figure 8.
Top panel: band area of the apolar CO component(points) and model fits (lines) as a function of time as measuredduring the isothermal desorption experiments. Bottom panel: nor-malized area of the polar CO band. scribed by the solution with the CO initially confined tothe lower part of the ice. Despite the relatively good fitthere are several effects which lead to large uncertaintiesin the extracted diffusion coefficient. First, some uncer-tainty arises from the determination of the thickness ofthe ice film. Our estimate depends on the densities of COand H O which will change due to heating and, especiallyfor water, depend strongly on the deposition methodand temperature (Stevenson et al. 1999; Dohn´alek et al.2003). The thickness estimate also depends on the bandstrengths, which are also influenced by temperature andby the mixing ratio (Bouwman et al. 2007). We believethat the thickness is the most uncertain parameter inthe model and the error bars in Figure 5 are based on a50% uncertainty of the ice thickness. This however, leadsto a systematic error which affects the pre-exponentialfactor, D , but not energy barrier E D . A related as-pect is the structural change in the water ice duringthe experiment which leads to compaction of the filmdue collapse of macropores (Bossa et al. 2012) and sub-sequently to trapping or release of CO (Bar-Nun et al.1985; Collings et al. 2003b). In a pure ice, this transitionis observed between 38 and 68 K (Jenniskens & Blake1994). Another source of error is the manner in whichthe band areas are extracted from the IR spectra. As canbe seen from Figure 7, the fit is not perfect. We attributethis mainly due to the large amount of H O in the cham-
YNAMICS OF CO IN AMORPHOUS WATER ICE ENVIRONMENTS Table 2
Parameters used to model the isothermal desorptionexperiments.
T D (Fitted) A (Fitted) s (Fitted) h d (K) (cm s − ) (cm − ) (cm − ) ( µ m) ( µ m)32 7 . × − . × − . × − . × − ber, which affects the baseline of the spectrum in the COstretch region. Differences between numerical integrationof the bands versus the fitting of Gaussians and variousmethods of baseline subtraction influence the extracteddiffusion energy barrier significantly. Given these con-siderations, we estimate the uncertainty in the diffusionbarrier to be 15 meV.The experimental results from Mispelaer et al. (2013)are also shown in Figure 5. Even though the authors alsomention several sources of errors in the data, it is reassur-ing to see that there is good agreement between the twoexperiments where the temperatures overlap. We find asomewhat steeper slope on the diffusion coefficient whichis reflected in a higher energy barrier of (26 ±
15) meVagainst (10 ±
15) meV from Mispelaer et al. (2013). Thekey difference between the two experiments is that westart by depositing a mixture of H O and CO insteadof a pure layer of CO. This procedure binds the COmolecules more strongly in the ice film and allowed usto perform isothermal experiments up to 50 K. Addi-tionally, the mixture of H O:CO provides a better ther-mal conduction with the H O overlayer which decreasesthe temperature gradient in the ice. The larger tem-perature range studied facilitates the extraction of thetemperature dependence of the diffusion coefficient andthe calculation of the energy barrier for diffusion. COMBINING SIMULATIONS ANDEXPERIMENTSThe computational and experimental predictions of thediffusivity of CO in an ASW environment can be com-pared in both Table 1 and Figure 5. Although the simu-lations describe a pure surface processes and the experi-ments measure bulk CO diffusion through a ASW layer,the two are comparable given the porous nature of thevapor deposited water ice, if we assume that the size ofthe cracks and macropores in the ASW are sufficientlylarge to interpret the CO diffusion as an effective sur-face process, along the pore walls. This assumption isreasonable because we start from a macroporous ice andthe diffusion rate of CO at these low temperatures ismuch faster than the reorganization rate in ASW, withwhich the porous structure can collapse (Mispelaer et al.2013). The diffusion of CO along the macropore walls isthen comparable to the diffusion in the AKMC simula-tions because the walls of the macropores will also con-tain nanopores, similar to those on the simulated ices.On a qualitative level, we see that there is good agree-ment between the simulations and experiments on twokey points. First and foremost, the diffusion coefficientsare in good agreement on an order of magnitude level.Secondly, both simulations and experiment show that a fraction of CO is immobilized by the water due to trap-ping in nanopores. We will discuss these points in moredetail below.The diffusion coefficients can be compared betweentheory and experiment but one should keep in mind thatin the experiments the coverage of CO molecules is sig-nificantly higher so the mobility is more strongly influ-enced by the weak CO CO interactions than in thesimulations. The simulations where the nanopores arefilled with CO are thus the most representative of theexperimental situation. This is also reflected by the flat-tening of the slope in Figure 5 when the CO coverageis increased. The diffusion barriers extracted from thesimulations with 3 or 6 of the pore sites filled vary be-tween the substrates but are all within 66 ±
20 meV. Thisis higher than the experimental value of 26 meV. Thestronger contribution from the CO CO interactions isone possible explanation for this discrepancy. Anothereffect which could play a role is, as mentioned before,the simultaneous compaction of the ASW film and theCO diffusion through it. At higher temperatures, the icebecomes more compact due to the closing of cracks andthe collapse of macropores. This results in a relativelylower mobility at higher temperatures and thus a flatterslope and lower diffusion barrier. The agreement betweenthe pre-exponential factors, D , is not as good. From thesimulations, this attempt frequency derives mainly fromthe vibrational excitations of the system and the valueswe find are close to the value of 10 − s − which is oftenused (Kellogg 1994). From the analysis of the experimen-tal results we have seen that the pre-exponential varieslargely, though systematically, with the ice film thicknessand conclusive values cannot be obtained. As input forastrochemical models, the pre-exponential factor has lesssignificance than the energy barrier for diffusion.The key importance of the morphology of the watersubstrate is seen in both the experiments and the simu-lations. The simulations show clearly that pores in theASW substrate, even if only sub-nanometer in size, havea critical influence on the mobility of adsorbed CO. Forthe CO to be mobile on the surface, it is a prerequisiteto have either no pore sites or to have them all filled. Inthe experiments, we see the polar CO band, correspond-ing to CO interacting with OH dangling bonds, remain-ing constant, or even increasing in strength, while thetotal amount of CO is decreasing. This indicates thatthis band corresponds to those CO sites with the highestbinding energy. This leads us to the conclusion that theCO signal from the polar band in the IR spectra corre-sponds to CO occupying sites similar to the, nano-sized,pore sites in the simulations.A final remark regarding the comparison between sim-ulations and experiments concerns the effect of the poten-tial energy functions used in the simulations. The some-what higher binding energies of CO on ASW compared tothe experimental value (see Figure 6) could point to anoverbinding in our H O CO potential. Although thediffusion coefficients are not derived from the absolutevalue of the potential but from the energy difference be-tween the minima and saddle points on the PES, it couldbe that the minima are affected more strongly than thesaddle points leading to a higher diffusion energy barrier. ASTROPHYSICAL SIGNIFICANCE2
Karssemeijer et al.
The AKMC simulations, as well as the isothermal des-orption experiments have demonstrated that nanoporesin ASW play an important role in the kinetics of ad-sorbed CO. These pores, where an adsorbed CO moleculecan interact with a large number of water molecules, arefound to have very high binding energies, leading to trap-ping of CO at low coverages. In the Monte Carlo sim-ulations, we observed that a single CO admolecule wastrapped in one of the nanopores for 98% of the simula-tion time, even at T = 50 K. This effectively immobilizedthe admolecule. Diffusion became more rapid once thepore sites were effectively removed from the ice by fill-ing them with CO molecules. In the experiments, thestrong binding energy of the nanopore sites was deducedfrom the increased intensity of the 2152 cm − , polar CO,band, associated with CO interactions with dangling OHbonds in the ice.The large influence from the morphology of the sub-strate on the CO mobility is also seen when comparing toour previous results on crystalline ice (Karssemeijer et al.2012). Simulations at the lowest coverage, with just asingle CO admolecule show that the diffusion coefficienton crystalline substrates is about four orders of mag-nitude larger than on amorphous surfaces. The energybarrier for diffusion, E D , on crystalline ice (49 meV) is50 % lower than in the amorphous case ( ∼
100 meV).The diffusion prefactor is however rather unaffected bythe morphology of the substrate.Given the amorphous character of interstellar dustgrain mantles, the presence of small pores will be a keyfactor determining formation rates on grain surfaces. Inthe early stages of dense cloud formation, before catas-trophic CO freezeout, nanopores are likely to trap carbonmonoxide molecules for very long times. This will affectthe surface chemistry because the pores can act as reac-tive sites in this case. Although the overall mobility ofreactants will be low, they will tend to get trapped inthe same places, giving more time for reaction to occur,which is especially favorable if there is a reaction barrierto overcome. This was found by Fuchs et al. (2009) forhydrogenation reactions in CO ices.The diffusion barrier of hydrogen atoms on ASW ismuch lower than that of the CO molecules (Perets et al.2005; Matar et al. 2008; Hama et al. 2012). Recent ex-periments by Hama et al. (2012) have shown that themajority of H atom binding sites on ASW are shal-low, with diffusion barriers ≤
22 meV. A small frac-tion of the sites was found to have higher diffusion bar-riers ( ≥
30 meV), which might correspond to nanoporesites. Based on the results of Minissale et al. (2013), alsothe diffusion barrier of atomic oxygen on ASW is lowerthan that of CO. This rapid diffusion of atomic hydro-gen and oxygen will lead to an efficient conversion of CO,trapped in pore sites at low coverages, to CH OH, CO and HCOOH (with OH as a potential intermediate). Ifthis conversion is efficient enough, there will be no COleft in the pore sites of the H O dominated ices. Thisis consistent with the non-detection of the 2152 cm − inastronomical spectra (Pontoppidan et al. 2005) and thesuggestion that CO is mixed with CH OH in dust grainmantles to account for the red component of the COband (Cuppen et al. 2011).Once more CO freezes out on the grain, the nano-sizedpores get filled and the remaining CO can diffuse much faster. At the same time however, the remaining COmolecules will also desorb more easily, because they aremore weakly bound. We computationally studied thisdecrease in binding energy as an increasing amount ofCO molecules was adsorbed on an ASW substrate andfound good agreement with experimentally determinedtrends. This good agreement also adds to the reliabilityof the simulations at low coverages, where experimentaldata are still scarce.To provide astrochemical modelers with necessary in-put parameters, we determined the energy barrier for dif-fusion from both the simulations and the isothermal des-orption experiments. The simulations were found to bein reasonable agreement with the experiments as well aswith similar experiments carrier out by Mispelaer et al.(2013). Even though there are many uncertainties, theseare to the best of our knowledge the only available dataon CO diffusion in ASW. The analysis has shown thatthe CO mobility and its binding energy are highly depen-dent on both the position on the ice surface as well ason the CO coverage. In this respect, amorphous surfacesare essentially different from crystalline substrates, whichshow much less inhomogeneity. Because diffusion is animportant parameter, modelers should try to include asmuch of these local variation as possible and avoid takingdiffusion barriers as fixed ratios of the binding energy.There are several approaches to include effects from in-homogeneity in the ice mantle in astrochemical models.In lattice gas KMC, this can be done by making the bind-ing energy and diffusion barrier height site-dependent toaccount for pore sites. This site-dependent approach wasused by Cuppen & Herbst (2005) to account for surfaceroughness in simulations of H formation on dust grainanalogs. In rate equation and master equation methods,one could include two populations of CO molecules, sim-ilar to the approach taken by Cuppen & Garrod (2011)for H . One population would represent the immobileCO molecules in the nanopores and the other the moremobile molecules on the surface. Based on the resultspresented in this paper we suggest surface binding ener-gies of 130 and 80 meV for the strong and weakly boundpopulations, respectively, and diffusion energy barriersof 80 and 30 meV. Another possibility to account for in-homogeneity is to include a direct dependence on surfacecoverage of the diffusion barriers and binding energies.In lattice KMC models this would be in the spirit of thework by Fuchs et al. (2009), where the sticking probabil-ity of impinging H atoms was given a dependence on theH surface coverage. CONCLUSIONSWe have studied the dynamics of CO in amorphouswater ice environments at low temperatures by means ofkinetic Monte Carlo simulations and isothermal desorp-tion experiments. The main conclusions of this analysisare the following:1. The CO mobility is highly dependent on the mor-phology of the ice. At the lowest coverage, thepresence of nanometer-size pores in ASW increasesthe energy barrier for diffusion to around 100 meV;twice the value of 50 meV for crystalline ice, whichdoes not contain pores.2. The surface coverage of CO on ASW critically in-
YNAMICS OF CO IN AMORPHOUS WATER ICE ENVIRONMENTS A. INTERACTION POTENTIALS
This appendix describes the three interaction poten-tials used in this work. All molecules are fully flex-ible and their internal motions are described by in-tramolecular potentials. Intermolecular interactions con-tain an electrostatic contribution from a set of pointcharges on each molecule and several functions de-scribing the van der Waals contributions. For theH O H O interactions, the TIP4P/2005f potentialwas used (Gonz´alez & Abascal 2011). This modelwas developed as a flexible version of the successfulTIP4P/2005 potential which gives a good description ofthe condensed phases of water (Abascal & Vega 2005;Vega & Abascal 2011). The H O CO potential, aswell as the CO intramolecular potential are described inKarssemeijer et al. (2012). The intermolecular CO CO
Table 3
Potential parameters for the CO COBuckingham potential.Interaction A ij B ij C ij (eV) (˚A − ) (eV˚A )C-C 361.4 2.835 33.45C-O 1517 3.543 15.19O-O 6370 4.252 10.55 potential is described below. Because it has not beenpublished before, we have included the details of the fit-ting procedure. For the sake of completeness, we alsogive the intramolecular part of the potential again.A.1. CO-CO Potential
Following up on the work by Vissers et al. (2003), thepotential energy surface of the CO dimer was calculatedfrom a set of interaction energies on a grid of geometries.When the C O bond length is fixed, four coordinatesdescribe the geometry of the system: R , θ A , θ B , and φ . The distance R is the length of the vector R fromthe center of mass of monomer A to that of monomerB. The angles θ are between R and the vectors r CO(A) and r CO(B) , which point from the C to the O atom inthe respective monomer. The dihedral angle φ is be-tween the planes spanned by ( R , r A ) and ( R , r B ). The4D grid consists of 7 θ A angles, 7 θ B angles, 6 φ angles,and 13 R values (3.5, 3.75, 4.0, 4.25, 4.5, 5.0, 5.5, 6.0,6.5, 7.0, 9.0, 15.0 and 20.0 ˚A). The angles were chosenin order to enable a spherical expansion of the interac-tion energy (see below). Calculations were performed forthree C O distances, i.e. molecule A was kept at theground state equilibrium distance r e = 1 .
128 ˚A whilemolecule B has r e , 1 . r e , and 0 . r e . The interaction en-ergy was found from CCSD(T) calculations using a stan-dard aug-cc-pVQZ (Woon & Dunning 1993) basis setwith the Boys-Bernardi counterpoise correction. All cal-culations were performed with the Molpro (Werner et al.2008) program.A spherical expansion of the potential (van Hemert1983) was made to analyze the potential energy surface.This expansion was used to generate contour plots ofthe potential. When both molecules are in one plane,there are two minima separated by a barrier. The low-est minimum, with an interaction energy of -16.7 meV,occurs with the two CO molecules parallel with θ A and θ B angles of 135 ◦ and -135 ◦ respectively, so with the twocarbon atoms closest together. The other minimum is at − . θ A and θ B at 60 ◦ and -60 ◦ respec-tively, with now the oxygen atoms closest together. Thecenter of mass distances for the minima are 4.5 and 3.7 ˚Arespectively.For the AKMC simulations, the spherical expansionparametrization is too expensive. Instead, the ab-initio interaction potential was parametrized as a site-site po-tential with electrostatic, exchange repulsion, dispersion,and intramolecular contributions: V CO CO = V el-st + V exch + V disp + V intra(A) + V intra(B) (A1)The electrostatic part, V el-st , contains interactions be-tween charges located on each atom and on the molecular4 Karssemeijer et al. centers of mass. The values of the charges are initiallychosen to exactly reproduce the dipole and quadrupolemoments at the specific r CO . The moments were takenfrom MCSCF/CCI calculations with the aug-cc-pVQZbasis set. The exchange repulsion and dispersion termsare expressed as a Buckingham potential between theatomic sites in the dimer: V exch + V disp = X i ∈ A X j ∈ B A ij exp ( − B ij r ij ) − C ij r ij . (A2)The parameters A ij and B ij were optimized in a leastsquares procedure where all interaction energies with avalue below 25 meV were included. The limit of 25 meVwas chosen in order to focus the optimization of the po-tential on the bound part. The values of the B ij pa-rameters were initially derived from the relation betweenstandard Lennard-Jones C and C coefficients and theparameters used here, as given by Lim (2009). The B ij parameters and the charges were then slightly adapted inorder to minimize the least squares standard deviation.All parameters for the Buckingham potential are listed inTable 3. It was found that the interaction energy couldbe well be represented only when the electrostatic termwas made dependent on the intramolecular distance r CO .More specifically, the changes in the charges were madeproportional to the changes in computed charge as de-rived from the ab-initio dipole and quadrupole moments: q i = q i exp ( − σ i ( r CO − r e )) . (A3)The charges q i are − . e on the C atom and − . e on the O atom. The values of σ i are 3 .
844 and 2 .
132 ˚A − for C and O, respectively. The parameters of the Buck-ingham potential remain independent of r CO . The in-tramolecular interactions within each monomer are de-scribed with a Morse potential: V intra = D e [1 − exp ( − γ ( r CO − r e ))] , (A4)where D e = 11 .
23 eV and r e = 1 .
128 ˚A are the experi-mental dissociation energy and equilibrium bond length.These agree to within 0.1% with the ab-initio calcula-tions. The γ parameter was fitted to reproduce the ab-initio bond length dependence of the potential energyand has a value of 2 .
328 ˚A − .A.1.1. Quality Of The Potential
Contour plots were also generated from the site-sitepotential. Given the simple form of the parametrized po-tential, the agreement between the full ab-initio contoursand the model contours is satisfactory. In the in-planecase there are again two minimum energy structures sep-arated by a barrier, the energy ordering of the minimais however reversed and the center of mass distances aresomewhat different. The model molecules are softer thanthe ab-initio molecules. This is in part due to the bias inthe selection of configurations used in the least squaresfit. Nevertheless, the use of these model parameters inmolecular dynamics (MD) calculations of crystalline COdoes not lead to large structural deviations from experi-ments, as becomes clear from the analysis below.The quality of the potential was tested by a seriesof MD simulations. These were based on the analyt-ical forces derived from the site-site potential, the ve- locity Verlet integration scheme, and a Berendsen ther-mostat (Allen & Tildesley 1989). As a first test, pureamorphous CO crystals consisting of 200, 300, 500, 800,and 1200 CO molecules were created. These crystals weregrown by adding CO molecules, one by one, to the previ-ous molecules. Each new molecule was positioned at 10 ˚Afrom the surface of the core formed by the molecules al-ready present. The four angles describing the initial ori-entation and center of mass position of the new moleculewith respect to the ones already present were determinedby a random number generator. The energetically mostfavorable final position of the new molecule was then de-termined with the simplex method. In this procedurethe positions of the molecules already present were keptfrozen. The energies in the simplex procedure were de-rived from the site-site potential. When all CO moleculeswere added the crystal was made to undergo temperaturecycling in an MD procedure, in 6 K steps from 0 to 30 Kand back to 6 K. At each temperature the crystal waskept for 200000 a.u. of time ( ∼ α -CO (P2 × × g CC (r) C−CAmorphousCrystallineExperimental g OO (r) O−OAmorphousCrystallineExperimental2 3 4 5 6 7 8 9 10 g C O (r) r (Å) C−OAmorphousCrystallineExperimental Figure 9.
Radial distribution functions, g ( r ), in arbitrary units,from 6 K molecular dynamics simulations of amorphous CO andcrystalline α -CO. Data from the experimental crystal structure of α -CO (Vegard 1930) is shown for reference. The structure of the various samples is clearest whenlooking at the radial distribution functions. Becausethere is little difference in the radial distribution func-tions for the various sizes we consider here only two
YNAMICS OF CO IN AMORPHOUS WATER ICE ENVIRONMENTS α -CO crystal structure. The distribu-tions were obtained by averaging over the structures dur-ing the last 100 steps in the 5 ps re-equilibration run at6 K. The broad distributions of the amorphous samplesshow clearly that these systems are still amorphous afterthe temperature cycling. The crystalline systems showmuch sharper peaks and also reflect a slightly higherdensity. The most notable difference between the twois the relatively short nearest neighbor O-O distance inthe amorphous crystals when compared with the crys-talline sample. Due to temperature broadening the splitnearest neighbor peak seen in the experimental C-O dis-tribution cannot be seen in the MD simulations. Energyminimizations of the α -CO structure do show this fea-ture, however the splitting is less pronounced. From thecrystalline sample, we calculated the density at 6 K as1.05 g cm − . This compares well with the density of1.03 g cm − derived from the X-ray based fcc unit celllength of 5.64 ˚A from Krupskii et al. (1973). The densityof the amorphous systems is about 2% lower than thatof the crystalline samples.From the MD data, we have derived the specific heatof the crystal and made an estimate of the binding en-ergy that would be obtained for an infinite crystal. Theaverage specific heat corresponds to 5 × − eV K − .By extrapolating the data to 0 K we estimate the co-hesive energies to be 81 meV for the amorphous systemand 84 meV for the crystalline one. The latter compareswell to the experimental value of 86 meV for crystallineCO (Kelley 1935). B. ZERO POINT ENERGY CORRECTIONS