© ESO, 2016

1. Introduction

Ices are a prominent material that covers many planets, moons and small bodies in the outer solar system and in other planetary systems (Gudipati & Castillo-Rogez 2013). In addition, small ice grains – and ice-covered silicate grains – are ubiquitous in the interplanetary and interstellar medium as dust particles (Draine 2003). The composition of the ices is complex and assumed to be a mixture of various species such as H₂O, NH₃, CO₂, and others (Bottinelli et al. 2010). The fate of these grains is governed by the further accretion of material by collisions with gas-phase molecules or other grains. In addition, they are exposed to a radiation environment that consists of both a photonic and a massive (ion) component (Baragiola et al. 2013). The latter is composed of the stellar winds (in planetary systems), magnetospheric ions (in the vicinity of planets), and cosmic radiation (everywhere). Such radiation has many effects on the grains; it may compact the material, erode it, induce phase transformations (to the amorphous phase or to other solid phases), and induce chemical transformations in the grain (Johnson 1990; Johnson et al. 2013). By ejecting material, such irradiation events may set chemically transformed material – also including highly reactive radicals – free by bringing it to the gas phase.

Since cosmic radiation presents a ubiquitous radiation source, a study of its consequences is most relevant. Due to their larger energy deposition, heavy nuclei are more relevant to study than proton irradiation, even though their fluxes are smaller (Shen et al. 2004). Energies have a broad maximum at 10–100 MeV/amu. Estimations of the flux density show that the flux of so-called heavy nuclei (12 ≤ Z ≤ 29) with energies between 0.1 and 10 MeV/amu is around 0.05 cm^-2 s^-1 and hence dominates over the solar wind flux already beyond the Mars orbit (Pilling et al. 2010b, Fig. 8).

The effects of energetic particle irradiation of ices has recently been studied using controlled laboratory experiments, in which ice samples are irradiated by energetic heavy ions; an overview is given by Allodi et al. (2013). Among the ices studied are water ice and mixtures with carbon dioxide and ammonia (Pilling et al. 2010b; Bordalo et al. 2013), ammonia (Pilling et al. 2010a), methanol (de Barros et al. 2011b), and methane ices (de Barros et al. 2011a). These studies report on extensive chemical transformation found in the irradiated samples.

Impact of MeV ions on matter leads to the formation of a so-called ion track (Fleischer et al. 1975). The physics is as follows: the swift ion delivers its energy to the electronic system of the target. The electrons convey their energy to the atomic system by a complex series of processes which depend in detail strongly on the molecular composition of the material. In consequence, the material is heated along the path of the ion, and thermodynamic processes (shock wave emission), chemical transformation, as well as material ejection can occur. Due to the complexity of the processes taking place in the electronic system, theoretical modeling of the processes occurring in an ion track usually starts with the time when atoms have obtained kinetic energy, and studies the subsequent processes (Johnson 1990; Toulemonde et al. 2006). Studies that incorporate the electrons in the analysis are comparatively rare and have only been performed for relatively simple materials such as LiF (Cherednikov et al. 2013).

Computer simulation on the basis of molecular dynamics (MD) offers an important theoretical tool to study the processes occurring in ion tracks. Previous simulations focused on ion tracks in plane surfaces. They started with studies of the sputtering of large organic molecules from tracks by Fenyö et al. (1990). Urbassek et al. (1994) reported on the influence of track radius and energy on the sputtering of a model Lennard-Jones solid. Later, Beuve et al. (2003) studied the influence of the spatial and temporal structure of the deposited-energy distribution on the sputtering process, again in a model solid; this study was later continued by Mookerjee et al. (2008).

Studies of ion tracks in grains are comparatively rare. Bringa & Johnson (2002) used MD and later Bringa & Johnson (2004) a thermal model to study sputtering from ion tracks in ice grains. This model was later expanded to calculate the amorphization of grains (Bringa et al. 2007)

Processes occurring in condensed-gas targets under ion impact were studied experimentally by David & Michl (1989) and theoretically by Urbassek & Michl (1987) and Balaji et al. (1995). Later simulation work concentrated on diatomic and triatomic species (Anders & Urbassek 2007, 2009, 2010). While much work concentrated on mild energy depositions in the ice (Johnson & Liu 1996), the consequences of a track in a molecular target, N₂, could also be studied (Kafemann & Urbassek 1993).

