Coulomb Spike Modelling of Ion Sputtering of Amorphous Water Ice

: The effects of electronic excitations on the ion sputtering of water ice are not well understood even though there is a clear dependence of the sputtering yield on the electronic stopping power of high-energy ions. Ion sputtering of amorphous water ice induced by electronic excitations is modelled by using the Coulomb explosion approach. The momentum transfer to ionized target atoms in the Coulomb ﬁeld that is generated by swift ion irradiation is computed. Positively charged ions produced inside tracks are emitted from the surface whenever the kinetic energy gained in the repulsive electrical ﬁeld is higher than the surface binding energy. For that, the energy loss of deep-lying ions to reach the surface is taken into account in the sputtering yield and emitted ion velocity distribution. Monte Carlo simulations are carried out by taking into account the interactions of primary ions and secondary electrons ( δ -rays) with the amorphous water ice medium. A jet-like anisotropic ion emission is found in the perpendicular direction in the angular distribution of the sputtering yield for normal incidence of 1-MeV protons. This directional emission decreases with an increasing incidence angle and vanishes for grazing incidence, in agreement with experimental data on several oxides upon swift ion irradiation. The role of the target material’s properties in this process is discussed.


Introduction
The understanding of electronic excitation processes in solids is still a challenging issue, as regards, for instance, the radiation effects by fission fragments for nuclear applications, such as inert matrix fuels for actinide transmutation [1], or by cosmic rays for space applications, such as space travel [2,3] and satellites [4]. Among these, ion emission from surfaces induced by swift ions is an important issue [5]. Ion sputtering is known to be a very useful tool for various applications, such as thin film-controlled deposition [6] or nano-patterning of surfaces [7]. Besides, a more specific interest was also given to sputtering of water ice for astrophysics [8][9][10]. Two regimes of sputtering were found: the nuclear collision process dominates for low energy ions, below 10 keV, whereas the electronic excitation process dominates for higher energy ions [11]. The former regime is well described by Sigmund's collisional model [8], whereas the latter one is not fully understood. The sputtering yield of water ice increases dramatically with the electronic stopping power in the latter case [12]. It deals with the transfer of electronic excitation to the kinetic energy of surface atoms which is a similar problem to ion track formation in solids. Two classical approaches were developed for tackling such problems: namely, the two-temperature model (TTM) or inelastic thermal spike (i-TS) model [13][14][15] and the Coulomb explosion model [16][17][18][19][20].
In the case of ionic-covalent solids such as crystalline LiF and UO 2 , huge sputtering yields were measured for swift heavy ion irradiations with a high electronic stopping power [21,22]. The broad isotropic angular distribution of the yield was interpreted on the basis of the i-TS model. However, this latter model could not account for the puzzling jetlike emission in the normal direction on top of this broad background. A similar directional sputtering was found for near-IR pulse laser sputtering of crystalline Al 2 O 3 on time scales shorter than 1 ps [23].
Molecular dynamics (MD) simulations on the 150 fs time scale of amorphous solid Ar have also shown that such a directional sputtering could occur perpendicularly to the sample plane [24]. Other MD simulations of ion sputtering have also predicted atomic emission for time scales shorter than the Debye characteristic time (τ D = ω D −1~1 ps) [18]. MD simulations of Coulomb explosion in silicon exposed to highly charged ions have shown that a shock wave is formed on a time scale shorter than 100 fs [18]. This is thought to occur when electrons are pumped out by the impinging ions from the surface by Coulomb interaction. In the case of ice, the H 3 O + ion emission upon 100 keV Ar ion irradiation at 80 K was related to charging effects [25]. There is clear experimental evidence of space charge formation and dielectric breakdown in ice [25][26][27].
Ion sputtering was widely investigated for low and high energy heavy ions in a wide range of solids [5]. Typical cross-over behavior was found for ice from the nuclear collision regime for low energy ions and the electronic excitation regime for high energy ions [11]. A huge increase of the sputtering yield versus stopping power was found for the latter process in water ice [8,11,12] and other ionic-covalent materials such as LiF, CaF 2 , UO 2 , and Y 3 Fe 5 O 12 [5,20,28]. Attempts were made to reproduce the variation of the total yield (Y tot ) integrated over all emission angles versus electronic stopping power (S e ) with the i-TS model for swift heavy ion irradiation [20,21]. The i-TS model was currently used to reproduce the variation of ion track radius versus S e for amorphizable as well as nonamorphizable solids [13]. A power law dependence of Y tot on S e was actually found for swift heavy ions with an exponent ranging between 2 and 4 for insulators [5,28], whereas a quadratic dependence was found for the lower energy ions [8]. The atomic collision process is well described by Sigmund's collisional model, whereas the electronic excitation process is not convincingly modelled.
In the case of LiF, for instance, the binding energy of surface atoms deduced from the i-TS model applied to the dependence of Y tot on S e is much lower than the actual value by a factor of about two [21]. Moreover, the anisotropic jet-like emission for crystalline materials such as LiF and UO 2 could not be modelled quantitatively. Note that this directional emission vanishes at high fluences for LiF in relation to the corrugation of the crystal surface [5]. However, this jet-like emission does not occur for amorphous silica (a-SiO 2 ) [29]. The isotropic component for LiF is found to originate from the near-surface layer, whereas the jet component may be produced by deeper layers [28]. These authors claimed that this directional sputtering derives from a hydrodynamic effect on the pressurized vapor produced at high temperature in the tracks [21,22]. However, no quantitative model was given for this special feature. Therefore, our aim was to investigate the relevance of a Coulomb spike effect on ion sputtering of amorphous water ice on very short time scales, by assuming a sustained electrical field induced by the ionizations inside ion tracks. It is a convenient way to test the likelihood of a Coulomb explosion in a solid to explain the sputtering yield. We show how ion sputtering with the jet-like emission can be explained by pure Coulomb field effects by using the RITRACKS computer code [2,3] devised for the liquid water medium that can be considered as a surrogate for amorphous water ice, as regards the ionization processes [30], electron pulse radiolysis [31], and the dielectric function [30,32]. In addition to ionization, electronic excitation, and dissociative electron attachment, RITRACKS considers vibrational and rotational excitations, which are characteristic of liquid water. Track structure calculations in ice require phonon excitation in place of the molecular excitations but, as discussed later, only ionizations are considered in this study; therefore, the differences of the excitation channels do not impact the applicability of RITRACKS. Hence, we have surmised that the sputtering of amorphous water ice can be modelled with the RITRACKS Quantum Beam Sci. 2023, 7, 7 3 of 16 code. The impact of target material properties in this process is also discussed at length. This study can serve as a basis for addressing the ion sputtering of other dielectric materials such as oxides upon electronic excitations.

