Interstellar ice embedded glycine response to H⁺/proton irradiation. A theoretical study

Abstract

The glycine molecule is the simplest amino acid molecule which is expected to be easily formed on ice grains in the interstellar medium. Assessing its stability under interstellar irradiation is therefore crucial in a prebiotic context. In this work we investigate the ionisation mechanism of ice-embedded glycine under H⁺ irradiation at kinetic energies of 1, 10, and 100 keV by studying charge and energy deposition dynamics. To this aim, real time-time dependent-density functional theory combined with a quantum mechanics/molecular mechanics approach is used. Charge distribution analysis reveals ebb-and-flow effects at both molecular and atomic scales, with collisions at nuclei inducing stronger polarization and dipole dynamics than impacts along chemical bonds. While electronic excitation occurs locally, glycine–water interactions can collectively enhance, reduce or leave unchanged the energy deposition, depending on the projectile trajectory and kinetic energy. The charge transfer between the glycine molecule and its surrounding environment does not exceed 0.4 electrons once the projectile is gone. When the projectile does not directly intersect the glycine backbone, the deposited energy on the glycine molecule is negligible and glycine remains neutral. Consequently, the degree of molecular perturbation is governed not only by the projectile's energy, but also by its microscopic trajectory and the surrounding ice environment.

---

1 Introduction

Ionising radiation as cosmic rays (CRs) are believed to initiate a complex and diverse chemistry in space. They ionise H₂ forming H₂⁺ which can initiate an efficient gas phase ion–molecule chemistry¹ in the interstellar dense clouds. They can also dissociate H₂ producing UV photons that can break apart simple molecules within interstellar ices, generating reactive radicals which can recombine to form complex organic molecules (COMs). Such a mechanism has been proposed to explain the potential formation of glycine—the simplest amino acid—under interstellar conditions.^2,3 Despite extensive efforts, glycine has yet to be detected in the gas phase of interstellar environments. An initial claim of detection was made in 2003,⁴ but this was later refuted.⁵ Nevertheless, glycine has been identified in extraterrestrial materials, including cometary dust from Comet 81P/Wild 2,⁶ the coma of Comet 67P/Churyumov–Gerasimenko,⁷ and more recently, in returned samples from the Ryugu⁸ and Bennu⁹ asteroids. These findings support the hypothesis that glycine formed in interstellar ices have survived incorporation into small Solar System bodies. However, they also prompt important questions regarding the molecule's stability and resistance to degradation under exposure to ionising radiation.

Various studies have considered the effects of ionising radiation on interstellar COMs. Comprehensive astrochemical modelling covering gas-phase and solid state chemistry in the interstellar medium, including processes enhanced by high-energy collisions, are available in Arumainayagam et al.¹⁰ and Renoud et al.¹¹ references. Experiments and theoretical studies simulating the effects of cosmic rays on interstellar ices have also been reported,^11–14 demonstrating the formation and evolution of COMs under astrophysically relevant conditions. Among specific studies focusing on glycine, Maclot et al.^15,16 examined its gas-phase fragmentation following collisions with low- and high-energy multiply charged ions. Furthermore, the irradiation of glycine in a water environment has been investigated experimentally in references.^17–19 Additionally, Portugal et al.²⁰ investigated the irradiation of crystalline glycine using Ni¹¹⁺ ions, shedding light on its stability and decomposition pathways under energetic ion bombardment.

On the other hand, theoretical studies in this area remain limited. The RT-TD-DFT ( functional theory) method is widely regarded as the most suitable approach for studying high-energy phenomena, particularly those that induce ultrafast electron rearrangement in colliding systems. This method is particularly effective for investigating how energy is deposited into a molecule during a collision with a high-energy particle, and how this energy is subsequently dissipated, often through electron emission. To date, only two studies have applied this method to astrochemical scenarios of glycine irradiation. The first one, conducted by Parise et al.,²¹ explored the response of a glycine molecule in the gas phase upon collision with a high-energy alpha particle (He²⁺), demonstrating the ionisation of glycine. The second study focused on glycine embedded in an interstellar ice, investigating its degree of ionisation and fragmentation following a collision with a high-energy heavy ion (Ni¹¹⁺) at 46 MeV.²²

Protons account for almost 90% of cosmic rays (CRs) that traverse the interstellar medium. Despite being the least energetic of these particles, protons are the most frequent ionising agents. However studying their impact on molecules experimentally presents significant challenges. In our former work,²² we showed that, among various possible chemical fragmentation pathways, the triply ionised glycine via collision with Ni¹¹⁺, dissociated spontaneously. In the case of collision with a proton, the efficiency of the ionisation dissociation channel is under-investigated. In the present study, we address this question using the RT-TD-DFT approach to gain a deeper understanding of both the ionisation process of glycine by H⁺ and the role played by the surrounding ice environment restricting our analysis to the electronic rearrangement. A proper treatment of the ensuing nuclear response would need non-adiabatic molecular dynamics method implementation, coupling nuclear and electronic motions, i.e. the Ehrenfest approach. Such studies are beyond the scope of this work.

2 Methodology

2.1 Models and computational details

The response of the glycine, embedded in interstellar ice, to the collision with H⁺ was studied with the RT-TD-DFT ²³ approach using the deMon2k software package,²⁴ developers version 6.1.6.²⁵