Recently, the so-called REAX potentials were developed (van Duin et al. 2001; Chenoweth et al. 2008; Liu et al. 2011) which model the breaking and re-formation of chemical bonds. These potentials are based on the bond-order idea that models the environment dependence of chemical bonds. The ionic contribution to bonding is taken into account by dynamically calculating the partial charges of each atom. These potentials thus allow one to study ion impact on molecular targets and to monitor the dissociations and molecular reactions occurring. These potentials were used by Anders & Urbassek (2013) to study collisions between ice clusters and single-ion impact into an ice target.

In the present paper, we use MD simulations based on the REAX potentials to investigate the processes occurring in an ion track in an ice grain. We use a specific mixture presented by Martins et al. (2013) and focus attention on an ion track that is characterized by an energy deposition of 1 keV/nm and a track radius of 5 Å. In addition, an ion track with doubled radius is studied. We report on shock wave emission from the track, ejection of atoms, molecules and clusters and on the chemical transformations in the target and in the ejecta.

2. Method

2.1. Ice grain and interaction potential

The ice grain consists of a mixture of H₂O, CO₂, NH₃, and CH₃OH with molar ratios of 9.10:8.00:3.75:1.00. This mixture was proposed by Martins et al. (2013) in their study of shock wave interaction with ices; a similar mixture has been used by us previously to study ice grain collisions (Anders & Urbassek 2013). Several physico-chemical parameters of the species are assembled in Table 1. We prepare a spherical grain with a radius of 89 Å by using PACKMOL (Martínez et al. 2009) to put a total of 225, 864 atoms – 28, 428 H₂O, 24, 992 CO₂, 11, 715 NH₃, and 3, 124 CH₃OH molecules – into the sphere and to generate an amorphous structure. The grain was relaxed for a time of 50 ps in order to remove excess kinetic energy and reduce the pressure to zero. After relaxation, the grain has a radius of 93 Å, an average density of 0.070 Å^-3 and a density in the center of 0.074 Å^-3. Its temperature is 1 K.

We use the REAX potential developed by Van Duin et al. (van Duin et al. 2001; Chenoweth et al. 2008; Liu et al. 2011) to model the intra- and intermolecular interactions in this ice mixture. Note that this potential also allows us to model molecular fragmentation and reactions. We use the potential parameters in its latest form (Monti et al. 2012), since this potential was optimized for the amino acid glycine, and this amino acid contains all relevant atom types for the present study. Towards higher energies we add the Ziegler-Biersack-Littmark (ZBL) potential (Ziegler et al. 1985) to model the repulsive short-range interaction. We have described the details of the combination of the REAX and ZBL potentials previously (Anders & Urbassek 2013).

2.2. Simulation procedure

The simulations are performed using the massively parallel open-source code LAMMPS (Plimpton 1995). It uses a variable time step which varies from initially 0.01 to finally 0.1 fs. The simulation runs until 49.3 ps after ion-track formation.

We form an initial “ion track” in the grain as follows. The ion is assumed to pass along a sphere diameter L = 182.3 Å through the center of the grain. In a cylindrical region with radius r = 5 Å around the ion path, we deposit an energy of dE/dx = 1 keV/nm, i.e., in total E_dep = L·dE/ dx = 18.23 keV is deposited in the grain. We note that such an (electronic) energy deposition corresponds, for example, to the stopping of a 3.5 MeV Fe ion in the ice mixture according to the SRIM code (Ziegler 2004), but is by no means the maximum deposition which is achieved for faster projectiles. The track radius is small compared to the grain radius. The energy is given to each atom in the geometrical ion-track region as kinetic energy; each of the N_track = 1101 track atoms receives on average a kinetic energy of E₀ = E_dep/N_track = 16.56 eV, but with random direction. In detail, the code samples the individual track atom velocities from a three-dimensional Maxwellian velocity distribution with temperature E₀/ 3k = 64.04 kK, where a specific heat of 3k per atom is assumed and k is Boltzmann’s constant, The subroutine “velocity track create” within the LAMMPS code ensures that the total momentum and the total angular momentum of the grain vanish. Figure 1 gives an impression of the grain containing the heated track shortly after track initiation.

We can calculate the equilibrium temperature, T_equ, that the grain will reach at long times after track formation via $\mathrm{3}\mathit{Nk}{\mathit{T}}_{\mathrm{equ}}\mathrm{=}\mathit{L}\mathrm{·}\mathrm{d}\mathit{E}\mathit{/}\mathrm{d}\mathit{x}\mathit{.}$(1)With the values given above, and the total number of atoms N in the grain, the equilibration temperature amounts to T_equ = 304 K. We note that this is a gross estimate since the energy used for bond breaking, energy losses due to phase changes (melting and evaporation) and particle ejection, as well as a possible nonequilibrium of the various degrees of freedom present in the molecular system are not taken into account. It shows, however, that we expect considerable heating of the grain by the ion track.