Method
We use the RITRACKS computer code [2,3] considering the liquid water medium at 0 • C. This Monte Carlo (MC) code can carry out the track-structure simulations of the primary ions on the basis of ionization cross section and excitation cross section of water. The secondary electrons produced by ionization reactions, the so-called δ-rays, and the valence electrons knocked out by these secondary electrons are also tracked considering reaction channels, namely, ionization, excitation, vibration-rotation, dissociative electron attachment, elastic collisions, and Bremsstrahlung until they reach 2 eV. The valence electrons knocked out by the electrons are also transported. The primary ions are transported until they reach the target surface defined as one of the setup parameters. On the other hand, electrons are tracked until they get thermalized. Unlike systematic fitting formulae [33,34], RITRACKS simulates all these reactions on an event-by-event basis by explicitly tracking every secondary electron and identifying the produced radiolytic species. Thus, RITRACKS can predict the spatial coordinates, timing, and transferred energy for each incident particle allowing the simulation of spatial and statistical fluctuation of the ionization reactions and excitation reactions. It should be noted that the RITRACKS always assumes that the target is liquid water. The cross sections and radiolytic species production are calculated on this premise.
The basic setup of the calculation performed by RITRACKS was a primary proton with kinetic energy of 1 MeV directed to a 500 nm-thick and infinite-wide water slab target. Every time the primary proton knocks out an electron from a water molecule, its spatial coordinate was dumped to an external file as a point of ionization. This feature was used by activating the "events" option of RITRACKS. Since the range of 1-MeV proton is much longer than 500 nm, it is reasonable to assume that the proton energy is constant. Calculation is conducted for 10 4 histories by using different random number seeds to achieve good statistical accuracy better than 10 %, except for insignificant minor events. This statistical accuracy is justified by the fact that the angular distributions of the ejected ions, the main scope of this paper presented later, have fluctuations in the order of 10%, which is small enough to analyze the distribution trends. The spatial coordinates of ionizations directly induced by the incident particle are the output at the end of all the histories and are used to calculate the electric potential between the resulted ions. We have considered that the produced ions are H 3 O + [25]. The kinetic energies of ions which are driven by the potential are calculated. The kinetic energy (E c ) owing to 500 neighboring positively charged ionized atoms is computed based on the spatial coordinates of ionizations. The direction of the ions' motion is calculated on the basis of the electric field exerted by the most neighboring ions. This calculation is justified by the fact that the electric potential is inversely proportional to the distance between charges while the electric force, which determines the direction of the kinetic motion, is inversely proportional to the square of the distance. Calculations for 2000 and 4000 neighboring ionized atoms show that the electric potential attributed to the ionized atoms beyond the 500th neighbor is negligible. Ionized atoms in the ion track are emitted from the surface if E c ≥ E b , where E b is the surface binding energy. A low value of E b = 450 meV was selected, as recommended by the literature for water ice [8]. Since no value is available for H 3 O + ions, the binding energy of water molecules is used. However, the dependence of E b on the direction of emission is not taken into account in the present simulations because the kinetic energy of the emitted H 3 O + ions is in the order of eV whereas E b is in the range of meV.
Contributions of deep layers are included with corrections to E c due to the energy loss (∆E) of H 3 O + ions to escape from the depth to the target surface, as computed with the SRIM2013 MC codes [35]. SRIM2013 cannot calculate the range of H 3 O + ions so the stopping power of 19 O was used instead. The calculations of projected range and range Quantum Beam Sci. 2023, 7, 7 4 of 16 straggling of 19 O + ions are plotted as a function of the kinetic energy below 100 eV ( Figure 1, which are calculated by fitting based on LSS theory [36,37]. The calculation by SRIM of the nuclear stopping power, which is dominant at low energies, is more accurate than that by the LSS theory [38]. We assumed that those ions are one of the main emitted charged species for high-energy projectiles [8,9,25,39]. The possible emission of neutral species is of course not addressed by the present model based on Coulomb interactions. Multiple ionizations in tracks may reduce the yield of H 3 O + ions for heavy ions with a higher stopping power [40] which are not taken into account in the present case for the lower density of ionizations produced by protons in tracks. Benchmark of RITRACKS code is performed against the radial distribution of deposited energy by ionization and electronic excitation for 1-MeV protons (p + ) in water and number of ionizations per primary proton ( Figure 2a) with reference to literature data [2,3]. The primary δ-ray spectrum is also computed (Figure 2b).
the SRIM2013 MC codes [36]. SRIM2013 cannot calculate the range of H3O + ions so the stopping power of 19 O was used instead. The calculations of projected range and range straggling of 19 O + ions are plotted as a function of the kinetic energy below 100 eV ( Figure  1 , which are calculated by fitting based on LSS theory [36,37]. The calculation by SRIM of the nuclear stopping power, which is dominant at low energies, is more accurate than that by the LSS theory [38]. We assumed that those ions are one of the main emitted charged species for high-energy projectiles [8,9,25,39]. The possible emission of neutral species is of course not addressed by the present model based on Coulomb interactions. Multiple ionizations in tracks may reduce the yield of H3O + ions for heavy ions with a higher stopping power [40] which are not taken into account in the present case for the lower density of ionizations produced by protons in tracks. Benchmark of RITRACKS code is performed against the radial distribution of deposited energy by ionization and electronic excitation for 1-MeV protons (p + ) in water and number of ionizations per primary proton (Figure 2a) with reference to literature data [2,3]. The primary δ-ray spectrum is also computed (Figure 2b).   (Figure 3). Spatial distributions of ionizations are displayed for three randomly selected proton impact events. The time evolutions of the real (ε') and imaginary (ε") parts of the dielectric constant show that the electric field is relaxed after~100 fs [41,42]. Within 100 fs, a kinetic energy up to 20 eV is imparted to H 3 O + ions which move less than 6.2 nm during this time step. The time step of ∆t = 100 fs is large enough to calculate the kinetic energy based on the electric potential without consideration for time-dependent kinematics. The real part of the dielectric constant (ε') is taken for calculating the electrostatic potential and electrical forces exerted on the ionized atoms for such a short time scale. The total yields (Y tot ) of produced and emitted H 3 O + ions are computed as a function of depth. The sputtering yield (Y) is also computed as a function of the emission angle (θ) for variable incidence angle (α).   (Figure 3). Spatial distributions of ionizations are displayed for three randomly selected proton impact events. The time evolutions of the real (ε') and imaginary (ε'') parts of the dielectric constant show that the electric field is relaxed after ~100 fs [41,42]. Within 100 fs, a kinetic energy up to 20 eV is imparted to H3O + ions which move less than 6.2 nm during this time step. The time step of t = 100 fs is large enough to calculate the kinetic energy based on the electric potential without consideration for time-dependent kinematics. The real part of the dielectric constant (ε') is taken for calculating the electrostatic potential and electrical forces exerted on the ionized atoms for such a short time scale. The total yields (Ytot) of produced and emitted H3O + ions are computed as a function of depth. The sputtering yield (Y) is also computed as a function of the emission angle (θ) for variable incidence angle (α). The time scales are summarized as below in Figure 4. The plasmon decay time is ~1-10 fs and exciton recombination time is ~10-100 fs. The cooling time of electronic cascade down to 1-keV kinetic energy is ~10-100 fs, and the cooling time from 1 keV to thermal energy is ~40 fs. The ion penetration takes place in ~100 fs, which is the end of electronic cascade and the time t = 0 for our calculations. The electrical field relaxation is achieved for ~200 ps, whereas the Debye relaxation time is ~1 ps and phonon emission time is ~1-10 ps.  The time scales are summarized as below in Figure 4. The plasmon decay time is 1-10 fs and exciton recombination time is~10-100 fs. The cooling time of electronic cascade down to 1-keV kinetic energy is~10-100 fs, and the cooling time from 1 keV to thermal energy is~40 fs. The ion penetration takes place in~100 fs, which is the end of electronic cascade and the time t = 0 for our calculations. The electrical field relaxation is achieved for~200 ps, whereas the Debye relaxation time is~1 ps and phonon emission time is~1-10 ps.
time-dependent kinematics. The real part of the dielectric constant (ε') is taken for calculating the electrostatic potential and electrical forces exerted on the ionized atoms for such a short time scale. The total yields (Ytot) of produced and emitted H3O + ions are computed as a function of depth. The sputtering yield (Y) is also computed as a function of the emission angle (θ) for variable incidence angle (α). The time scales are summarized as below in Figure 4. The plasmon decay time is ~1-10 fs and exciton recombination time is ~10-100 fs. The cooling time of electronic cascade down to 1-keV kinetic energy is ~10-100 fs, and the cooling time from 1 keV to thermal energy is ~40 fs. The ion penetration takes place in ~100 fs, which is the end of electronic cascade and the time t = 0 for our calculations. The electrical field relaxation is achieved for ~200 ps, whereas the Debye relaxation time is ~1 ps and phonon emission time is ~1-10 ps.  This time scale shows that the total cooling time of the electron cascade (~150 fs) to the cutoff energy of 2 eV is longer than the time step used in the simulations, ∆t = 100 fs. Moreover, the range of δ-rays is larger than the track radius of~10 nm (Figure 2a) in our simulations. This means that the electron gas is not thermalized at this stage and spread Quantum Beam Sci. 2023, 7, 7 7 of 16 away from the track core during ∆t. Simulations of ionization cascades in both liquid water and amorphous water ice actually found that the number of secondary electrons increases as a function of time up to a saturation in~100 fs with only small differences between the liquid and solid phase [30]. A backward motion of the thermalized electrons to the track core might take place on a longer time scale because they are strongly trapped, as discussed below.
The mobility of the solvated electron (e aq − ) in amorphous water ice can be considered as very low due to self-trapping into polaronic states on polar H 2 O molecule clusters [43]. Two-photon pump probe laser excitation of water ice has shown that electron localization occurs in~100 fs and self-trapping in~1 ps [44]. Calculation of the electron transport in liquid water was computed by a time-dependent MC method [45]. The temporal distribution of the spatial probability distribution of secondary electrons showed that electrons are trapped at the radial distance of 10 nm from the parent cations within 300 fs [45]. Trapping on localized levels in ice may even take place within a shorter time (~20 fs) for two-photon pulse laser excitation [46]. As a result, it is most unlikely that the ejected electrons would travel back to the track core from where they were kicked out during ∆t. Therefore, the Coulomb field inside the track is not cancelled out on this time scale, and the produced H 3 O + ions are not neutralized or shielded by the negative species which are sitting far away outside the tracks. The processes of neutral species production are not in the scope of this paper, as explained above.
Therefore, our assumption that the electric field remains unchanged until the ions are fully accelerated is justified. The H 3 O + ions are accelerated in the timescale of 100 fs whereas it takes 100 fs-1 ps for the electric field to be relaxed owing to the trapping of electrons by the polarization of water molecules. RITRACKS can calculate transport of secondary electrons but owing to this time sequence, in which the secondary electrons are scarcely transported in the first 100 fs, the electrons were not taken into account to calculate the electric field and the consequent ion kinetic energy because they remain in the penumbra of the proton track.
Violation of the energy conservation is checked carefully. The deposited energy is distributed amongst four channels: (i) the emitted ions for E c − ∆E ≥ E b , (ii) the hot trapped ions in the target for E c − ∆E < E b , where E c and E b are the kinetic energy and binding energy of H 3 O + ions, respectively, (iii) the energy deposited during the transport of ions (∆E), and (iv) the δ-rays that have been ejected from atoms inside the target. This last channel of energy deposition is stored in the attractive Coulomb interaction of the δ-rays with the ionized atoms in the track. The spatial coordinates of the ionizations are randomly sampled based on the ionization cross sections; therefore, ions are occasionally produced within a small distance yielding an electric potential in the order of GeV. Apparently, such events violate the energy conservation law and they are disregarded. This corresponds to the mean distance between water molecules, i.e., 0.31 nm.
The energy stored in the ion path is the sum of those contributions. The second and third channels are released as thermal energy and phonon emission for ∆t ≥ 1 ps. This will eventually induce a thermal spike effect for these longer time scales. As will be further seen in our results, the trapped ionized atoms in the second channel are produced with high kinetic energies.

Results
The kinetic energy and velocity spectra of sputtered H 3 O + ions integrated over the whole target volume are plotted for α = 60 • (Figure 5a,b). Broad in-depth (initial) and surface distributions of ions are observed as a function of E c with a peak at low energy for the initial energy spectrum (dotted curve, Figure 5a). A clear downward shift of the energy and velocity distributions is seen from the initial distribution to the surface distribution due to slowing down of ions from the depth of target to the surface (Figure 5b). The related depth dependence of the emission efficiency or yield per incoming proton and per Å for α = 60 • also exhibits a strong decay for the emitted ions whereas it is constant for the Quantum Beam Sci. 2023, 7, 7 8 of 16 produced ions (initial depth distribution, dotted curve, Figure 5c). No ion emission occurs beyond a depth~10 nm. The total yield of H 3 O + ion emission deduced from Figure 5a is of Y tot = 0.149 (ion/primary proton). Given the high range of kinetic energies (Figure 5a), this emission yield is weakly depending on the selected value of E b . Hence, the overall behavior of energy spectra and angular distribution will not be very much impacted by this value. Note that larger sputtering yields of ice at 10 K were measured for 95-MeV Xe ions, yet for a much higher stopping power of S e~8 × 10 3 keV µm −1 [47].
The angular distribution of yield shows a typical jet-like anisotropic emission in the normal direction (θ = 90 • ) for the normal incidence (α = 90 • ) ( Figure 6). The plotted yields show the integrated yield over a polar angle variation ∆θ corresponding to a ribbon detector with radius of 2 cm and width of 1 cm (Figure 7), like in experiments [21,22] Positive ions in the track are expelled perpendicularly to the target surface. The peak intensity decreases for decreasing α to 60 • and is completely blurred out for near grazing incidence, i.e., α = 19 • (Figure 6). A sharp Lorentzian profile is used for fitting the central part of this peak for α = 90 • , and a Gaussian profile with a standard deviation of~50 • for α = 60 • . Those fits are only meant to show the different shapes of the angular distributions. The width of these angular distributions is clearly larger than the angular straggling (≤15 • ) of emitted H 3 O + ions (Figure 1b). The full θ-ϕ 2D angular distributions are displayed for α = 0 • , α = 19 • , and α = 60 • (Figure 8a-c), where θ is the polar angle, α is the incident angle of the primary ion beam with respect to the target surface, and ϕ is the azimuthal angle in the target plane (defined in Figure 7). It is seen that the jet-like emission is only occurring for α = 90 • . For the latter incidence value of 90 • , a high density of blue dots is found near θ = 90 • , corresponding to sputtered H 3 O + ions collected on the ribbon detector (Figure 8c), whereas for α = 19 • , the emission is directed toward θ = 0 • in the target plane (Figure 8a). An intermediate behavior is found for α = 60 • , in agreement with the angular distributions ( Figure 6). and velocity distributions is seen from the initial distribution to the surface distribution due to slowing down of ions from the depth of target to the surface (Figure 5b). The related depth dependence of the emission efficiency or yield per incoming proton and per Å for α = 60° also exhibits a strong decay for the emitted ions whereas it is constant for the produced ions (initial depth distribution, dotted curve, Figure 5c). No ion emission occurs beyond a depth ~10 nm. The total yield of H3O + ion emission deduced from Figure 5a is of Ytot = 0.149 (ion/primary proton). Given the high range of kinetic energies (Figure 5a), this emission yield is weakly depending on the selected value of Eb. Hence, the overall behavior of energy spectra and angular distribution will not be very much impacted by this value. Note that larger sputtering yields of ice at 10 K were measured for 95-MeV Xe ions, yet for a much higher stopping power of Se ~8 × 10 3 keV µ m −1 [47]. The angular distribution of yield shows a typical jet-like anisotropic emission in the normal direction (θ = 90°) for the normal incidence (α = 90°) ( Figure 6). The plotted yields show the integrated yield over a polar angle variation Δθ corresponding to a ribbon detector with radius of 2 cm and width of 1 cm (Figure 7), like in experiments [21,22] Positive ions in the track are expelled perpendicularly to the target surface. The peak intensity near θ = 90°, corresponding to sputtered H3O + ions collected on the ribbon detector ( Figure  8c), whereas for α = 19°, the emission is directed toward θ = 0° in the target plane ( Figure  8a). An intermediate behavior is found for α = 60°, in agreement with the angular distributions ( Figure 6).   Figure 7. Geometry of the sputtered ion ribbon detector after references [21,22]. Ions impact detector ribbon are counted with identification of the detection angle (θ): θ is the polar an tween the target surface and the vector from the ribbon detector center to the detection poin ribbon, α is the incidence angle of the primary ion beam with respect to the target surface, a the azimuthal angle in the target plane. Figure 7. Geometry of the sputtered ion ribbon detector after references [21,22]. Ions impacting the detector ribbon are counted with identification of the detection angle (θ): θ is the polar angle between the target surface and the vector from the ribbon detector center to the detection point on the ribbon, α is the incidence angle of the primary ion beam with respect to the target surface, and ϕ is the azimuthal angle in the target plane. Figure 7. Geometry of the sputtered ion ribbon detector after references [21,22]. Ions impacting the detector ribbon are counted with identification of the detection angle (θ): θ is the polar angle between the target surface and the vector from the ribbon detector center to the detection point on the ribbon, α is the incidence angle of the primary ion beam with respect to the target surface, and φ is the azimuthal angle in the target plane.

Coulomb Spike Sputtering of Amorphous Water Ice
The present simulations show that ion sputtering induced by ionizations does occur during the time period in which the electric field is retained as explained above. During this short time step, a sufficient kinetic energy and momentum are transferred to the target ions to escape out from the surface. The kinetic energy spectra (Figure 5a) are in the same range as that of H3O + ion emission with a sharp peak at ~35 eV, showing supra-thermal emission induced by 100-keV Ar ion irradiation at 45° incidence angle on an ice film of 420-nm thickness at 80 K [25]. Even though there is some contribution of sputtering induced by nuclear collisions in the latter case (~1 atom per incident ion, as computed with

Coulomb Spike Sputtering of Amorphous Water Ice
The present simulations show that ion sputtering induced by ionizations does occur during the time period in which the electric field is retained as explained above. During this short time step, a sufficient kinetic energy and momentum are transferred to the target ions to escape out from the surface. The kinetic energy spectra (Figure 5a) are in the same range as that of H 3 O + ion emission with a sharp peak at~35 eV, showing supra-thermal emission induced by 100-keV Ar ion irradiation at 45 • incidence angle on an ice film of 420-nm thickness at 80 K [25]. Even though there is some contribution of sputtering induced by nuclear collisions in the latter case (~1 atom per incident ion, as computed with the SRIM2013 code), this cannot account for such a supra-thermal ion emission in the latter case.
The jet-like effect decreases for higher ion incidence angle ( Figure 6), in agreement with experimental data on several insulating materials [5,21,22,28]. It is to be noted that, for such short time scales, neither thermally activated atomic motion can occur, nor any phase change such as melting or boiling, unlike for the i-TS model. No such assumptions are made since phonon emission/absorption processes are not required. Likewise, the major role of secondary electrons was shown for desorption of H 2 O molecules upon fs laser irradiation of liquid water, with a minor effect of phonons [48].

Role of Target Properties
Some challenging questions remain on the role of the material's properties in this ion emission process. Actually, the jet-like anisotropic ion sputtering occurs for crystalline oxides, such as UO 2 , or LiF [21,22], which are not amorphizable, but not for an amorphous one, such as a-SiO 2 [29]. However, this behavior also occurs for an amorphizable oxide such as Y 3 Fe 5 O 12 . The analysis based on a thermal effect is thus questionable for these refractory solids with similar high melting temperatures but showing quite different behaviors. The choice of E b is indeed an important issue since it may govern the yield. However, this does not change drastically the overall behavior, except for the value of Y tot . Surface charge effects were put forward to account for the narrow jet-like emission for LiF [49] and UO 2 [50].
The screening and/or neutralization of the positive ionic charge by the ejected secondary electrons are also important issues. As mentioned above, this is quite unlikely in the case of water ice owing to the strong trapping of ejected electrons on water molecule clusters [48], far outside the track cores. On the basis of X-ray emission spectroscopy data on a SiO 2 aerogel irradiated with 11.4 MeV u −1 Ca ions, authors have claimed that core hole neutralization occurs on a time scale~1 fs [51]. However, this strongly depends on the electron mobility which can be quite low in dielectric materials. This effect is certainly enhanced upon room temperature irradiations. Moreover, positive hole mobility is always very low in any case. The case of a-SiO 2 is quite puzzling since fast neutralization was claimed by some authors on such a short time scale [51].
To summarize, based on simple approximations, the present simulations show that a directional ion sputtering with high-velocity ion emission from the surface can occur due to a Coulomb spike effect for very short time scales. No thermal and/or hydrodynamic effects are taking place at such short time scales. Thermal spike may occur for ∆t > 1 ps and generate an isotropic sputtering background. Such a model can also serve as an input for the TS model by considering the ions with high kinetic energies trapped inside the target material. These ions will be subsequently thermalized for ∆t > 1 ps and induce heating of the target material surface, as explained above in section II. The kinetic energy distribution of these ions corresponds to a temperature of 3 × 10 4 K by using a Maxwell-Boltzmann statistical distribution ( Figure 9). However, accurate fitting with such a kind of distribution is not feasible since these ions are not thermalized at this time step as for the pump-probe laser sputtering of Al 2 O 3 [23]. Our model may be useful for understanding processes at longer time scales, such as hillock formation on surfaces of oxides after swift heavy ion-irradiation [52]. The extension of this study to the sputtering of ionic crystals induced by highly ionizing swift heavy ion irradiation is contemplated as a next step.  We are aware that our modelling, which is quite rough and simple, must be considered as preliminary. However, these MC simulations can serve as an input for a self-consistent MD simulation that calculates both the action of the Coulomb forces on the ions and the stopping of the ions by the neutral molecules in the frozen ice.

Conclusions
MC simulations of ion sputtering in amorphous water ice are carried out on the basis of the track-structure calculation by the RITRACKS computer code under 1-MeV protons for a time scale of 100 fs. Ionized atoms are emitted from a depth of ~10 nm to the target surface by taking into account the energy loss in the matter. A directional emission of We are aware that our modelling, which is quite rough and simple, must be considered as preliminary. However, these MC simulations can serve as an input for a self-consistent MD simulation that calculates both the action of the Coulomb forces on the ions and the stopping of the ions by the neutral molecules in the frozen ice.

Conclusions
MC simulations of ion sputtering in amorphous water ice are carried out on the basis of the track-structure calculation by the RITRACKS computer code under 1-MeV protons for a time scale of 100 fs. Ionized atoms are emitted from a depth of~10 nm to the target surface by taking into account the energy loss in the matter. A directional emission of positively charged H 3 O + species is found in the angular distribution normal to the target surface for normal incidence of 1-MeV p + due to the Coulomb field effect. This anisotropic jet-like emission is blurred out for oblique incidence of projectiles. Our calculation based on track-structure calculation suggests that the directional emission of charged species from irradiated materials is a result of the electric repulsion among the atoms ionized by the incident radiation. Moreover, our calculation showed that trackstructure calculation, originally developed for radiobiology, is a unique and useful tool to understand the phenomena attributed to the electronic and molecular configurations in the material exposed to radiation.
One of the major limitations of the method used in this study is the target material. Owing to this limitation, the calculation was performed for liquid water. By using the generalized track-structure calculation code [53], the idea of this study can be applied to other materials and verified with greater variety of experimental data.

Data Availability Statement:
The data that support the findings of this study are available from the corresponding author upon reasonable request.