The ice-embedded glycine model used in this work has been taken from our previous simulations.²⁶ It consists of a neutral canonical glycine molecule (N-Gly) trapped in a non periodic water cluster with 20 Å cubic box, dimensions obtained from the optimization of the low density amorphous (LDA) ice through classical dynamic simulations,²⁷ and named in this paper as N-Gly-LDA. Only the trans–trans (tt) N-Gly isomer was considered in accordance with experimental evidence²⁸ as the tt isomer can easily be formed on interstellar ices. In addition, experimental²⁸ and our theoretical²⁶ investigations showed that this isomer transforms into its zwitterionic structure only at temperatures above 200 K, not relevant for interstellar ice conditions.

We employed an additive quantum mechanics (QM)/classical mechanics (MM) scheme with electrostatic embedding.^29,30 The QM (DFT) region includes the tt N-Gly, surrounded by 62 water molecules, as illustrated in Fig. 1a. The N-Gly structure and the surrounding LDA model (i.e. without the N-Gly molecule) are illustrated in Fig. 1b and c, respectively. The MM region (282 water molecules) was described by the TIP3P force field,^31–34 which includes Coulomb effects inside the MM part. The Lennard-Jones potential was employed to describe the QM-MM and MM interactions and the polarization of the QM region by the MM environment was captured by including the Coulomb interactions between MM atomic charges and QM electronic densities. The glycine QM-MM parameters were taken from the AMBER-FF99SB force field.^31–33 TIP3P QM-MM parameters were assigned to the DFT treated water molecules. Several studies pointed out the importance of using polarizable force fields (MMpol),^23,35,36 nevertheless, in this study we included a wide size QM subsystem, allowing to treat more accurately the interactions between either glycine and water molecules, or between water molecules in the glycine solvation shells. Furthermore, energy deposition and charge migration are expected to be local processes, therefore interactions are likely to occur in the vicinity of the collision area.

Fig. 1 Models used in this study: (a) tt-N-Gly-LDA molecular system with 72 atoms treated at the QM (DFT) level (balls representation) and 282 water atoms (lines representation) treated at the MM level; (b) tt-N-Gly gas phase molecule; (c) surrounding LDA ice model that contains 62 H₂O molecules described with DFT (ball representation) and 282 molecules described with MM (lines representation). The LDA ice model in (c) has an internal cavity due to the elimination of the tt-N-Gly molecule. The atoms are represented with the following color code: O in red, H in white, C in cyan, and N in dark blue.

The Perdew–Becke–Ernzerhof (PBE) exchange–correlation functional³⁷ was employed and the electrons were described with triple-ξ basis functions (TZVP).³⁸ Automatically generated auxiliary basis set, GEN-A2,³⁹ was used to fit the electronic density. The functional and the basis sets were chosen after an extensive acquisition of benchmark results on the deposited energy, because it is sensitive to the chosen basis set.^40,41 The benchmark results are reported in Fig. S1–S4 and their detailed analysis is provided in Section S1 “Benchmark study”.

The self consistent field (SCF) tolerance and the auxiliary density convergence criterion were restricted to 10⁻⁹ and 10⁻⁶ a.u., respectively. An adaptive grid for the exchange correlation (XC) functional numerical integration was selected with a tolerance of 10⁻⁸ a.u. and the charge convergence criterion was adjusted to 10⁻⁶ a.u. The Hirshfeld density partitioning was used to provide charge dynamics analysis.⁴²

2.2 Electron dynamics simulations

The H⁺ colliding particle was modeled as a point charge particle with a constant velocity and a constant charge of +1. Three different irradiation energies E_H⁺ = 1, 10 and 100 keV have been investigated. These projectile kinetic energy values are in the range of the average energies, established for stellar winds and low-energy cosmic rays in the interstellar medium.^10,13,43 At each simulation, the projectile position was initialized 13 Å from the center of mass (CoM) of glycine (placed in the center of the N-Gly-LDA system), ensuring its starting position lies outside the ice box and avoiding any interaction between the projectile and the system at the early stages of the collision process. Due to computational cost limitations, the number of trajectories was reduced to four representative ones covering typical situations in which the proton's trajectory passes through a bond, an atom, or near the glycine molecule. This limited number of trajectories prevents the derivation of statistical properties (e.g., cross sections), but provides insights into the orders of magnitude of the deposited energies and the underlying physical processes such as the glycine ionization. Fig. 2 shows collision trajectories striking the N-Gly molecule at the center of the C–C covalent bond (labelled as CC and colored in orange), at the nitrogen atom of the amino functional group (N, purple color) and at the oxygen atom of the carboxylic functional group (OH, green color). The last chosen collision site, labelled as DISP-5 (see Fig. 2), depicts the projectile trajectory passing through the water-ice environment avoiding a direct collision with the N-Gly molecule. In this latter scenario, H⁺ was displaced by 5 Å from the N-Gly-LDA CoM and the projectile only passes through water molecules. This trajectory has been chosen to capture the electron perturbation caused by H⁺ collision on the water-ice environment and its diffused effects on N-Gly molecules. For easier reading, the suffix “col” will be added to each projectile orientation label, i.e. CCcol, Ncol, OHcol and DISP-5-col, to avoid misunderstanding between notation of a specific atom and collision trajectory.