Our simulation volume is a cuboid that dynamically grows to include all atoms. In order to allow for efficient parallelization, all particles (i.e. molecular clusters, isolated molecules and atoms) that reach a certain boundary (here: a sphere of 120 Å radius) are not followed further in the simulation, but their positions and their final velocities, as well as their assignment to molecular clusters, are recorded for further evaluation. Thus, we constrain the simulation volume and allow a reasonable assignment of processors to simulation subvolumes. Atoms that reach the boundary are not followed further in the simulation, but their positions and their final velocities, as well as their assignment to molecular clusters, are recorded for further evaluation.

2.3. Detection algorithms

We use several detector algorithms to identify important quantities of our simulation. They are described in the following.

Local temperatures are evaluated for each atom as the local average of the kinetic energy in a sphere with radius of r_c = 8.51 Å, containing roughly 180 atoms; in order to subtract the directed energy of motion, the kinetic energies are evaluated in the center-of-mass system of the sphere. Similarly local densities are determined as the average around each atom in a sphere with radius of r_c (Colla & Urbassek 1997). We note that the concept of “local temperature” is often used in MD simulation, but may be open to some doubt due to the small atomic ensembles used, and also because of the question whether these ensembles are thermally equilibrated. However, the local temperature may be considered as representing the average local kinetic energy density.

Ejection from the grain is identified by using Stoddard’s (1978) cluster detector. It identifies all atoms that belong to the same cluster by using a distance criterion, 10 Å. This scheme allows for an efficient dynamic grouping of the atoms in the simulation volume into clusters. The largest cluster is identified as the grain; the other groups are either isolated atoms, molecules, or molecular clusters, and all contribute to the ejection yield.

The cluster detector by Stoddard (1978) is also used to identify molecules from its atomic constituents. It uses a distance criterion to identify all atoms belonging to the same molecule. Then a smaller distance d is used as the distance criterion. We use d to denote the maximum distance between two atoms such that they are assigned to the same molecule. In view of the fact that molecules may be vibrationally excited the distance must be chosen larger than the equilibrium bond length. We chose a value of d = 1.7 Å for C, N, and O atoms; while for H we adopted 1.3 Å (Mücksch et al. 2014). This scheme allows for an efficient dynamic grouping of all atoms in the irradiated target into molecules. We have used it previously in several instances (Anders & Urbassek 2007, 2009, 2013; Anders et al. 2013; Mücksch et al. 2014).

3. Results

3.1. Time evolution of density and temperature

Figure 2 visualizes the fate of the grain after track formation in the form of snapshots at selected times. The points where the track penetrates the grain surfaces lead to abundant ejection. While the equatorial surfaces of the track remain cold until the very last snapshot displayed, the radius of the heated zone around the track axis broadens with time as a result of heat conduction. Emission has slowed down in the last snapshot, and the ion track has become considerably depopulated.

Because of the sudden formation of the ion track, a pressure wave passes through the grain. Its characteristics are plotted in Fig. 3, where density and temperature profiles are shown as a function of radial distance r to the track axis. The data apply to a 60 Å high equatorial slice, and have been averaged in azimuthal direction and over the height of the slice.

Within 0.5 ps the density in the center of the track decreases to 25% of the original value. The material has moved sideways and the motion of a compression peak towards the grain surface is seen. The height of the compressive peak decreases with time and its amplitude broadens. When the compression wave reaches the surface – at around 3.2 ps – the surface bulges out. At later times, the outer surface undergoes radial oscillations. In the meantime, the center of the track has become more or less completely emptied of material; but material starts to move in from the sides and may close the track again on time scales large compared to our simulation time.

The temperatures reach exorbitantly high values immediately after track formation. These temperatures express the high energy densities present in the track core. After around 1 ps, the temperature drops below 5000 K. These high temperatures are considerably above the critical temperatures of all species considered in this work, see Table 1; this fact explains the high molecular emission observed in the snapshots of Fig. 2. At later times, the temperature in the track center continuously decreases; but even at the end of the simulation, it has values of around 500 K, sufficient to drive abundant evaporation at least for the species with the smallest critical temperature, CO₂.

Heat is transported out of the central track region both by diffusion and by the compression wave discussed above. Beyond a radius of around 30 Å from the track center, the temperatures are only slowly increasing and are below 300 K even at the end of the simulation. As the final snapshot of Fig. 2 shows, the surface temperatures around the equator remain cold and do not influence emission.

3.2. Ejection

In Fig. 2 we have already shown that the ion track induces abundant ejection. In total, 3573 atoms in 1278 molecules are ejected. The early ejection is caused by the intense energy deposition in the track itself and may be termed (prompt) “sputtering”; it is due to the high pressure building up in the track, which drives material not only sideways, see Fig. 3, but also in axial direction out of the track region. The emission at later times is not driven by pressure, but is due to the high temperatures prevailing in the track center, and also in the surrounding regions, see Fig. 3. It may hence be termed “evaporation”.

We can follow molecule emission by monitoring the number of initial molecules that have been lost from the grain, Fig. 4a. Note that molecules can disappear from the grain also by chemically transforming to other species. We therefore also provide Fig. 4b which gives the number of molecules lost solely due to reactions.

For all species with the exception of H₂O reactions occur quickly, within the first 2 ps after ion track formation. At later times the molecule numbers only slightly change, with a decreasing tendency. Water forms a striking exception by showing a non-monotonous behavior. Beyond 2 ps, the number of water molecules increases again, as a result of recombination processes. We think that this peculiar behavior is caused for two reasons: (i) H⁺ as the dissociation product of water is highly mobile, and can quickly find its recombination partners; (ii) since H₂O is the majority species in the ice mixture, the probability of finding the suitable recombination partners is high.

The total number of molecules lost from the grain, Fig. 4a, follows well the number of reactions monitored in Fig. 4b. For water, however, we see a decrease setting in after around 5 ps, which counteracts molecule restoration by recombination. This decrease is caused by the long-term continuous molecule evaporation. Evaporation is also seen for NH₃, but less for methanol. It is particularly pronounced for CO₂, which features the highest loss rate from the grain.

The total number of original molecules ejected is given in Table 2. These numbers are correlated with the (inverse of the) critical temperatures of these species, see Table 1. CO₂ with the smallest critical temperature of 304 K features the highest ejection yield, followed by NH₃ and CH₃OH in the inverse sequence of their critical temperatures. Water is again an exception; despite the high critical (and also triple-point) temperature, emission is substantial, and occurs over a large period of time, see Fig. 4a. Again, we correlate this behavior with the high recombination propensity of H⁺. Since recombination generates heat, it locally heats the grain enabling the water molecules to be ejected. The correlation between the termination of recombination reactions – at 10 ps, see Fig. 4b – and of water ejection at roughly the same time – see Fig. 4a – corroborates this argument. We note that the enthalpies of mixing in the ice mixture may also play a role in determining which species is most easily emitted. In particular, since water is the dominant species, hydration energies may influence emission, as hydrated species need to break their hydration bonds before emission; otherwise they may be emitted as part of a hydrated cluster, see the discussion in Sect. 3.4 below.

We conclude that grain heating by the track can be decisive for the losses if the final grain temperature is in the range of the critical temperatures.

In summary, Fig. 4 demonstrates that our simulation time is sufficient to model the chemical transformations in the molecular species, and also the prompt sputtering processes, but not necessarily the evaporation processes. Since the temperature in the grain, heated by the ion track, is high compared to the triple point temperatures of the species, Table 1, evaporation will continue on long time scales until evaporation cooling will stop evaporation. The grain temperature at the end of the simulation amounts to 168 K.

3.3. Energy and angle distributions

We summarize important characteristics of the ejection process in Fig. 5. These data are for particles that have reached at least a distance of 120 Å from the grain center. The data are taken at the end of the simulation, such that further changes in the data may be assumed to be small. Such changes may be affected by various processes: further ongoing evaporation; collisions among the ejecta; or break-up of emitted clusters.