Fig. 2 Four simulated collision trajectories for H⁺-N-Gly-LDA irradiation: C–C bond (CCcol), nitrogen atom (Ncol), carboxylic oxygen (OHcol), and displaced at 5 Å from the glycine center of mass (DISP-5col), where H⁺ passes only through the water environment. The atoms are represented with the following color code: O in red, H in white, C in cyan, and N in dark blue.

The electron density was propagated to a total of 7 fs dynamics employing the magnus second order propagator coupled to a predictor-corrector scheme with a propagation time step (t_p) of 1 as. The choice of the propagation time step depends on the projectile velocity.²⁵ Several test calculations – t_p = 0.01, 0.1 and 1 as – were performed beforehand. The results, presented in Fig. S5 indicate high similarity between the deposited energy profiles for the three time steps. Therefore, t_p = 1 as was chosen for the present electronic dynamics simulations as a good compromise between accuracy and computational cost.

2.2.1 Energy deposition and charge analysis. The electronic deposit energy (E_d) has been computed according to the following expression:²⁵where E_tot corresponds to the total energy of the system, expressed as the sum of the DFT total energy and the interaction energy between projectile and the N-Gly-LDA system (E_tot = E^DFT + E_proj). E^DFT_scf refers to the ground state energy obtained from the SCF procedure and ϕ_proj, Z_A and R_A correspond, respectively, to the potential generated by the projectile, the atomic charge and the coordinate vector of A nucleus.

With the aim to compare and decompose the irradiation effects on the deposit energy in N-Gly-LDA molecular subsystems, a set of electron dynamics simulations were carried out separately on the gas-phase N-Gly molecule and on the LDA ice model using the same computational and dynamics simulation details as used for the entire system. These trajectories are shown in Fig. S6. The N-Gly structure was extracted from the N-Gly-LDA system without further geometry optimization. Similarly, no initial geometry optimization was performed on the LDA model (in Fig. 1c), so that the proton encounters the same molecules with similar orientation during collision. It must be noted that the cavity created from the N-Gly molecule elimination was preserved and that the level of treatment (QM or MM) of each water molecule was consistent with its description in the full system simulations.

3 Results

In this section, the electronic dynamics within the N-Gly-LDA system induced by different proton collision trajectories are reported. We first address the electron density migrations (also referred to as charge migration), followed by the electronic energy deposition and the examination of the water–glycine collective contributions to E_d.

3.1 Charge migration in N-Gly-LDA system under H⁺ ion collision process

The charge fluctuations were studied for the three H⁺ collision energies, i.e. 1, 10 and 100 keV. The three charge profiles presenting similar behaviors, we only show the results at 10 keV in Fig. 3, the charge evolutions at 1 and 100 keV being reported in Fig. S7 and S8, respectively. For easier reading, only the most significant charge variations are shown in Fig. S9–S11. The subplots a.1, b.1, c.1 and d.1 in Fig. 3 depict the charge exchange between the glycine molecule and its water environment for the four projectile trajectories. Similar behavior is observed for the three collisions involving a trajectory crossing the glycine molecule (a.1, b.1 and c.1). In the initial stage of the simulation, when the proton projectile approaches the glycine molecule, the electrons migrate toward the impact point. After the projectile leaves the glycine region and its Coulomb attraction vanishes in this area, the excess of negative charge on the glycine molecule creates strong Coulomb repulsion that ejects electrons in the ice environment, leaving the glycine molecule positively charged. This ebb-and-flow effect has already been pointed out for other molecules in liquid water impacted by a positive charge (see for instance ref. 35 and 44).

Fig. 3 Charge evolution in N-Gly-LDA molecular system (upper panels a.1, b.1, c.1 and d.1) and of each N-Gly atom (lower panels a.2, b.2, c.2, d.2) under 10 keV H⁺ collision with the four impact sites: (a) CCcol, (b) Ncol, (c) OHcol, and (d) DISP-5col. The color codes are the following: dark green for the carboxylic oxygen (O_H), light red for carboxylic hydrogen (H_O), gray for oxygen, dark brown for carbon from the carboxylic group (C₁), yellow for the α-carbon (bonded to amino group) (C₂), light green and light purple for H₁ and H₂, respectively, purple for N, dark blue and light blue for amino group hydrogens (HN₁ and HN₂).

We now examine the lower panels (a.2, b.2, and c.2), which show the charge fluctuations of each glycine atom over the full simulation time, with insets zooming into the first 2.5 fs. For a collision centered on the C–C bond (CCcol, Fig. 3(a.2)), both carbon atoms attract electrons in a comparable manner. The two hydrogen atoms of the methylene group, H₁ and H₂, display opposite trends, creating a dipole oriented toward the projectile: before impact, H₁ loses electron density while H₂ gains it, and this pattern reverses once the projectile passes the C–C bond. After the projectile exits the glycine region (after around 1 to 2 fs), pronounced intramolecular charge dynamics persist, with atomic charge variations reaching about 0.05 e⁻ (see inset zoom in Fig. 3(a.2)).