Figure 5a displays the distribution f(E) of the kinetic energy of the ejected particles, E. To be precise, E denotes the translational center-of-mass energy of all ejected molecules (or isolated atoms) and f(E)dE equals the (normalized) number of molecules with energy in the interval [ E,E + dE ]. The distribution has a maximum in the range of E = 0.075–0.10 eV. The distribution can be fit by a Maxwellian, $\mathit{f}\mathrm{\left(}\mathit{E}\mathrm{\right)}\mathrm{=}\mathit{CE}{\mathit{e}}^{\mathrm{-}\mathit{E}\mathit{/}\mathit{k}{\mathit{T}}_{\mathrm{fit}}}\mathit{,}$(2)where C is a constant, k is Boltzmann’s constant, and T_fit is a fit temperature. As shown in Fig. 5a, a fit with T_fit = 1093 K represents the data quite satisfactorily in the energy range of E ≤ 0.5 eV; a number of particles are emitted with higher energies and are not represented well by the Maxwellian fit function. Fits to Maxwellian distributions have been attempted already previously when describing data of emission from ion impact on condensed-gas targets (Haring et al. 1984). We note that if we detect the energy distribution at an earlier time, 21 ps, see Fig. 5c, we obtain a slightly hotter distribution, with a maximum in the range of E = 0.125–0.15 eV and a fit temperature of T_fit = 1946 K; note that the Maxwellian fit works for a wider energy range – up to 1.5 eV – than for the later time. The apparent cooling of the distribution results from the late emission of CO₂ molecules with a low temperature whose contribution decreases the average temperature of the entire ejected ensemble.

The temperature in this fit, Eq. (2), should not be assigned an immediate physical relevance; rather it denotes an average temperature. Since emission occurs over a long time, and the temperature in the emission region is constantly decreasing, the emission distribution might be thought of being composed of a series of (weighted) Maxwellians. Still, the relatively high temperature T_fit – consistent with the high energy of the most probable emission events – indicates the importance of early emission events, where the central track temperature is still high.

Analogously, we display in Fig. 5b the angular distribution of emitted molecules. Let ϑ denote the emission angle of the center-of-mass velocity of an ejected particle measured with respect to the track axis; positive angles denote the upper, and negative angles the bottom hemisphere. f(cosϑ)dcosϑ equals the (normalized) number of particles with cosine of the emission angle in the interval [ cosϑ,cos(ϑ + dϑ) ]. We note that an isotropic emission from the surface of the grain corresponds to a constant function, f(cosϑ) = const. Our data, Fig. 5b, show that emission is far from isotropic; it is rather concentrated on emission along the angles ϑ = 0° and 180°; that is, along the track axis. The difference in emission in the upper and lower hemispheres are due to thermal fluctuations in the emission process.

We can fit the angular distribution with an overcosine distribution, $\mathit{f}\mathrm{\left(}\cos \mathit{\vartheta }\mathrm{\right)}\mathrm{=}\mathit{a}\mathrm{+}\mathit{b}\mathrm{|}\cos \mathit{\vartheta }{\mathrm{|}}^{\mathit{\gamma }}\mathit{.}$(3)The constant term, a corresponds to isotropic emission, to evaporation from throughout the grain surface, while the second term models a jet component. We can provide a satisfactory fit with a = 0.12, b = 1.4, and γ = 2.65. The ratio of the integrals under the two contributions to Eq. (3), [ b/ (γ + 1) ] :a = 3.2, shows that the jet component is dominant in the emission distribution. Such peaked distributions have been found previously for emission from ion tracks, both in simulations (Urbassek et al. 1994) and in experiments (Toulemonde et al. 2002), and have been interpreted as being due to an emission jet along the ion track axis, such as found here.

3.4. Cluster emission

The data assembled in Table 2 also compares the ejection yield of isolated molecules with that bound in larger molecular clusters. We see that NH₃ and CH₃OH are often emitted as parts of a larger cluster, while H₂O – and even more so CO₂ – are ejected as single molecules. This partitioning is clearly connected to the importance of the evaporation contribution to the total emission yield. CO₂ and also H₂O show a high fraction of evaporation contributing to the emission, see Fig. 4a; since evaporation occurs via single molecules, it increases the single-molecule contribution to the ejection yield. NH₃ and CH₃OH, on the other side, feature no evaporation, due to their high critical temperatures. In consequence, their single-molecule contribution is suppressed.

This analysis implies that early emission occurs with a large percentage through molecular clusters. Such a feature is well known for sputtering under high-energy density conditions, the so-called spikes, that occur under heavy ion and in particular cluster ion impacts (Colla et al. 2000; Kissel & Urbassek 2001; Urbassek 2008; Sandoval & Urbassek 2015).

Figure 6 displays an overview over the molecular clusters emitted from the grain up to the end of the simulation. All clusters in which at least two molecules are combined are included. In an attempt to identify the molecular content, different colors are used for the different molecular species. Quantitatively, we find that the astonishingly high fraction of 35% of the molecules contained in clusters are fragments and product molecules. Accordingly the fractions of original molecules are reduced to 31% for CO₂, 26% for H₂O, and only 7% for NH₃, and 1.2% for CH₃OH. Thus, CO₂ has become the majority component, evidently due to the intense CO₂ evaporation occurring throughout the simulation. NH₃, and in particular CH₃OH are depleted in the clusters, since they leave the grain preferentially in the form of fragments rather than as intact molecules.

In detail, we find 127 molecular dimers, 22 trimers, and nine quadrumers. The largest clusters found by us, containing five or more molecules, have the composition:

• (H₂O)₄(CO₂)₂(NH₃)₂(OH),

• (H₂O)₆(CO₂)₂(NH₃),

• (H₂O)₄(CO₂)₂(NH₂)(NH₄),

• (H₂O)₄(CO₂)(NH₃)(NH₄),

• (H₂O)₄(CO₂)(NH₃), and

• (H₂O)₄H.

So these largest clusters feature a prominent contribution of H₂O, more than the average values reported above. We argue that this is caused by the large hydration binding energy that H₂O binds to other molecules.

3.5. Chemical transformations

We assemble the species found in the grain at the end of the simulation in Fig. 7a in a form analogous to mass spectra. The original molecular species have been omitted in the plot. The figures demonstrates that a rich variety of fragment and product molecules have been formed, ranging in intensity from one molecule to more than a hundred. The most abundant species are NH₂ (mass 16 amu), OH (17 amu) and CO (28 amu). Similarly, Fig. 7b displays the species found in the ejecta, again by omitting the original molecules whose emission characteristics are assembled in Table 2. Again OH and CO are prominent, while NH₂ is somewhat decreased. Table 3 provides all species which were produced with at least three specimens.

Let us first discuss the fragments and products evidently deriving from water. While OH is abundant, only very few (three) isolated H atoms are found; this is due to the high reactivity of H in the irradiated ice mixture. H₂ is found to a higher percentage, and so are O₂ and H₃O. H₂O₂ is found only rarely, and no O₃ was detected.

In the ejecta, isolated H is relatively frequent and even more abundant than H₂. This is due to the fact that during the early hot sputtering phase, dissociation products such as H (and also OH) are emitted and then find no collision partners to recombine or react with. Similarly, single O atoms are found frequently in the ejecta, but also reacted O₂.

Similarly, NH₃ fragments – in particular NH₂, but also a few NH – are found both in the grain and in the ejecta. However, the product NH₄ is found enriched in the grain, but only rarely in the ejecta. Further frequent reaction products are NO, both in the grain and in the ejecta, and HNO, which is enriched in the ejecta, since it is highly reactive and all but one HNO species in the grain have reacted.

Furthermore, carbon containing molecules are found. Here CO is the major fragment, both in the grain and in the ejecta. But, interestingly, the second in abundance is formaldehyde (H₂CO), which derives from the dissociation of methanol. Also the methanol fragment CH₂OH – and the methoxy group H₃CO, which we can differentiate by inspection of its molecular geometry – are abundantly produced.

We conclude that fragments are found enriched in the ejecta, while product molecules tend to form – and remain – in the irradiated grain.

Since the simulation keeps track of the identity of atoms, we can follow the chemical reactions in detail. We describe a few interesting results below:

• O₂ is formed mostly from CO₂ fragmentation (46 out of 76 O atoms); the next largest source is H₂O (27 O atoms);

• H₂ originates mostly from all the H-containing species, in the order of H₂O (59 out of 96 H atoms), NH₃ (29) and CH₃OH (8);

• the 37 formaldehyde (H₂CO) molecules produced take their C and O atoms mostly from methanol (33 C and 27 O); the remainder of the C atoms stem from CO₂ and the remainder of the O atoms originate from H₂O and CO₂;

• by far the most recombinations leading to the same product molecule occur in H₂O, where 1304 recombination reactions were monitored; the next species is NH₃, where only 36 recombinations occur.

In the following, we compare our simulation results with those of laboratory astrophysics experiments, in which irradiation by cosmic rays is replaced by swift-heavy ion impacts. The formation of fragments and reaction products is observed spectroscopically in the irradiated ice targets, such that we can compare with our results in the grain, but not with the ejecta. While the specific ice mixture we use in our simulations has not been studied, several ices and also ice mixtures have been investigated experimentally.