In the collision at the nitrogen atom (Ncol, Fig. 3(b.2)), the nitrogen atom accumulates ∼0.9 electrons. Most of this electronic density is drawn from the water molecules (∼0.6 e⁻) and the amino hydrogen HN₂ (∼0.3 e⁻). The projectile gives birth to a dipole dynamics: HN₂ initially acquires electrons, as it lies in the direction of the incoming projectile, before being electronically depleted after the projectile leaves the N-Gly molecule. The opposite behaviour is observed for the second amino hydrogen (HN₁). This dipole oscillation persists throughout the simulation, as reflected in the charge dynamics of HN₁ and HN₂. Notably, the nitrogen atom retains a net electron loss of about 0.15 e⁻ due to the ebb-and-flow effect at the N-atom until the end of the simulation (see inset zoom in Fig. 3(b.2)). Therefore, roughly half of the electrons transferred from glycine to water originate from the nitrogen atom.

The third projectile trajectory through the carboxylic oxygen O_H (OHcol, Fig. 3(c.2)) also induces a significant negative charge accumulation (∼0.8 e⁻) at the impact site O_H. Roughly 0.5 electrons are drawn from the surrounding water, while the remaining charge originates from intramolecular migration, mainly from the other glycine oxygen atom. Once the projectile passes beyond N-Gly, the Coulomb repulsion ejects electrons from the impacted oxygen atom, revealing that the ebb-and-flow behavior occurs at the atomic level in addition to the molecular scale. In the remaining part of the simulation (see inset zoom in Fig. 3(c.2)), charge dynamics with an amplitude of about 0.1 e⁻ are observed predominantly on the carboxyl group as follows from the more pronounced O_H, C₁ and O charge fluctuations.

In the last trajectory (DISP-5col, Fig. 3(d.2)), the projectile does not collide with the N-Gly molecule, instead it passes entirely through the surrounding ice environment. As a result, the charge transfer between the glycine molecule and the nearby water molecules is minimal. It produces charge fluctuations of less than 0.01 e⁻ (see inset zoom in Fig. 3(d.2)), indicating that the projectile's effect in this trajectory is essentially negligible.

The above analysis of the charge migration at 10 keV remains valid across the 1 and 100 keV energy range. The charge evolutions at 1 and 100 keV are shown in Fig. S7 and S8, respectively. The main difference lies in the quantitative magnitude of the accumulated (on LDA) and depleted (on N-Gly) charges. They range from 0 to approximately 0.1 e⁻ at 1 keV (Fig. S7) and from 0 to about 0.4 e⁻ at 100 keV. The higher projectile energy induces a stronger polarization of both N-Gly and LDA. At 100 keV, the dipole dynamics of the impacted glycine atom and its bonded atoms are more pronounced and require a longer time to relax (see Fig. S8). In contrast, a projectile traversing only water (DISP-5col) does not produce significant charge exchange between the glycine molecule and its environment, even at 100 keV (Fig. S8).

3.2 Deposit energy analysis in N-Gly-LDA system under H⁺ collision

Monitoring the electronic potential energy allows to track the energy deposited in the system. The simulations allow to capture the projectile energy transferred to the electronic degrees of freedom only, which is reasonable considering that, for collision energies of ∼10 keV and above, the electronic stopping power dominates over the nuclear stopping power, i.e. the energy is deposited primarily into the electronic subsystem rather than into nuclear motion.^45–47 In other words, only the electron density is perturbed by the energy transferred during irradiation processes while nuclei remained fix along the simulation time.

The energy deposition profile for the specific case of the 10 keV H⁺ CCcol trajectory is represented in Fig. 4 as a function of the projectile position d_projectile. The latter is preferred as it makes the comparison of dynamics performed for different collision energies easier.

Fig. 4 Energy deposition in the N-Gly-LDA system by the 10 keV projectile going through the CC bond (CCcol) as a function of the travelled distance (d_projectile). The snapshots show the position of the projectile (orange ball), corresponding to each peak location extracted at: 0.8 fs (10.5 Å) where ΔE_d(H₂O) ∼ 11.0 eV are deposited in the water environment; 1.0 fs (13.5 Å) when the projectile hits the glycine C–C bond and ΔE_d(Gly) ∼ 17 eV are deposited; 1.3 fs (17.5 Å) where ΔE_d(H₂O) ∼ 10 eV are deposited locally in one of the H₂O molecules colliding with H⁺; 1.5 fs (20.5 Å) where ΔE_d(H₂O) ∼ 1.0 eV are deposited locally in the water molecules.

It can be seen that the energy remains constant below 7 Å and above 23 Å, which corresponds to subparts of the trajectory where the projectile travels through the MM region. The energy deposition profile has several peaks, which are associated to geometrical configurations where the projectile approaches the molecules, as shown in the snapshots in Fig. 4. The second peak at 13.5 Å corresponds to the collision with the glycine molecule at the center of the CC bond. The other peaks occur when the projectile collides or passes at the vicinity of water molecules (see snapshots at t = 0.8, 1.3 and 1.5 fs in Fig. 4). We note here that the simulations describe the collective electronic response of the N-Gly-LDA QM region, and are not treated as a sum of individual independent molecule responses. However, the structuring of the deposited energy profile in successive peaks indicates that the energy deposition remains predominantly local. This allows for a semi-quantitative estimation of the energy transferred in a given molecule. In particular, we derive the energy deposited in the glycine electronic subsystem ΔE_d(Gly) ∼ 17 eV from the difference between the deposit energy after and before the peak at 13.5 Å (see Fig. 4).

The energy deposition curves are reported in Fig. 5 for all the investigated collision trajectories (CCcol, Ncol, OHcol and DISP-5-col) and projectile kinetic energies (E_H⁺ = 1, 10 and 100 keV). The energy deposition profiles present very well marked peaks for trajectories with 1 and 10 keV collision energies passing through an atom (Ncol and OHcol). A slightly less pronounced peak is observed for the bond crossing trajectory (CCcol), whereas no corresponding peaks are visible in the DISP-5col trajectory. At 100 keV, the energy deposition appears much more continuous and it even becomes difficult to identify peaks corresponding to the energy depositions within the glycine molecule. At 1 keV, the energy deposited in N-Gly (see ΔE_d(Gly) in Table 1) is rather low, in particular with respect to the glycine and water ionisation potential (estimated to be in the 9–10 eV range^48–52 for glycine and about 12–13 eV for water⁵³). Interestingly, about 0.4–0.5 electrons are ejected from the glycine molecule; however, its neutral character is rapidly recovered, with only a minor ebb-and-flow effect resulting in a final net negative charge of about 0.1 for the Ncol and OHcol trajectories in Fig. S7.

Fig. 5 Energy deposition in the N-Gly-LDA molecular system under H⁺ collision for 1 (top), 10 (middle) and 100 keV (bottom). Each panel shows four curves for CCcol (orange line), Ncol (purple), OHcol (green) and DISP-5col (pink), impact sites. In order to simultaneously compare results at three collision energies, the x-axis has been defined as the distance travelled by projectile H⁺ through the N-Gly-LDA slab (as a product of simulation time and projectile velocity).

Table 1 Deposited (transferred) energy in N-Gly (ΔE_d), computed from the deposited energy profiles of the CCcol, Ncol, OHcol and DISP-5col trajectories of H⁺ projectile at 1, 10 and 100 keV. Values in parentheses were computed from simulations of gas-phase N-Gly

Collision energy (keV)	ΔEd(Gly) (eV)
	CCcol	Ncol	OHcol	DISP-5col
1	1.8 (2.1)	2.8 (7.4)	1.3 (5.0)	—
10	15.2 (14.4)	21.9 (28.2)	23.7 (25.0)	—
100	32.2	63.4	62.1	—

---

The energy transferred in the electronic glycine subsystem increases with the projectile collision energy (Table 1). The ΔE_d(Gly) values are larger when the proton impacts an atom (OHcol and Ncol trajectories leading to similar values) rather than a bond (CCcol). At 10 keV for CCcol, Ncol and OHcol trajectories, the deposit energy in the glycine subsystem is larger than the glycine and water ionisation potentials and the glycine molecule carries a positive net charge between 0.2–0.4 at the end of the simulations (see Fig. 3). When the projectile crosses only the water region (DISP-5col trajectories), the deposited energy in the glycine electronic subsystem is so small that it was not even possible to determine it from the deposit energy profiles.

We finally compare the former results with the values of ΔE_d(Gly) obtained from simulations of the gas-phase N-Gly at 1 and 10 keV, reported in Table 1. The energy transferred to the electronic subsystem of N-Gly is larger for the gas-phase molecule than for N-Gly embedded in ice, with the only exception being the CCcol trajectory at 10 keV.

We now examine the additive nature of the energy deposition process, specifically the role of glycine–water collective effects. To this end, an additional set of simulations was performed using the same collision energies and projectile trajectories, but including either only the glycine molecule or only the water environment (QM and MM parts), each fixed at the geometry of the full N-Gly-LDA system. These latter calculations were performed at 1 and 10 keV energies because the deposited energy profiles are better structured (sharper peaks) than at 100 keV. Fig. 6 left and right panels report the energy deposition for these simulations along with their sum and the deposit energy for trajectories across the full system (orange lines). For CCcol trajectories, the energy deposited in the full system is slightly higher than the sum of the energies deposited in the separate subsystems (see panels a in Fig. 6); in other words, the collective response of glycine and water electrons enhances the energy deposition. The opposite behavior is observed for the Ncol trajectory (see panels b in Fig. 6), where collective effects instead reduce the deposited energy. For the OHcol and DISP-5col trajectories (Fig. 6 panels c and d), collective effects have essentially no influence: summing the energies deposited in the isolated glycine and water subsystems provides a good approximation to the energy deposited in the full system. In summary, collective glycine–water interactions may enhance, reduce, or leave unchanged the energy deposition, depending on the projectile trajectory and kinetic energy. The largest effect occurs for the Ncol trajectory at 1 keV, where collective interactions lead to a reduction of ∼35% of the deposited energy.

Fig. 6 Energy deposition decomposition for 1 keV H⁺ (left panel) and 10 keV H⁺ collision (right panel) for: (a) CCcol, (b) Ncol, (c) OHcol, and (d) DISP-5col. N-Gly-LDA molecular system (orange curve), isolated N-Gly molecule (red), LDA bulk collision (blue) and the combination of individual N-Gly + LDA deposited energy components (black dash line).

3.3 General discussion

To qualitatively assess the susceptibility of glycine to radiation-driven decomposition, we remind the previously established dissociation thresholds of 74 kcal mol⁻¹ (3 eV) for C–C cleavage and 52 kcal mol⁻¹ (2.25 eV) for C–N bond breaking computed for a neutral glycine molecule.²² These barriers drop when the glycine molecule is singly or multiply charged,²² consistently with the work of Maclot et al.¹⁵