1. The processes occurring in a frozen methanol target after impact of energetic particles using laboratory astrophysics techniques are studied by de Barros et al. (2011b). They find that the induced CH₃OH dissociation gives rise to the formation of molecular species, particularly H₂CO and CH₄, but also CH₂OH, CO, CO₂, HCO and HCOOCH₃. In their Fig. 7, they propose a pathway for the methanol dissociation and the successive chemical reactions. OH abstraction generates CH₃, and by H⁺ addition: CH₄. On the other hand, consecutive H abstraction generates H₂CO and eventually CO, from which O addition leads to CO₂. The abundant formation of formaldehyde and of CH₂OH parallels our results; however, we find only little methane (two molecules). CO is formed in our simulation probably more as a CO₂ fragment than by multiple-step reactions from methanol.

2. Ammonia ice is studied by Bordalo et al. (2013). They find new species: hydrazine (N₂H₄), diazene (N₂H₂), molecular hydrogen (H₂), and nitrogen (N₂). In water-ammonia ice mixtures they identify additionally nitrous oxide (N₂O), nitrogen oxide (NO), nitrogen dioxide (NO₂), and hydroxylamine (NH₂OH). They find that ammonia destruction is enhanced in the presence of H₂O in the ice. These studies are extended by Pilling et al. (2010a) to mixtures of water ices with ammonia and/or CO. Among the novel species, they find HNCO, N₂O, OCN⁻, and NH$\begin{array}{c}\mathrm{+}\\ \mathrm{4}\end{array}$. While H₂ is amply formed in our simulations – but may also derive from methanol fragments – we observe no N₂ or any other of the molecular products involving two N atoms. Indeed, we checked whether free N atoms were produced at some step in the simulation, but found none; this explains why no N₂ is created. However, the reactions with O are also seen in our simulations, as we observe NO abundantly, and also (rarely) NO₂. Interestingly, also NH₂OH is found in our simulations, albeit more abundantly in the ejecta. NH$\begin{array}{c}\mathrm{+}\\ \mathrm{4}\end{array}$ is a product commonly seen in our simulation.

3. Ice mixtures of H₂O:CO₂ were studied by Pilling et al. (2010b). They find as major products CO, H₂O₂, and to a lesser extent O₃, H₂CO₃ and CO₃. We find CO as the most abundant fragment in our simulations, and also H₂O₂, but only in small amounts. However, no O additions leading to molecules with three O are observed. However, we observe a few molecules containing two C atoms, among which are C₂O, C₂O₂, and (COH)₂.

As a result of this comparison to available experimental studies we conclude the following:

• The molecules found in our study parallel many of those found in experiment. While this appears trivial for the fragments identified, it is reassuring for the true product molecules encountered.

• However, we do not observe N₂ or any molecules containing two N atoms. While we do observe O₂ – mainly as a dissociation product of CO₂ – we do not observe O₃ or any other molecules containing three O atoms. However, we observe molecules containing two C atoms, which must be true reaction products. We surmise that this lack of reaction products is caused by the small density of radicals formed and the short time available to molecules for reacting. In experiments, larger radical densities may build up since for the ion fluences used – up to 2 × 10¹³ cm^-2 in Pilling et al. (2010a) – several ion tracks may overlay each other.

• When searching for possible interesting and astrophysically relevant organic molecules containing both C and O, possibly repeated impacts into a grain should be simulated. We note, however, that comparable laboratory experiments only found traces of such molecules. In a methanol: ammonia ice mixture, Caro et al. (2014) identify a variety of organic species caused by energetic impacts, containing both C and N; as an example, they find urea NH₂CONH₂ (an amide). However, many of the compounds are only formed during warm-up of the ice. This paper discusses the probable reaction pathways leading to these products.

3.6. Changed ion-track radius

While the total energy deposition – or stopping power – of a cosmic ray ion of known energy can be predicted quite accurately, the exact modeling of the radial energy deposition profile is open to some arbitrariness. In order to obtain some information about the effects of the radial energy deposition profile, we performed a second simulation (simulation 2) where we doubled the track radius, 10 Å, but kept an identical energy deposition.

In simulation 2, with increased track radius, the initial temperatures reached in the track are smaller than in the standard simulation, while the final temperatures reached after equilibration of the grain is the same via Eq. (1). Here we observe similar chemical transformations as in our standard simulation: in fact, as Table 3 shows, the identical species are produced as in the standard simulation, with only very few exceptions. This feature corroborates our chemical analysis and shows that the chemical reactions proceeding after ion track formation are reproducible in the simulation. In the species that are produced only rarely – less than three molecules – some fluctuations exist between the two simulations. The total number of chemically transformed species amounts to 1099 molecules in the simulation 2 as compared to 1122 in the standard simulation. Of these, 385 (445) are emitted, while 714 (667) are retained in the grain in simulation 2 (standard simulation). We conclude that for the amount of chemical transformations, it is mainly the total deposited energy that counts and to a lesser extent the maximum temperatures achieved in the track.