In a former work,²² we showed that highly charged Ni¹¹⁺ ions at much higher energies (46 MeV) could lead to pronounced ionisation (glycine was shown to reach a triple ionisation state when directly impacted), opening the way to the spontaneous dissociation. The simulations performed in the present work show that fragmentation routes involving multiply charged glycine states can be disregarded for H⁺ collisions within the investigated energy window (1–100 keV) as the glycine was shown to be only partially ionized (reaching a maximum loss of 0.4 e⁻).

The proton impacts at 10 and 100 keV can, in principle, deposit sufficient energy inside the glycine molecule to exceed the neutral dissociation threshold energies (see Table 1). A collision energy of 1 keV appears to be a limit case as only a single trajectory (Ncol) deposited more than the threshold dissociation energy for C–N breaking. At all the investigated energies, the simulations indicate that the amount of deposited energy is governed not only by the incident energy, but also by the specific trajectory of the incoming projectile. When the projectile traverses only the water environment (DISP-5col trajectory), it deposits almost no energy on glycine, as inferred from the electronic deposit energy profiles and ΔE_d(Gly) in Table 1, neither induces notable electronic excitations in glycine or water subsystems as concluded from the charge dynamics analysis in Section 3.1.

One should stress that only the energy deposition processes have been simulated in the present work, preventing drawing conclusions on glycine stability and the ice protecting role. Indeed, crucial points lie in the description of subsequent energy flows between electronic and nuclear degrees of freedom and in the efficiency of (electronic and vibrational) energy spreading between the energy deposition region, i.e. close to the cosmic ray trajectory, and the ice environment, which could be viewed as an energy absorbing bath.

Indeed previous studies have demonstrated that, following electronic excitation, gas phase glycine can undergo intersystem crossing, internal conversion, relaxation through conical intersections, processes that can drive photochemical reactions such as bond cleavage.^54,55 We, therefore, do not exclude the role that could be played by such electronic excited states, if populated, either directly during the energy deposition or through subsequent electronic relaxation. In consequence, it cannot be completely ruled out that glycine molecules embedded in ice can be destroyed by protons with 10 or 100 keV kinetic energies.

4 Conclusions

RT-TD-DFT simulations were carried out to investigate the electronic response of a glycine molecule subjected to collisions with proton cosmic rays at incident energies of 1, 10, and 100 keV and four distinct trajectories. The simulations reveal that proton trajectories intersecting glycine induce a minor net charge transfer from the molecule to the surrounding amorphous ice matrix. This effect is negligible at 1 keV (∼0.1 e⁻) and increases to ∼0.4 e⁻ at 100 keV. Analysis of the charge dynamics reveals an ebb-and-flow effect on migration charge at both molecular scale (between the glycine molecule and its surrounding water environment) and atomic (within the glycine molecule) scale. Collisions of the projectile with nuclei (carboxylic oxygen and nitrogen) induce stronger charge polarization and enhanced dipole dynamics in the affected molecular fragments than collisions occurring along the C–C chemical bond.

Unlike collisions with a highly charged heavy ions Ni¹¹⁺ at 46 MeV, proton trajectories intersecting glycine lead only to weak ionisation, making fragmentation through multiply charged states unlikely. However, since this work only considers the energy deposited along four distinct projectile trajectories, no definitive conclusions can be drawn regarding the stability of glycine embedded in ice. Furthermore, a subsequent redistribution of energy from the ice environment to the excited states of glycine could further promote bond breaking.

From an astrochemical perspective, these findings, though semiquantitative, carry significant implications. Interstellar ices, that are known reservoirs of simple and complex organic species, are continuously processed by cosmic rays. If glycine remains neutral, the presence of high energy dissociation barriers requires that a significant amount of the deposited energy can be finally located in a dissociation mode, and one could expect that, before this happens, the deposited energy is rapidly redistributed towards the ice environment. Glycine multiple ionisation offers a way to lower barriers on dissociation channels but such multiple ionisation was not observed in H⁺ collisions in the 1–100 keV energy range and, in the case of Ni¹¹⁺ ions at much higher energies (46 MeV) only when the projectile passes through the molecule. These clues are in line with the scenario proposed by Portugal et al.,²⁰ in which glycine embedded within cold interstellar ice mantles may retain chemical stability against cosmic ray processing over extended timescales. The present study strengthens the view that water-rich ices could enhance the survival probability of glycine during exposure to astrophysical radiation fields.

Author contributions

MYB wrote the draft and carried out dynamics simulations, TM edited the draft and carried out part of the dynamics simulations. MR and DT edited the draft. TM, DT and MR supervised the project. The manuscript was written with contribution from all the co-authors. All the co-authors read and approved the manuscript.

Conflicts of interest

There are no conflicts to declare.

Data availability

All data will be given by the authors upon request.

Supplementary information (SI): benchmark study on basis functions and simulation time-step; supplementary models; distribution of charge migration at E_kin = 1, 10 and 100 keV. See DOI: https://doi.org/10.1039/d6cp00651e.

Acknowledgements

All authors acknowledge the Agence Nationale de la Recherche (ANR) (project ANR-19-CE29-0011-02 RUBI). TM acknowledges access to the HPC resources of CCRT/CINES/IDRIS, which was granted under the allocation A0130807369 by GENCI. MYB acknowledge HPC mesocenter CALMIP (UAR CNRS 3667, project p18009) for allocation of computer resources. All authors address special and kind thanks to Dr A. de la Lande for providing deMon2k software package developers version 6.1.6 and a complete tutorial of usage.