However, the total number of original species emitted is increased in simulation 2 for all species, except NH₃, see Table 2. Also the total ejection yield has risen from 3573 atoms in the standard simulation to 4314 atoms. This happens because prompt sputtering, and to a greater extent long-term evaporation, is effective over a larger area, thus increasing the total ejection flux. We note that the gas-flow model of condensed-gas sputtering (Urbassek & Michl 1987), which describes atom emission from a region that has been energized above the critical temperature, predicts only a slight dependence of the yield on the track temperature, while below the phase transition temperature, the Arrhenius law (exp−(U/kT)) introduces a strong dependence.

The overall reduced temperature in simulation 2 leads to a slight reduction in the effective temperature of ejecta, which amounts to only 989 K compared to 1093 K in simulation 1. More pronounced is the reduction in the maximum kinetic energy of the ejecta from 16.3 eV (simulation 1) to 6.7 eV (simulation 2). Also the degree of clustering in the emission process has been increased, and we now find 147 molecular dimers, 23 trimers and 10 quadrumers. In these the proportion of H₂O and CO₂ has risen by 41% and 24%, respectively, as compared to simulation 1. The largest clusters (five or more molecules) found are

• (H₂O)₂₄(CO₂)₅(NH₃)₇(CH₃OH)(OH)₂(NH₂)(CH₃O),

• (H₂O)₁₄(CO₂)₄(NH₃)₂(CH₃OH)(NH₂)₂,

• (H₂O)(CO₂)₂(OH)₂(CO)₂(NH),

• (H₂O)2(CO₂)₂(NH₃)(NO),

• (CO₂)₃(NH₃)(H₃O).

Here in particular the largest clusters contain H₂O over-proportionally.

In summary, a doubled track radius – at identical total energy deposition – reduces the central temperature by a factor of around four. Changes in the chemical transformations and in the properties of ejecta are small, in the order of 10–20%. Most pronounced is an increase in the fraction of clustered ejecta.

4. Summary

We present a MD simulation study of the consequences of the impact of a cosmic ray ion into an astrophysical ice grain. The use of a REAX potential, which allows the breaking and forming of chemical bonds, permits us to study not only the mechanical and thermal effects of the irradiation but also the chemical transformation of the molecules present in the original ice mixture. We focus on the case of an energy deposition of 1 keV/nm into a H₂O:CO₂:NH₃:CH₃OH mixture.

The heated track emits a shock wave into the surrounding grain. While the temperatures in the track are initially exceedingly high, they quickly drop as the shock wave penetrates the grain. Upon reflection from the grain boundaries, the wave excites radial oscillations of the entire grain.

As a consequence of the strong energy deposition, abundant molecule fragmentation is observed. We find CO, OH, and NH₂ as the most abundant fragments.

The fragments are found to undergo subsequent chemical reactions. These are partly recombination reactions, restoring the original molecules. However, new products are also formed. The most abundant product molecules are formaldehyde (H₂CO), NH₄, H₂, H₃O, and O₂.

The ion track abundantly sputters matter from the grain; at later times also strong evaporation adds to the mass loss from the grain. The ejecta are composed of both original and chemically transformed molecules. Besides isolated atoms and molecules, molecular clusters are emitted.

The number of all original molecules continuously decreases after the impact. However, H₂O forms an exception as the number of molecules recovers again after a few ps. This is caused by the particularly strong tendency of OH⁻ and H⁺ to recombine again after dissociation, and is caused by the strongly polar nature of water.

The energy distribution of the ejecta is characterized by a Maxwellian energy distribution with temperature of around 1100 K. Its angular characteristics are governed by a jet-like structure aligned with the ion-track axis.

Fragments are found enriched in the ejecta, while product molecules tend to form – and remain – in the irradiated grain.

For the number of chemical transformations, it is mainly the total deposited energy that counts and less the maximum temperatures achieved in the track. For the emission yields, however, grain heating by the track can be decisive for the losses if the final grain temperature is in the range of the critical temperatures.

Acknowledgments

We acknowledge financial support by the Deutsche Forschungsgemeinschaft within project Ur 32/27-1. Simulations were performed at the High Performance Cluster Elwetritsch (RHRK, TU Kaiserslautern, Germany).

References

All Tables

All Figures