References

1. I. R. Cooke and I. R. Sims, ACS Earth Space Chem., 2019, 3, 1109–1134.
2. R. T. Garrod, Astrophys. J., 2013, 765, 60.
3. T. Suzuki, L. Majumdar, M. Ohishi, M. Saito, T. Hirota and V. Wakelam, Astrophys. J., 2018, 863, 51.
4. Y.-J. Kuan, S. B. Charnley, H.-C. Huang, W.-L. Tseng and Z. Kisiel, Astrophys. J., 2003, 593, 848.
5. L. E. Snyder, F. J. Lovas, J. M. Hollis, D. N. Friedel, P. R. Jewell, A. Remijan, V. V. Ilyushin, E. A. Alekseev and S. F. Dyubko, Astrophys. J., 2005, 619, 914.
6. J. E. Elsila, D. P. Glavin and J. P. Dworkin, Meteorit. Planet. Sci., 2009, 44, 1323–1330.
7. K. Altwegg, H. Balsiger, A. Bar-Nun, J.-J. Berthelier, A. Bieler, P. Bochsler, C. Briois, U. Calmonte, M. R. Combi and H. Cottin, Sci. Adv., 2016, 2, e1600285.
8. C. Potiszil, T. Ota, M. Yamanaka, C. Sakaguchi, K. Kobayashi, R. Tanaka, T. Kunihiro, H. Kitagawa, M. Abe, A. Miyazaki, A. Nakato, S. Nakazawa, M. Nishimura, T. Okada, T. Saiki, S. Tanaka, F. Terui, Y. Tsuda, T. Usui, S.-I. Watanabe, T. Yada, K. Yogata, M. Yoshikawa and E. Nakamura, Nat. Commun., 2023, 14, 1482.
9. D. P. Glavin, J. P. Dworkin and C. M. O. Alexander, et al., Nat. Astron., 2025, 9, 199–210.
10. C. R. Arumainayagam, R. T. Garrod, M. C. Boyer, A. K. Hay, S. T. Bao, J. S. Campbell, J. Wang, C. M. Nowak, M. R. Arumainayagam and P. J. Hodge, Chem. Soc. Rev., 2019, 48, 2293–2314.
11. J. Renoud, S. Indrajith, A. Domaracka, P. Rousseau, P. Moretto-Capelle, B. A. Huber and J.-P. Champeaux, Phys. Chem. Chem. Phys., 2020, 22, 5785–5796.
12. E. Herbst and R. T. Garrod, Front. Astron. Space Sci., 2022, 8, 789428.
13. H. Rothard, A. Domaracka, P. Boduch, M. E. Palumbo, G. Strazzulla, E. F. Da Silveira and E. Dartois, J. Phys. B: At., Mol. Opt. Phys., 2017, 50, 062001.
14. F. Matias, N. E. Koval, P. de Vera, R. Garcia-Molina, I. Abril, J. M. B. Shorto, H. Yoriyaz, J. J. N. Pereira, T. F. Silva, M. H. Tabacniks, M. Vos and P. L. Grande, Phys. Rev. Lett., 2025, 135, 148003.
15. S. Maclot, D. G. Piekarski, A. Domaracka, A. Méry, V. Vizcaino, L. Adoui, F. Martn, M. Alcam, B. A. Huber, P. Rousseau and S. Daz-Tendero, J. Phys.Chem. Lett., 2013, 4, 3903–3909.
16. S. Maclot, D. G. Piekarski, R. Delaunay, A. Domaracka, A. Méry, V. Vizcaino, J.-Y. Chesnel, F. Martn, M. Alcam, B. A. Huber, L. Adoui, P. Rousseau and S. Daz-Tendero, Eur. Phys. J. D, 2014, 68, 149.
17. P. A. Gerakines, R. L. Hudson, M. H. Moore and J.-L. Bell, Icarus, 2012, 220, 647–659.
18. S. Pilling, L. A. Mendes, V. Bordalo, C. F. Guaman, C. R. Ponciano and E. F. da Silveira, Astrobiology, 2013, 13, 79–91.
19. B. Maté, I. Tanarro, R. Escribano, M. A. Moreno and V. J. Herrero, Astrophys. J., 2015, 806, 151.
20. W. Portugal, S. Pilling, P. Boduch, H. Rothard and D. P. P. Andrade, Mon. Not. R. Astron. Soc., 2014, 441, 3209–3225.
21. A. Parise, A. de la Lande, T. Marino and N. Russo, Phys. Life Rev., 2020, 32, 114–116.
22. D. Talbi, S. Christodoulou, M. Yusef-Buey, M. Rapacioli and T. Mineva, Life Sci. Space Res., 2026, 49, 88–93.
23. X. Wu, J.-M. Teuler, F. Cailliez, C. Clavaguéra, D. R. Salahub and A. de la Lande, J. Chem. Theory Comput., 2017, 13, 3985–4002.
24. G. Geudtner, P. Calaminici, J. Carmona-Espíndola, J. M. del Campo, V. D. Domínguez-Soria, R. F. Moreno, G. U. Gamboa, A. Goursot, A. M. Köster and J. Reveles, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2012, 2, 548–555.
25. D. Tolu, D. Guillaumont and A. de la Lande, J. Phys. Chem. A, 2023, 127, 7045–7057.
26. M. Yusef-Buey, T. Mineva, D. Talbi and M. Rapacioli, Phys. Chem. Chem. Phys., 2024, 26, 2414–2425.
27. P. Ghesquière, Rôle des glaces interstellaires dans la complexitémoléculaire de l'espace: modélisation par les méthodes de la chimie théorique, 2015.
28. S. Ioppolo, G. Fedoseev, K. J. Chuang, H. M. Cuppen, A. R. Clements, M. Jin, R. T. Garrod, D. Qasim, V. Kofman, E. F. van Dishoeck and H. Linnartz, Nat. Astron., 2021, 5, 197–205.
29. A. de La Lande, A. Alvarez-Ibarra, K. Hasnaoui, F. Cailliez, X. Wu, T. Mineva, J. Cuny, P. Calaminici, L. López-Sosa and G. Geudtner, Molecules, 2019, 24, 1653.
30. J. D. Samaniego-Rojas, L.-I. Hernández-Segura, L. López-Sosa, R. I. Delgado-Venegas, B. Gomez, J.-C. Lambry, A. de la Lande, T. Mineva, J. Alejandre and B. A. Zúñiga-Gutiérrez, Multiscale Dynamics Simulations: Nano and Nano-bio Systems in Complex Environments. QM/MM with auxiliary DFT in deMon2k, The Royal Society of Chemistry, 2021, vol. 1, pp. 1–54.
31. V. Hornak, R. Abel, A. Okur, B. Strockbine, A. Roitberg and C. Simmerling, Proteins: Struct., Funct., Bioinf., 2006, 65, 712–725.
32. W. D. Cornell, P. Cieplak, C. I. Bayly, I. R. Gould, K. M. Merz, D. M. Ferguson, D. C. Spellmeyer, T. Fox, J. W. Caldwell and P. A. Kollman, J. Am. Chem. Soc., 1995, 117, 5179–5197.
33. W. S. Ross and C. C. Hardin, J. Am. Chem. Soc., 1994, 116, 6070–6080.
34. J. Ȧqvist, J. Phys. Chem. Lett., 1990, 94, 8021–8024.
35. A. de la Lande, S. Denisov and M. Mostafavi, Phys. Chem. Chem. Phys., 2021, 23, 21148–21162.
36. G. Donati, A. Wildman, S. Caprasecca, D. B. Lingerfelt, F. Lipparini, B. Mennucci and X. Li, J. Phys. Chem. Lett., 2017, 8, 5283–5289.
37. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1998, 80, 891.
38. N. Godbout, D. R. Salahub, J. Andzelm and E. Wimmer, Can. J. Chem., 1992, 70, 560–571.
39. A. M. Köster, J. U. Reveles and J. M. del Campo, J. Chem. Phys., 2004, 121, 3417–3424.
40. B. Gu, B. Cunningham, D. Muñoz Santiburcio, F. Da Pieve, E. Artacho and J. Kohanoff, J. Chem. Phys., 2020, 153, 034113.
41. I. Maliyov, J.-P. Crocombette and F. Bruneval, Phys. Rev. B, 2020, 101, 035136.
42. F. L. Hirshfeld, Theor. Chim. Acta, 1977, 44, 129–138.
43. A. M. Bykov, A. Marcowith, E. Amato, M. E. Kalyashova, J. M. D. Kruijssen and E. Waxman, Space Sci. Rev., 2020, 216, 42.
44. A. Alvarez-Ibarra, A. Parise, K. Hasnaoui and A. de la Lande, Phys. Chem. Chem. Phys., 2020, 22, 7747–7758.
45. Y. Yao, D. C. Yost and Y. Kanai, Phys. Rev. Lett., 2019, 123, 066401.
46. C. Race, D. Mason, M. Finnis, W. Foulkes, A. Horsfield and A. Sutton, Rep. Prog. Phys., 2010, 73, 116501.
47. J. F. Ziegler, J. Appl. Phys., 1999, 85, 1249–1272.
48. P. Cannington and N. S. Ham, J. Electron Spectrosc. Relat. Phenom., 1983, 32, 139–151.
49. T. P. Debies and J. W. Rabalais, J. Electron Spectrosc. Relat. Phenom., 1974, 3, 315–322.
50. V. Zaretskii, V. Sadovskaya, N. Wulfson, V. Sizoy and V. Merimson, Org. Mass Spectrom, 1971, 5, 1179–1182.
51. H. J. Svec and G. A. Junk, J. Am. Chem. Soc., 1967, 89, 790–796.
52. M. Slifkin and A. Allison, Nature, 1967, 215, 949–950.
53. F. Ambrosio, Z. Guo and A. Pasquarello, J. Phys. Chem. Lett., 2018, 9, 3212–3216.
54. C. Guerra, Y. A. Rodríguez-Núñez and A. E. Ensuncho, Chem. Phys. Chem., 2024, 25, e202300655.
55. J. Nirasok, P. Panajapo, P. Promma, P. Suwannakham and K. Sagarik, J. Photochem. Photobiol., A, 2023, 436, 114354.

---

Footnotes

† Current affiliation: Laboratoire de Recherche Spécialisé dans l'Analyse et l'Architecture des Systémes (LAAS-CNRS), Toulouse, France. E-mail: maysa.yusef-buey@laas.fr

‡ These authors contributed equally to this work.

---