# Ultrafast proton transfer in a photoionized glycine by a mixed quantum–classical and quantum dynamics

Kossi Kety, Jesús González-Vázquez, Piero Decleva, Jorge Delgado, Alicia Palacios, Fernando Martín, Daniel Peláez, Fabien Gatti, Raluca Cireasa, Loïc Joubert-Doriol, et al.

► To cite this version:

Kossi Kety, Jesús González-Vázquez, Piero Decleva, Jorge Delgado, Alicia Palacios, et al. Ultrafast proton transfer in a photoionized glycine by a mixed quantum–classical and quantum dynamics. The Journal of Chemical Physics, 2025, 163 (16), pp.164122. 10.1063/5.0288179. hal-05365978

HAL Id: hal-05365978

https://hal.science/hal-05365978v1

Submitted on 14 Nov 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Ultrafast Proton Transfer in a Photoionized Glycine by a Mixed Quantum-Classical
and Quantum Dynamics

Kossi Kety,¹ Jesus González-Vázquez,² Piero Decleva,³ Jorge Delgado,² Alicia

Palacios,² Fernando Martín,² Daniel Peláez,⁴ Fabien Gatti,⁴,⁵ Raluca Cireasa,⁴ Loïc
Joubert-Doriol,¹ and Alexander O. Mitrushchenkov¹, a)

1) Université Gustave Eiffel, Laboratoire de Modélisation et Simulation Multi-échelle MSME, F-77454 Marne-la-Vallée, France

²) Departamento de Química, Modulo 13, Universidad Autónoma de Madrid, 28049 Madrid, Spain.

³) CNR IOM and Dipartimento di Scienze Chimiche e Farmaceutiche,
Università di Trieste, 34127 Trieste, Italy

⁴) Université Paris-Saclay, Institut des Sciences Moléculaires d'Orsay ISMO, UMR CNRS 8214, F-91405 Orsay, France

⁵) Université Bourgogne Europe, CNRS, Laboratoire Interdisciplinaire Carnot de Bourgogne ICB UMR 6303, F-21000 Dijon, France

(Dated: 14 November 2025)
We have theoretically investigated the ultrafast intramolecular hydrogen transfer in the glycine molecule after ionization, as observed by R. Cireasa *et al*. [The journal of physical chemistry letters 9 (2018) 6012-6016], following excitation of the molecule with an XUV attosecond pulse train of 1.5 fs duration. In the experiment, the interaction of the glycine molecule with the XUV pulse creates a superposition of electronic states, whose dynamics is coupled to the nuclear one.

We employed the static exchange RAS-DFT (Restricted Active Space Density Functional Theory) correlated approach, as implemented in the Tiresia code [P. De-cleva et al. Continuum electronic states: The tiresia code. Molecules, 27(6):2026, 2022.], to evaluate ionization probabilities. Coherence effects were studied through quantum dynamics simulations using the Multi-Layer Multi-Configuration Time-Dependent Hartree (ML-MCTDH) method on a vibronic coupling Hamiltonian model. Our findings indicate that, for the pulses used in Cireasa et al experiment, electronic coherence dissipates very rapidly, in less than 3 fs. Consequently, we performed simulations starting from single electronic states. In addition, we described the long-term coupled electron-nuclear dynamics using the Trajectory Surface Hopping (TSH) method. Our results reveal that hydrogen transfer predominantly occurs when the active state reaches the cationic ground state. Charge analysis confirms that this process corresponds to a proton transfer.

<footnote>a) Electronic mail: Authors to whom correspondence should be addressed: fabien.gatti@u-bourgogne.fr, raluca.cireasa@universite-paris-saclay.fr, alexander.mitrushchenkov@univ-eiffel.fr</footnote>
# I. INTRODUCTION

With the advent of ultrashort attosecond laser pulses, a new field of Chemistry has emerged, based on the ability to trigger and track electronic motion on its natural time scale. Since changes in the electronic density are responsible for breaking and forming chemical bonds, these new techniques could enable control of chemical processes at their most fundamental level. From the experimental perspective, this is possible by using a the pump-probe scheme. In this regard, the advent of femtochemistry-pioneered by the 1999 Nobel Prize in Chemistry Ahmed Zewail¹- has allowed scientists to follow nuclear motion in real time. As mentioned above, the ability to generate attosecond pulses now makes it possible to track not only nuclear motion but also electronic motion.

In Ref. 2, Cireasa et al. performed a pump-probe experiment on the glycine molecule in the gas phase. These time-dependent measurements revealed the presence of a hydrogen migration process occurring at about 50 fs. They observed that the hydrogen migration can occur following different pathways. More precisely, the molecule was irradiated using an ultrafast light pulse of wavelengths lying in the extreme ultraviolet range (10–100 eV) to induce ionization (estimated overall duration of 1.5 fs). To probe the dynamics, a second infra-red (IR) radiation was employed (estimated 5-fs duration centered at 750 nm). The photo-fragments yields were measured as a function of the time delay between the two pulses. Similar experiments on phenylalanine revealed ultrafast oscillations in the temporal evolution of the photo-fragment yield, interpreted as a signature of electronic coherence.³. For this type of experiments, theoretical modeling is essential to better understand the experimental measurements. The interest in studying the glycine molecule lies primarily in the fact that it is the simplest amino acid, the monomers that build up proteins. This amino acid is often used as a model to interpret the chemical and physical processes that arise in more complex amino acids and to study molecules of biological interest. For example, understanding its photo-ionization is crucial for predicting the effects of oxidative damage in proteins due to ionizing radiation. Additionally, glycine serves as a prototype for studying the fundamental process of proton transfer.

A full quantum description of the interplay between the electronic and nuclear motions is an extremely challenging task, especially if an ionization process and a motion of large amplitude for the nuclei (here the proton transfer) are involved. This work does not aim to
provide a complete quantum description of the aforementioned experiment on the glycine molecule. Rather, when ultrashort pulses are used to initiate a chemical reaction like H⁺/H transfer, we seek to assess the decoherence time of the electronic motion and to observe the H⁺/H transfer.

Several studies have been conducted on the possible isomerization reaction of the glycine radical cation in the gas phase by exploring the potential energy surface using either the a semi-empirical MNDO method⁴ or density functional theory (B3LYP-BHLYP) and coupled cluster CCSD(T) methods⁵⁻⁷. These studies have identified four independent pathways through which isomerization can occur. However, all four pathways require overcoming energy barriers higher than that of the C–C fragmentation channel, making this fragmentation the more favorable process.

These previous works have found that the enol [NH₂ – CH – C(OH)₂] is the most stable isomer of the glycine radical cation. The first study of this reaction using molecular dynamics was conducted by D. Shemesh et al.⁸ on the semi-empirical PM3 potential energy surface. In addition to confirming the mains results already available in the literature, this study observed a hydrogen migration from the amino group to the carbonyl group occurring prior to dissociation approximately within 320 fs. A conformer-specific dynamic was investigated using the same method⁹ revealing a water loss occurring in 3 ps following a hydrogen migration in about 2.97 ps.

It is worth noting that ab initio molecular dynamics simulations on the glycine cation have been conducted in a very different experimental context¹⁰. Specifically, previous simulations investigated the stability of singly charged glycine molecules formed through collisions between neutral glycine and low-energy multiply charged ions in the gas phase. In contrast, in Ref. 2, the glycine molecule is ionized by the pump pulse rather than by ion collisions.

In this work, we seek to determine whether proton (or hydrogen) transfer can also be observed in Ref. 2 and to investigate the duration of the induced electronic coherence. In this context, important simulations have been conducted to describe the multi-electron wavepacket propagation in glycine using an *ab initio method*¹¹. They investigated the response of the electronic cloud to sudden removal of one of its electrons. They observed, in very specific cases (i.e. when the electron is removed from certain orbitals and certain conformers), charge migration occurs within the molecule within a few femtoseconds. In these simulations, the nuclei were held fixed. An important question is whether the corresponding electronic
quantum interferences could be used to guide hydrogen migration. The interplay between the electronic and nuclear motions in glycine has been studied in Ref. 12 using a combination of a multi-reference static-exchange method and a surface hopping approach. Their simulations suggest that electronic coherence associated with subsets of cationic states can be imprinted in specific fragmentation channels. However, this study focused on the dication (double ionization with an XUV probe pulse). In addition, they did not observe any hydrogen migration.

We assume that decoupling of the photoelectron from the ion core is very fast, so that it does not influence the ion dynamics 13, and compute ionization amplitudes in first order perturbation theory. Standard quantum chemistry methods, typically used to describe electronic excited bound states, fail in this context. Instead, approaches capable of accurately describing electronic continuum states are required. To this end, we have employed the static-exchange RAS-DFT approach available in the Tiresia code¹⁴ to model the electronic continuum.

Regarding the coupled electron-nuclear dynamics, due to the very large number of electronic states involved and the size of the molecule which exponentially explodes the computational time, we have worked in two steps. First a quantum simulation based on a linear vibronic coupling model and the Multi-Layer Multi-Configuration Time-Dependent Hartree (ML-MCTDH) method¹⁵⁻¹⁸ including all the nuclear degrees of freedom has been achieved to analyze the electronic decoherence. Second, we used a semi-classical method namely the trajectory surface hopping where the nuclear wave function is represented by a bunch of independent classical trajectories and the electronic motion are described by solving the quantum mechanical Schrödinger equation. This method has already been used with success on the glycine molecule with a focus on the fragmentation pathways¹².

The paper is structured as follows. First, we detail the methodology and computational procedures employed for the MCTDH and semiclassical calculations. Next, we present and discuss the results. Finally, we summarize our findings and outline potential avenues for future research on this and related systems.
## II. COMPUTATIONAL METHODS

### A. Photoionization observables

To accurately evaluate the continuum states of the singly ionized glycine molecule, we employ the single-channel approximation within the close coupling expansion¹⁴,¹⁹ of the wave function, as implemented in the Tiresia code¹⁴. The molecular continuum states are written as an antisymmetrized product of $N-1$ particles bound eigenstates $\Phi_F^{N-1}$, which describe the final state and the continuum orbitals $\phi_{EFlm}$, which must ensure the correct asymptotic conditions of scattering states²⁰:

$$ \Psi_{EFlm}^{N}(\mathbf{r}) = \mathcal{A}[\Phi_{F}^{N-1}(\mathbf{r}_{1}, \dots, \mathbf{r}_{N-1})\phi_{EFlm}(\mathbf{r}_{N})]. \quad (1) $$

In the above expression, $E$ represents the photon-electron energy, while $l$ and $m$ denote the quantum angular momentum numbers. In the correlated single-channel approach, the Dyson orbital is defined as the overlap between the electronic ground-state wave function of the neutral molecule $\Phi_I^N$, and the electronic wave function of the cation, $\Phi_F^{N-1}$, where $I$ stands for the initial state and $F$ for the final state:

$$ \phi_{FI}^{D}(\mathbf{r}_{N}) = \langle \Phi_{F}^{N-1} | \Phi_{I}^{N} \rangle = \quad (2) $$

$$ \sqrt{N} \times \int \Phi_{F}^{*N-1}(\mathbf{r}_{1}, \dots, \mathbf{r}_{N-1}) \Phi_{I}^{N}(\mathbf{r}_{1}, \dots, \mathbf{r}_{N}) d\mathbf{r}_{1} \dots d\mathbf{r}_{N-1}. \quad (3) $$

The many-particle dipole transition moments can then be expressed as:

$$ D_{EFlm\gamma} = \langle \Psi_{EFlm}^{N} | D_{\gamma} | \Phi_{I}^{N} \rangle \approx \langle \phi_{EFlm} | d_{\gamma} | \phi_{FI}^{D} \rangle, \quad (4) $$

where we have neglected a small term involving the conjugated Dyson orbital¹⁴, and denoted $\gamma$ as the dipole operator component. Hence, the dipole transition moments reduce to single-particle transition matrix elements between the Dyson and the continuum states. Consequently, in the correlated single-channel approach, computing the dipole matrix elements requires evaluating the Dyson orbitals $\phi_{FI}^D$, and the continuum orbitals $\phi_{EFlm}$. This calculation is carried out in several steps. The geometry of the electronic ground state minimum of the most abundant conformer of the glycine molecule is optimized using the Gaussian09 package at the DFT/B3LYP level of theory with 6311+g(3df, 2p) basis set. To evaluate the continuum orbital, first a Kohn Sham DFT calculation is performed with the
Amsterdam Density Functional (ADF) package²¹ with the LB94 functional to obtain the ground-state density of glycine. The LB94 functional is chosen because it accurately describes the asymptotic behavior of the exchange potential, with a proper long range coulomb tail²². The calculated ground-state density is initially expressed in terms of Gaussian-type orbitals, which are unsuitable for describing photoionization processes, particularly the rapid oscillations of continuum states¹⁴,²³ To address this, the density is used as an initial guess to construct the Hamiltonian matrix in a new B-spline LCAO basis set, from which both the bound and continuum orbitals are extracted. The new primitive basis functions consist of B-splines radial functions multiplied by spherical harmonics, which describe the angular component¹⁴,²⁰:

$$ \chi_{ilm}(r, \theta, \phi) = \frac{1}{r} B_i(r) Y_{lm}(\theta, \phi). \qquad (5) $$

The LCAO basis set consists of a large one-center expansion (OCE) located at the molecule's center of mass, containing 137 B-splines within a sphere of radius $R_{max}$ = 30 a.u. and a spherical harmonics expansion with maximum angular momentum of $L_{max}$ = 16. The OCE is supplemented by small off-center (OC) expansions positioned over the nonequivalent nuclei, with radii ranging from 0.8 a.u to 1.2 a.u.. The maximum angular momentum is set to $l_{max}$ = 1 for hydrogen atoms and $l_{max}$ = 2 for heavier atoms. The photoelectron continuum states are obtained by solving the generalized eigenvalue problem using the Galerkin method¹⁴.

The Dyson orbitals are computed at the state-averaged multiconfigurational self-consistent field (SA-MCSCF) level using the Open Molcas software package²⁴. A restricted active space (RAS) self consistent field (SCF) is applied to both neutral and cationic glycine, considering up to 11 states. The inactive space consists of the first eight (8) Hartree-Fock (HF) orbitals, while the three active subspaces are defined as follows: the RAS1 includes twelve (12) HF orbitals; RAS2 contains no orbital; RAS3 contains nine (9) valence orbitals. The maximum occupation number of the RAS3 subspace is set to two, implying an equal maximum number of holes in the RAS1 for neutral glycine and a maximum of three holes for the cationic state. The configuration accounts for two-hole-one-particle (2h1p) and three-hole-two-particle (3h2p) excitations, defining a RAS(23,21). The MCSCF orbitals are then projected onto a B-spline LCAO basis for expansion. The validity of this approach has been previously verified in Ref. 12. Finally, the photoionization cross section ($\sigma_i$) is evaluated as a function of the dipole transition moments for the eleven lowest states ($i = 1...11$) of the
glycine cation.

$$ \sigma_{\alpha}(\omega) = \frac{4\pi^2\omega}{\hbar c} \sum_{lm} |\mathbf{D}_{EFIlm\gamma}|^2. \qquad (6) $$

We have recalculated the cross sections to compare them with those of Ref. 12. We observe a good agreement and these cross sections start to be very small above 25 eV. For our purpose, what is important is the ionization probabilities for a specific experimental pulse. The ionization probability for a given channel α can be expressed as a function of the cross-section:

$$ P_{\alpha} = \frac{c}{4\pi^2} \int \frac{1}{\omega} \sigma_{\alpha}(\omega) I(\omega) d\omega, \qquad (7) $$

where I is the intensity of the experimental pulse as a function of the photon energy. The numerical evaluation of the above integral was performed using an interpolation of experimental data. The first 11 values of $P_{\alpha}$ (normalized) are given in Tab. I.

Figure 1 illustrates the cubic spline interpolation of the experimental spectrum, $I(\omega)$, of the XUV attosecond pulse train (APT), which has a duration of 1.5 fs and was used in the experiments of Castrovilli et al.².

<figure><img src="image_2.png" /><figcaption>FIG. 1: Fitting of the experimental spectrum $I(\omega)$. The XUV radiation was obtained by high-order harmonic generation in krypton gas².</figcaption></figure>
<table><caption>TABLE I: Values of the population probabilities, $P_\alpha$, for the first 11 electronic states.</caption><thead><tr><th>P<sub>α</sub></th><th></th></tr></thead><tbody><tr><td>P<sub>1</sub></td><td>0.104</td></tr><tr><td>P<sub>2</sub></td><td>0.111</td></tr><tr><td>P<sub>3</sub></td><td>0.124</td></tr><tr><td>P<sub>4</sub></td><td>0.093</td></tr><tr><td>P<sub>5</sub></td><td>0.120</td></tr><tr><td>P<sub>6</sub></td><td>0.083</td></tr><tr><td>P<sub>7</sub></td><td>0.090</td></tr><tr><td>P<sub>8</sub></td><td>0.078</td></tr><tr><td>P<sub>9</sub></td><td>0.103</td></tr><tr><td>P<sub>10</sub></td><td>0.092</td></tr><tr><td>P<sub>11</sub></td><td>0.000001</td></tr></tbody></table>

## B. The Vibronic Coupling Hamiltonian Model

To investigate the electronic coherence effects, we constructed a vibronic coupling Hamiltonian model and performed quantum dynamics simulations using the Multilayer Multiconfiguration Time-Dependent Hartree (ML-MCTDH) method.

Adiabatic electronic states exhibit discontinuities in their first derivatives with respect to nuclear coordinates due to degeneracies occurring at certain nuclear conformations, such as conical intersections of potential energy surfaces. This divergence in the derivative, known as non-adiabatic coupling, prevents the use of an adiabatic basis for quantum dynamics simulation. To address this issue, a (quasi)-diabatic electronic basis can be introduced through a suitable unitary transformation$^{25,26}$. A common approach to building a vibronic coupling Hamiltonian is to assume a diabatic electronic basis and expand the diabatic potentials in a Taylor series around the equilibrium geometry$^{27-29}$. The resulting diabatic Hamiltonian model is then given by:

$$H = H^{(0)} + W^{(1)} + W^{(2)}, \qquad (8)$$

where the zeroth order term $H^{(0)}$ represents the ground-state Hamiltonian, which consists of the sum of the nuclear kinetic energy and an approximate ground-state harmonic potential
with the ionization energies added at the equilibrium geometry. The matrix elements of $\mathbf{H}^{(0)}$ can then be expressed as:

$$H_{ij}^{(0)} = (T_N + V_N^{\text{harm}}) \delta_{ij} + (E_i - E_0) \delta_{ij}, \qquad (9)$$

$E_0$ represents the ground state energy and $E_i$ denotes the energy of the ith-state at the equilibrium geometry.

$$\mathbf{T}_N + \mathbf{V}_N^{\text{harm}} = \sum_k \frac{\omega_k}{2} \left( -\frac{\partial^2}{\partial Q_k^2} + Q_k^2 \right), \qquad (10)$$

with $\omega_k$ the frequency of the k$^{th}$ normal mode $Q_k$. The matrix element of the first order term $\mathbf{W}^{(1)}$ are given by 30:

$$W_{ij}^{(1)} = \sum_k [\nabla \mathbf{H}]_{ij}^{ij} Q_k. \qquad (11)$$

If the gauge condition is imposed by setting the adiabatic basis identical to the diabatic one at the Franck-Condon point, it can be shown that the gradient of the Hamiltonian is given by:

$$[\nabla \mathbf{H}]_{k}^{ij} = \begin{cases} (\nabla E_i)_{Q_k} & , \text{ if } i = j \\ (E_i(Q_k) - E_j(Q_k)) \langle \varphi_i(Q_k) | (\nabla \varphi_i)_{Q_k} \rangle & , \text{ if } i \neq j. \end{cases} \qquad (12)$$

In this work, we have neglected the off-diagonal matrix elements of the second-order term and write it as follows:

$$W_{ij}^{(2)} = \sum_k \frac{1}{2} [\nabla^2 E_i]_{Q_k} Q_k^2 \delta_{ij}. \qquad (13)$$

The energies are computed using a RASSCF calculation within the Molpro software package. The gradient of the energies and of the wave functions is obtained by the finite difference method. The overlaps between the displaced geometries were evaluated in Molpro.

### C. Quantum Dynamics Simulation

With the Hamiltonian constructed above, our goal is to solve the time-dependent Schrödinger equation (TDSE). This is achieved using the multilayer multiconfiguration time-dependent Hartree (ML-MCTDH) method$^{15–18}$, an extension of the usual MCTDH approach$^{31–35}$. To handle a larger system than those with the standard method, a time-dependent basis known as single particle functions (SPFs) is introduced. These functions
form a smaller basis compared to the primitive basis, giving rise to the MCTDH ansatz for the wave function:

$$ \Psi(q_1, \dots, q_f, t) = \sum_{j_1=1}^{n_1} \dots \sum_{j_f=1}^{n_f} A_{j_1 \dots j_f}^{1}(t) \prod_{\kappa=1}^{f} \varphi_{j_{\kappa}}^{(1;\kappa)}(q_{\kappa}, t). \quad (14) $$

The SPFs are expanded in terms of time-independent primitive basis set:

$$ \varphi_{j_{\kappa}}^{(1;\kappa)}(q_{\kappa}, t) = \sum_{i_{\kappa}=1}^{N_{\kappa}} A_{i_{\kappa}j_{\kappa}}^{(2;\kappa)}(t) \chi_{i_{\kappa}}^{(\kappa)}(q_{\kappa}). \quad (15) $$

This approach reduces the exponential scaling of the wave function and can be viewed as a two layers scheme. The scaling base can be further reduced by applying mode combination, where the physical coordinates $q_1, \dots, q_f$ are grouped into logical coordinates $Q_1^1, \dots, Q_p^1$. Each logical coordinate consists of one or several physical coordinates: $Q_{\kappa}^1 = \{q_{1_{\kappa}}, \dots, q_{d_{\kappa}}\}$. As a result the SPFs become multidimensional, and the wave function ansatz can be expressed as:

$$ \Psi(Q_1^1, \dots, Q_p^1, t) = \sum_{j_1=1}^{n_1} \dots \sum_{j_p=1}^{n_p} A_{1; j_1 \dots j_p}^{1}(t) \prod_{\kappa=1}^{p} \varphi_{j_{\kappa}}^{(1;\kappa)}(Q_{\kappa}^1, t). \quad (16) $$

An MCTDH-like ansatz can then be used to expand the multidimensional time-dependent SPFs $\{\varphi_{j_{\kappa}}^{(1;\kappa)}(Q_{\kappa}, t)\}$, effectively introducing a third layer of expansion coefficients. However, combined mode SPFs are often too large to be efficiently propagated, so they are further divided into smaller groups of logical coordinates. This introduces additional layers of coefficients with manageable sizes. The hierarchical expansion is ultimately truncated by expressing the SPFs of the final layer as a direct product expansion in terms of time-independent primitive basis functions, similar to eqn. (15). This procedure can be repeated as many times as necessary; forming the ML-MCTDH method. The multilayer expansion of the wave function can be visualized using ML-tree diagrams, as described in the supplementary materials. According to recent work by David Mendive-Tapia et al.³⁶⁻³⁸, this type of hierarchical tree is the most efficient for vibronic coupling Hamiltonians, provided that the modes are arranged according to their correlations.

We employed a Hermite discrete variable representation (DVR) for the primitive basis functions and included up to $N = 11$ electronic states. All the 24 internal degrees of freedom have been included using the normal coordinates of the electronic ground state. The electronic degree of freedom, representing the electronic states are separated from the
vibrational degrees of freedom by the first layer of the ML-MCTDH hierarchy. The ML-MCTDH equations of motion were integrated using the variable mean-field technique with a fifth-order Runge-Kutta integrator and maintaining an accuracy of $10^{-7}$. The initial wave function was constructed as product of single particle Gaussian functions for each mode. All ML-MCTDH calculations were performed using the Heidelberg MCTDH package$^{39}$. Convergence with respect to the number of SPFs and primitive basis functions was verified by monitoring the population of the natural orbitals. After 15fs of propagation, the highest population of the least occupied natural orbital was approximately $4 \cdot 10^{-3}$. As initial states, we have taken the weights corresponding to the ionization probabilities extracted from the experimental pulse in Table I. Since the relative phases of the created electronic wave packets are unknown, they were chosen randomly. Several sets of random phases were tested, and the differences in the evolution of the purity (see below) were found to be negligible.

A way to quantify coherence in a quantum system is by evaluating the purity, defined as the trace of the square of the density matrix$^{40}$:

$$ \mathcal{P} = \mathrm{Tr}[\rho^2]. \qquad (17) $$

We plotted the purity as a function of time and verified its convergence with respect to the number of SPFs and primitive basis functions. Convergence is considered achieved if the observable remains unchanged when increasing the number of SPFs and primitive basis functions. We examined both the impact of including the second-order term in the Hamiltonian—by comparing the linear and quadratic vibronic coupling Hamiltonians—and the effect of the initial phase between electronic states on purity. The latter was found to be negligible.

## D. Surface Hopping

The long-time dynamics induced by the APT is simulated using a semi-classical method, specifically the trajectory surface hopping. In this approach, the electronic motion is described by the quantum mechanical Schrödinger equation, while the nuclear motion by the classical Newton’s equations. The electronic wave function $|\Psi[\vec{R}(t); \vec{r}, t]\rangle$, which depends parametrically on the nuclear positions $\vec{R}$ at each time $t$, is the solution of the time dependent
<figure><img src="image_3.png" /><figcaption>FIG. 2: Purity as a function of time.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>FIG. 3: Root means square displacement of atoms as a function of time.</figcaption></figure>

dent equation:

$$i\hbar \frac{\partial |\Psi[\vec{R}(t); \vec{r}, t]\rangle}{\partial t} = \hat{H}_{eff}[\vec{R}(t); \vec{r}] |\Psi[\vec{R}(t); \vec{r}, t]\rangle, \quad (18)$$

where the effective Hamiltonian $H_{eff}$ includes nucleus-nucleus, electron-nucleus, electron-electron interactions, and the electronic kinetic energy. To solve the time-dependent Schrödinger equation (18), the electronic wave function is expanded in basis functions $|\phi_{\alpha}[\vec{R}(t); \vec{r}]\rangle$, which are taken as eigenfunctions of the effective Hamiltonian (adiabatic states):

$$|\Psi[\vec{R}(t), \vec{r}, t]\rangle = \sum_{\alpha} c_{\alpha}(t) |\phi_{\alpha}[\vec{R}(t); \vec{r}]\rangle. \quad (19)$$

This leads to the equation for the time evolution of the coefficients:

$$\frac{\partial c_{\beta}(t)}{\partial t} = -\sum_{\alpha} \left\{ \frac{i}{\hbar} H_{\beta\alpha}[\vec{R}(t)] + K_{\beta\alpha}[\vec{R}(t)] \right\} c_{\alpha}(t), \quad (20)$$

where $H_{\beta\alpha}$ represents the adiabatic Hamiltonian, and the second term $K_{\beta\alpha}$ is the non-adiabatic coupling. The equation (20) is then solved using the Runge-Kutta method.

To account for the initial geometry spreading and velocity of the nuclear wave packet in the ground vibrational state, a set of nuclear trajectories, defined by the position $\vec{R}(t)$ and the velocity $\vec{v}(t)$ of the nuclei, are considered. Each nuclear trajectory is propagated using the well-known velocity Verlet algorithm, evolving on the adiabatic potential energy surfaces.

The non-adiabatic transitions at avoided crossings of the nuclear quantum wave packet are approximated by instantaneous hops between adiabatic states, utilizing the Tully's
fewest switches algorithm$^{41}$. The well-established coherence problem in surface hopping algorithm$^{42,43}$ is addressed by the scheme proposed by Truhlar and coworkers$^{44,45}$. Energy conservation is maintained by rescaling the nuclear velocity after each successful hop. To characterize the vibrational ground state of the nuclear wave packet, a set of 100 initial geometries $\vec{R}(0)$ and corresponding velocities $\vec{v}(0)$ are obtained by sampling the ground state phase space of the glycine molecule using the harmonic oscillator Wigner quasiprobability distribution$^{46}$.

For the independent TSH calculations, the initial coefficients $c_0$ are all set to one. Regarding the electronic structure method, we performed restricted active space (RAS) self consistent field (SCF) calculations for both neutral and cationic glycine, as implemented in the Firefly package$^{47,48}$. The inactive and active subspaces are identical to those used in the first part to compute the Dyson orbitals. To accurately calculate the energy gradient with respect to the nuclear coordinates, we employed an orbital optimization procedure and performed a state-average calculation that includes up to eleven states.

For each nuclear configuration, one trajectory is run for each electronic state $\alpha$. We have included up to twenty (20) electronic states in the superposition of states created by the pump in the cation for the RAS, with the first ten states selected to run the nuclear trajectories.

### III. RESULTS AND DISCUSSION

#### A. Electronic coherences

Figure 2 presents the purity as a function of time for the quadratic model, while Figure 3 displays the corresponding root mean square displacement (RMSD) of atoms with respect to their equilibrium positions. Within 5 fs, the purity drops below 0.2, approaching its minimum value, as it must satisfy the relation: $\frac{1}{N} \le P \le 1$, where N = 11 in our case. Meanwhile, the nuclei remain nearly stationary, with a maximum RMSD of just 0.7 Å. This result demonstrates that decoherence occurs before significant nuclear motion sets in. This rapid decoherence justifies treating the electronic states independently at the start of the dynamics. A full quantum treatment of the subsequent coupled electron-nuclear dynamics is very difficult due to the high number of degrees of freedom—particularly the construction of diabatic potentials for large-amplitude nuclear motions. For the long-time dynamics, we
thus adopt a mixed quantum-classical approach.

Previous works found that some electronic coherences might persist after sudden ionizing
a specific orbital<sup>11,49</sup> and after ionizing with a XUV broadband attosecond pulse<sup>12</sup>. However,
it is critical to emphasize that these works either do not account for decoherence mecha-
nisms or treat them in a highly approximate manner—limiting their ability to predict the
longevity of the observed dynamics. The former<sup>11,49</sup> are molecular-driven (through electron
correlation and relaxation) and the latter are, at least initially, light-induced due to the
≈ 18eV bandwidth encompassing many excited electronic states. Kuleff and Cederbaum<sup>11</sup>
have simulated the electronic dynamics by employing the direct propagation of the multi-
electron wave packet of the ionized system for a fixed nuclear frame. Depending on the
ionized orbitals, they observed ultrafast charge migration with cycles ranging from 7 to 20
fs, much longer charge migration times or no migration at all. In ref 12, the broadband
laser induces a coherent superposition of cation states most of them containing a delocalised
hole. The dynamics of the initial hole density and the charge fluctuations at different atomic
centers of the molecule were calculated using a time dependent DFT-Ehrenfest dynamics
(TDDFT-ED) method to account for the electron-electron correlations, nuclear motion and
nonadiabatic effects. The initially created coherences appear to be only slightly affected
in the first 10 fs when including the nuclear motion, the effect becoming more visible after
15 fs. Nevertheless, in these calculations the initial geometry distribution corresponding to
the initial nuclear wavepacket is ignored. The inclusion of the distribution may prevent the
observation of the coherences. Like in ref.<sup>12</sup>, we have used a broadband ultrashort laser
pulse that coherently excites a large manifold of states. All the cationic excited states (up
to 20 eV) considered in the above studies, were included in the vibronic coupling Hamil-
tonian model as the attosecond pulse train used in our calculations encompasses all these
states. The vibronic coupling model has two main limitations: it remains accurate only for
short times and it cannot describe correctly motions of large amplitude as those involved
in a H<sup>+</sup>/H transfer. However, since decoherence occurs very rapidly and we observe no
difference between the linear and quadratic models (with the quadratic model being valid
for longer times than the linear one), we are confident that this description of decoherence
is accurate. The loss of the coherence is so fast that it is unlikely that it can influence the
transfer..

In ref⁵⁰, the broadband laser induces a coherent superposition of cation states most of them
containing a delocalised hole. The dynamics of the initial hole density and the charge fluctuations at different atomic centers of the molecule were calculated using a time dependent DFT-Ehrenfest dynamics (TDDFT-ED) method to account for the electron-electron correlations, nuclear motion and nonadiabatic effects. The initially created coherences appear to be only slightly affected in the first 10 fs when including the nuclear motion, the effect becoming more visible after 15 fs. Nevertheless, in these calculations the initial geometry distribution corresponding to the initial nuclear wavepacket is ignored. The inclusion of the distribution may prevent the observation of the coherences. Like in ref.⁵⁰, we have used a broadband ultrashort laser pulse that coherently excites a large manifold of states. All the cationic excited states (up to 20 eV) considered in the above studies, were included in the vibronic coupling Hamiltonian model as the attosecond pulse train used in our calculations encompasses all these states. The vibronic coupling model has two main limitations: it remains accurate only for short times and it cannot describe correctly motions of large amplitude as those involved in a H⁺/H transfer. However, since decoherence occurs very rapidly and we observe no difference between the linear and quadratic models (with the quadratic model being valid for longer times than the linear one), we are confident that this description of decoherence is accurate. The loss of the coherence is so fast that it is unlikely that it can influence the transfer dynamics.

To conclude this section, we mention the work of Ref. 51, which also studied the stability of quantum coherence. However, their study included only two electronic states so that their conclusions are not applicable to the present case.

## B. Initial conditions and migration pathways

Figures 4 and 5 show the equilibrium geometries of the glycine conformers used in this study and that were optimized at the DFT-B3LYP/(3df,2p) level of the theory. In both cases, the heavy atoms lie in a plane and NH₂ group exhibits a bifurcated orientation and hydrogen bonds with the carbonylic or hydroxylic oxygen. In turn, the lowest lying cation conformer has a geometry where the N-H bonds are in the molecular plane⁶. This change requires rotation about the C-N bond and leads to a significant shortening of the H-bond from 2.859 Å (at BHLYP level)⁵ to 2.099 Å (at B3LYP level)⁶. Such a geometry and strengthening of the H-bond should favour the H⁺/H transfer from N to O. The adiabatic ionization
potential of the cation ground state is calculated at 9.0 eV (at CCSD(T) and BHLYP)⁵ and determined experimentally at 9.09 eV. But, vertical ionization in the Franck-Condon region requires a higher energy as determined experimentally by photoelectron spectroscopy at 10 eV⁵² and calculated at DFT-B3LYP/6-311++G(d,p) level of theory at 9.82 eV⁷. In the Table VII of SI we give the vertical ionization potentials (IPs) calculated in this work. The lowest ionization potential stands at 9.95 eV and was calculated using CCSD(T) method with the triple-zeta correlation-consistent basis set cc-pVTZ in the Gaussian software package. To model the distribution of the initial nuclear wavepacket we considered a temperature of T = 0 K and a Wigner distribution consisting of 100 geometries. Since we assume a statistical distribution in view of the fast decoherence time, for each nuclear configuration we launched a trajectory for each of the 10 lowest electronic cationic states. A total of 1010 trajectories have been launched for the lowest lying conformer (conformer 1) and 330 trajectories for the other conformer (conformer 2). The analysis of the interatomic distances between the hydrogen atoms and the nitrogen or the carbon atoms to which they are bonded reveals cases where a H atom breaks away to form a new bond with another atom. These intra-molecular H migrations, observed on 17 trajectories for conformer 1 and on 2 trajectories for the conformer 2, follow two main pathways: a migration from the amino group to the oxygen of the carbonyl group (pathway 1) and a migration from alpha carbon atom $C_α$ to the nitrogen atom (pathway 2), in agreement with the pathways reported previously⁴,⁶,⁷. A rare situation occurs in the case of conformer 1 where the hydrogen atom migrates from $C_α$ to the second carbon, resulting in a release of a hydroxyl group (OH). For the other conformer, the migration of H from the nitrogen to the oxygen of the carboxyl group is followed by a release of a water $H_2O$ molecule. The results are reported in Tables I and II of the SI and in Figs.6-11 in the main article and Figs. 4-33 in SI, for conformer 1 and in Figs.34-37 in SI, for conformer 2.

### C. Trajectories analysis

In the absence of migration, hydrogen atoms oscillate around their equilibrium positions at the initial site. However, as shown in Figures 7 and 10, at a specific time, The interatomic distance between one of the hydrogen (H) atoms and the initial site (nitrogen or carbon alpha, $C_α$) increases, while its distance to the final site (either an oxygen or a nitrogen
atom) decreases. As this change becomes more pronounced, the H atom migrates toward the final site. Figures 8 and 11 illustrate this process through successive snapshots of the geometry modifications during H migration from nitrogen (N) to oxygen (O) and from $C_α$ to N, respectively.

Migration is said to occur when the H atom is closer to the final site than to the initial site. Among trajectories exhibiting migration along Pathway 1, three cases show the H atom returning to the initial site (the amino group). Dynamics movies confirm that the migrating H atom forms a bond with the carbonyl oxygen in Pathway 1 and with the nitrogen atom in Pathway 2.

Upon analyzing all trajectories with migration, three key observations emerge:

(a) Migration times vary significantly across trajectories, ranging from 37 fs to 96 fs (Pathway 1) and 26 fs to 91 fs (Pathway 2). On average, migration in Pathway 1 takes longer than in Pathway 2.

(b) Only a narrow range of the 100 initial geometries—those closely resembling the equilibrium structure—lead to migration.

(c) Most critically, for all trajectories except three, H migration occurs on the ground-state potential energy surface (PES) of the cation, to which the molecule relaxes prior to migration. This is evident from successive hops between PESs.

The first observation aligns with the longer initial distances between the H atom and the final migration sites in Pathway 1 (initially 2.8 Å, corresponding to the H bond length) compared to Pathway 2. However, the PES landscape after ionization drives the dynamics and determines both the probability and timing of migration. For instance, in Pathway 1, two trajectories starting from the same geometry but on different electronic states (the 5th and 9th, separated by ~3 eV) exhibit only a modest 20 % difference in migration time due to energetics and nonadiabatic couplings. Conversely, for four trajectories launched on the same electronic state (the 6th state) but with different initial geometries and internal energies, the fastest and slowest migration times differ by a factor of three. Similar trends are observed for trajectories on Pathway 2 (Table II in SI). In Pathway 1, while four out of nine geometries leading to migration have dihedral angles similar to the neutral equilibrium structure, the migrating H atom typically has either the longest bond or the smallest dihedral
angle (H-C-C-O). The fastest transfer occurs when the dihedral angle is smallest. For the two trajectories launched from the same geometry but different electronic states, the same H atom migrates, further emphasizing the importance of initial geometry. These observations suggest that initial geometry and internal energy significantly influence migration time.

The second observation (b) is surprising given the structural differences between the neutral and cation ground states equilibrium geometries. However, geometries leading to migration often feature a more planar orientation or a shorter distance between the migrating H and O (or a more favorable arrangement). A rotation about the C-N bond transforms the neutral ground-state geometry into one resembling the cation's, with N-H bonds in the molecular plane. Such a geometry strengthens the O-H bond, facilitating migration via Pathway 1. Analysis of the O-H bond length reveals that it decreases well before the N-H bond increases, likely due to rotation about the C-N axis. In some cases, the O-H bond begins to increase before complete relaxation to the ground state, suggesting that torsional vibrations may drive this change. Snapshots from Figures 8 indicate that a more planar arrangement is achieved before migration. The NH₂ twist and CNH torsion motions were not explicitly considered in our simulations (performed at T=0 K). At higher temperatures, initial vibrational excitation could accelerate migration, particularly if the system starts in a favorable geometry where the NH₂ group lies in the COOH plane, with one H atom pointing toward the O atom via a hydrogen bond. Excitation of specific vibrational modes, such as NH₂ twist and CNH torsion (with frequencies of ~ 280 cm⁻¹ and ~ 230 cm⁻¹, respectively), could also induce such geometries. These modes are likely excited under experimental conditions (T=443–463 K) used in ref. 2. Simulations at realistic temperatures would enable exploration of these effects on dynamics and migration.

Regarding the third observation (c), we have computed the minimum energy path (MEP) of the isomerization reaction using the nudge elastic band method⁵³ as implemented in ORCA software package⁵⁴,⁵⁵ to explain why hydrogen transfer occurs when the active state reaches the cationic ground state. The optimized geometries along the MEP are used to compute excited-state energies at the SA-RASSCF level, providing an estimate of the energy barriers in excited states. Detailed calculations are provided in the Supporting Information. The results show that proton transfer reaction barriers are higher in excited states than in the ground state. This confirms that hydrogen migration is more likely to occur in the cationic ground state. However, these results must be interpreted with caution because the
simulations were performed at 0 K, and higher temperatures could yield different outcomes.

Migration was observed on the first excited state in three trajectories (Traj 73-6, Traj 79-6, and Traj 56-10), all of which exhibit minimal relaxation until just before migration. Interestingly, these trajectories show shorter migration times despite slower relaxation. In Pathway 1, seven out of nine trajectories start from cation states in the energy range 13.5–17.5 eV, with four from the same state (16.57 eV). For Pathway 2, all but one trajectory start from higher-energy states (16.5–18.8 eV), with half above 18 eV.

Pathway 1 was experimentally observed after photoionization at 60 eV using narrow-band synchrotron radiation, detected through the fragment ion HCOOH+, which cannot exist without prior hydrogen transfer to the carbonyl O. PEPICO measurements identified the energy range for this process as 12.25–20 eV, peaking around 13.75 eV⁵⁶. This matches the migration region observed in simulations and corresponds to the excitation of several closely spaced 1h configuration cation states involving ionization of the 3a", 2a", 1a" and 14a', 13a', 12a', and 11a' orbitals. Sudden ionization of these orbitals leads to hole migration over different timescales due to electron correlations¹¹,⁴⁹. For example, ejecting an electron from orbital 14a' creates a hole spread across the molecule, with highest density on the NH₂ group¹¹. Within ~ 3.5 fs, the hole transfers to orbital 13a', localized mainly on the COOH end, and returns to the initial state in ~ 7 fs.

<figure><img src="image_5.png" /><figcaption>FIG. 4: Most abundant conformer. Global minimum, *E* = 0 eV. Its energy point zero is equal to 2.168 eV.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>FIG. 5: Second most abundant conformer. Δ*E* = 0.07 eV. The geometries are optimized at the DFT-B3LYP/(3df,2p) level.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>FIG. 6: Energy diagram illustrating the system dynamics for a trajectory “TRAJ2-5” initiated in the 5th excited state. The solid red line represents the total energy of the system, while the black curve corresponds to the active state, showcasing the sequential transitions (hops) from the 5th excited state down to the ground state.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>FIG. 7: Variation of the N-H and O-H distances: Case of a trajectory starting in the 5th excited state. The plot highlights the gradual decrease in the $O_{10}-H_5$ distance and an increase in the $N-H_5$, illustrating the migration of the hydrogen atom within 72 fs, consistent with the first migration pathway.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>FIG. 8: Snapshots of the geometries during the dynamics showing the migration path 1.</figcaption></figure>

## D. Migration Probability

We define the probability for a hydrogen to migrate as:

$$P_m(t) = \frac{\sum_{\alpha,k} P_{\alpha}^k \Theta_{\alpha}^k}{\sum_{\alpha,k} P_{\alpha}^k}, \qquad (21)$$

where the sum is taken over the states $\alpha$ and the geometries $k$, $P_{\alpha}^k$ is the ionization probability defined before in equation (7). For computational simplicity, we have used the probabilities of the equilibrium geometry as the ionization probabilities for all the geometries obtained
<figure><img src="image_10.png" /><figcaption>FIG. 9: Energy diagram during the dynamics: Case of a trajectory ‘TRAJ14-7” starting in the 7th excited state. The solid red line represents the total energy of the system, while the black curve corresponds to the active state, illustrating the successive hops from the 7th state down to the ground state.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>FIG. 10: Variation of the C-H and N-H distances: Case of a trajectory starting in the 7th excited state. The figure shows a decrease in the N-H distance and a simultaneous increase in the C1-H3), indicating the migration of the hydrogen atom in 31 fs, consistent with the second migration pathway.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>FIG. 11: Snapshots of the molecular geometries during the dynamics, illustrating the migration along Path 2.</figcaption></figure>

by Wigner sampling. $\Theta_{\alpha}^k$ is given by:

$$ \Theta_{\alpha}^{k} = \begin{cases} 1 & \text{if } d_{N-H_5} < d_{H_5-O_{10}}, \\ 0 & \text{otherwise.} \end{cases} \qquad (22) $$

For each trajectory exhibiting hydrogen migration, we calculate its probability at the moment of migration. Figure 12 illustrates the migration probabilities for each individual pathway, while Figure 13 displays the total intra-molecular hydrogen migration probability,
which is the sum of the probabilities for both pathways.

It is evident that, in general, hydrogen migration in the most abundant conformer of the glycine cation is more likely to occur via the second pathway. The maximum migration probability for this pathway is approximately 6×10⁻³, compared to 4×10⁻³ for the first pathway. Regardless of the pathway, the highest observed probability for hydrogen migration is 10⁻².

In contrast, Ref. 9 reports significantly lower migration probabilities for pathway 1, specifically 8.0×10⁻⁶ % and 1.8×10⁻⁵ % for conformers 1 and 2, respectively. These discrepancies may arise from differences in the calculation models used, as well as the much lower excess energy considered in those studies (on the order of 5 eV).

Experimentally, the migration probability for pathway 1 was estimated from the mass spectrum obtained following ionization with the APT method. This was determined by the relative intensity of the HCOOH⁺ peak compared to the COOH⁺ peak, yielding a value of 9.4 %². However, this value corresponds only to specific migration channels detectable through these fragmentation pathways. To compare with the theoretical values presented here, all ionization and dissociative ionization channels indicative of migration would need to be summed and their relative intensity compared to the total mass spectrum signal.

Since not all such channels could be identified, a direct comparison is not feasible. Nevertheless, based on the magnitude of the HCOOH⁺ peak, it can be inferred that the experimental migration probability is likely smaller than the theoretical value of 0.4 % reported here.

In the graphs presented in Figures 12 and 13, each data point corresponds, in the experimental context, to the number of molecules that exhibited migration and were detected at a specific moment in time. However, unlike the experimental measurements, the graph in Figure 12 lacks sufficient statistical robustness. Despite this limitation, and by drawing an analogy between the two, we can still compare the time dependence of Figure 12 (Pathway 1) with the red curve in Figure 2 of reference 2. Two key differences emerge from this comparison:

* In the experiment, migration appears to initiate at very short timescales, following a delay of just a few femtoseconds (as seen in Figure 2 of Reference 2).

* In the simulations, the migration probability seems to reach a plateau more rapidly.
This observation aligns with the trajectory analysis, which indicates that migration occurs quickly after significant nonadiabatic relaxation to the lowest electronic states.”

Given that proton migration is a rare event, as also observed experimentally, a significantly larger number of simulations would be required to determine whether shorter migration times could be detected. Furthermore, it is plausible that an even more extensive dataset of trajectories—far exceeding those analyzed in this study—might uncover a small number of exceptionally fast migration events, potentially occurring in excited electronic states. Additionally, caution must be exercised when comparing our results with the experimental data from Ref. 3, as several approximations were employed in the simulations. For instance, the simulation does not account for the experimental temperature, which limits the range of initial conditions and geometric configurations considered. These factors could influence the outcomes and should be taken into account when interpreting the results.

<figure><img src="image_13.png" /><figcaption>FIG. 12: Migration probability as a function of time for the two main paths of hydrogen migration.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>FIG. 13: Total hydrogen migration probability as a function of time in the photo-ionized glycine.</figcaption></figure>

E. Charge Analysis

We performed a charge analysis to determine whether it is an hydrogen (H) or a proton (H⁺) that migrates. We evaluated the charge densities for each state by performing the same RASSCF calculation as before, but using the Molpro Software package³⁰,⁵⁷. We tested the RASSCF in Molpro and obtained the same energy as that obtained in Firefly for the
trajectory surface hopping. The average charge density $\rho_m$ is then obtained by summing the densities of each state, weighted by the corresponding population $p_i$ as in the following equation:

$$ \rho_m(\mathbf{r}) = \sum_{i}^{\text{states}} p_i \rho_i(\mathbf{r}), \qquad (23) $$

The mean charge density is used to evaluate the net atomic charges by the Density Derived Electrostatic and Chemical (DDEC6) approach$^{58}$. For each migration path, we calculated the net charge on the migrating hydrogen atom as a function of time (each 10 fs). The results are displayed in Figures 14 and 15 for the first migration pathway, for different starting geometries and electronic states. One can first observe a decrease of the net charge on the hydrogen, followed by a significant jump close to 1 (0.89 for Figure 14 and 0.8 for Figure 15) when the hydrogen changes its electronic environment and migrates. Similar results are shown in Figures 16 and 17 for the second migration pathway, where the same behavior of the net charge of the migrating hydrogen is observed. Due to the rapid increase of the charge, approaching 1, during the hydrogen transfer, one can conclude that, the process is a proton transfer, regardless of the starting configuration, starting electronic state, or the migration pathway.

<figure><img src="image_15.png" /><figcaption>FIG. 14: Net atomic charge on the migrating hydrogen during the dynamics: case of path 1, 5th starting state.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>FIG. 15: Net atomic charge on the migrating hydrogen during the dynamics: case of path 1, 9th starting state.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>FIG. 16: Net atomic charge on the migrating hydrogen during the dynamics: case of path 2, 7th starting state.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>FIG. 17: Net atomic charge on the migrating hydrogen during the dynamics: case of path 2, 4th starting state.</figcaption></figure>

IV. CONCLUSION

In Ref. 11, Kuelff and Cederbaum investigated hole charge migration following the ionization of three glycine conformers in which the nuclei were not allowed to move and predicted that electronic coherences can survive for a long time. In this work, we performed full-dimensional quantum simulations, which revealed that, for the pulses used in the experiment reported in Ref. 2, nuclear motion induces extremely rapid electronic decoherence (on the order of less than 3 fs). While our model cannot accurately describe large-amplitude motions like proton transfer, the decoherence is so fast that it does not influence the proton transfer process. We employed a surface-hopping approach to simulate non-adiabatic dynamics, starting from the first 11 electronic states. Proton transfer was observed in some trajectories for both conformers. Notably, proton migration consistently occurred when the system reached the electronic ground state of the cation, accompanied by a significant release of vibrational energy. The initial geometries leading to migration were similar, consistent with our use of a Wigner distribution centered around the equilibrium geometry. This behavior appears to be statistically driven.

When comparing our results to the experimental data in Ref. 2, care must be taken due to the approximations made in our simulations. However, some similarities are evident:
a small fraction of molecules undergo proton migration, and two distinct pathways are observed. Overall, the migration process appears slower in our simulations compared to the experiment, where migration seems to initiate at very short timescales.

Given that proton transfer is experimentally observed to be a rare event, increasing the trajectory count would generate more transfer events and could yield shorter migration times. More importantly, incorporating temperature effects would better replicate experimental conditions while potentially reducing migration times. However, it should be noted that elevated temperatures tend to disrupt quantum coherence and may impede quantum control of the process.

However, it is worth noting that higher temperatures tend to dissipate quantum coherence and hinder quantum control of the process.
V. SUPPORTING MATERIAL

The supporting material includes the description of the "tree" in the ML-MCTDH cal-
culations, the comparison of the purity and root mean square displacement of atoms for the
linear and quadratic vibronic coupling Hamiltonian models, the list of geometries leading to
the migration paths, the description of the Wigner distribution, and the ionization energies
of the first channels of the glycine molecule.

VI. ACKNOWLEDGEMENTS

We acknowledge financial support from the French National Research Agency (ANR) for
the AttoChemistry project (ANR-21-CE29-0005). The authors also express their gratitude
to the CNRS International Research Network (IRN) 'MCTDH' for its financial support. KK
would like to acknowledge the receipt of the short-term scientific mission grants from the
Attochem COST Action (CA18222) which enabled his stays at the Autonomous University
of Madrid and at the University of Trieste. JGV, AP and FM acknowledge the Minis-
terio de Ciencia e Innovación projects PID2022-138288NB-C31 and PID2022-138288NB-
C32 (MCIN/AEI/ 10.13039/501100011033/FEDER UE), the "Severo Ochoa" Programme
for Centres of Excellence in R&D (CEX2020-001039-S), and the "María de Maeztu" Pro-
gramme for Units of Excellence in R&D (CEX2018-000805-M). We thank the Mare Nostrum
Supercomputer of the Red Española de Supercomputación (BSC-RES) and the Centro de
Computación Científica de la Universidad Autónoma de Madrid (CCC-UAM) for providing
computer time.

REFERENCES

<footnote n="1">¹Ahmed H Zewail. Femtochemistry: Atomic-scale dynamics of the chemical bond. *The Journal of Physical Chemistry A*, 104(24):5660–5694, 2000.</footnote>

<footnote n="2">²MC Castrovilli, A Trabattoni, P Bolognesi, P O Keeffe, L Avaldi, M Nisoli, F Calegari, and R Cireasa. Ultrafast hydrogen migration in photoionized glycine. *The journal of physical chemistry letters*, 9(20):6012–6016, 2018.</footnote>
<footnote n="3">³Francesca Calegari, D Ayuso, Andrea Trabattoni, Louise Belshaw, Simone De Camillis, S Anumula, F Frassetto, L Poletto, A Palacios, P Decleva, et al. Ultrafast electron dynamics in phenylalanine initiated by attosecond pulses. *Science*, 346(6207):336–339, 2014.</footnote>

<footnote n="4">⁴Gisbert Depke, Nikolaus Heinrich, and Helmut Schwarz. On the gas phase chemistry of ionized glycine and its enol. a combined experimental and ab initio molecular orbital study. *International journal of mass spectrometry and ion processes*, 62(1):99–117, 1984.</footnote>

<footnote n="5">⁵L Rodriguez-Santiago, M Sodupe, A Oliva, and J Bertran. Intramolecular proton transfer in glycine radical cation. *The Journal of Physical Chemistry A*, 104(6):1256–1261, 2000.</footnote>

<footnote n="6">⁶Silvia Simon, Mariona Sodupe, and Juan Bertran. Isomerization versus fragmentation of glycine radical cation in gas phase. *The Journal of Physical Chemistry A*, 106(23):5697–5702, 2002.</footnote>

<footnote n="7">⁷Sylvain Maclot, Dariusz Grzegorz Piekarski, Rudy Delaunay, Alicja Domaracka, Alain Méry, Violaine Vizcaino, Jean-Yves Chesnel, Fernando Martín, Manuel Alcamí, Bernd A Huber, et al. Stability of the glycine cation in the gas phase after interaction with multiply charged ions. *The European Physical Journal D*, 68:1–7, 2014.</footnote>

<footnote n="8">⁸Dorit Shemesh, Galina M Chaban, and R Benny Gerber. Photoionization dynamics of glycine: The first 10 picoseconds. *The Journal of Physical Chemistry A*, 108(52):11477–11484, 2004.</footnote>

<footnote n="9">⁹D Shemesh and RB Gerber. Different chemical dynamics for different conformers of biological molecules: Photoionization of glycine. *The Journal of chemical physics*, 122(24), 2005.</footnote>

<footnote n="10">¹⁰Sylvain Maclot, Dariusz Grzegorz Piekarski, Alicja Domaracka, Alain Meéry, Violaine Vizcaino, Lamri Adoui, Fernando Martín, Manuel Alcamí, Bernd A Huber, Patrick Rousseau, et al. Dynamics of glycine dications in the gas phase: ultrafast intramolecular hydrogen migration versus coulomb repulsion. *The Journal of Physical Chemistry Letters*, 4(22):3903–3909, 2013.</footnote>

<footnote n="11">¹¹A I Kuleff and L S Cederbaum. Charge migration in different conformers of glycine: The role of nuclear geometry. *Chemical Physics*, 338:320–328, 2007.</footnote>

<footnote n="12">¹²Jorge Delgado, Manuel Lara-Astiaso, Jesús González-Vázquez, Piero Decleva, Alicia Palacios, and Fernando Martín. Molecular fragmentation as a way to reveal early electron dynamics induced by attosecond pulses. *Faraday Discussions*, 228:349–377, 2021.</footnote>
<footnote n="13">¹³Lara-Astias Manuel, Ayuso David, Tavernelli Ivani, Piero Decleva, Alicia Palacios, and Fernando Martín. Decoherence, control and attosecond probing of xuv-induced charge migration in biomolecules. a theoretical outlook. *Faraday Discussions*, 194:41–59, 2016.</footnote>

<footnote n="14">¹⁴Piero Decleva, Mauro Stener, and Daniele Toffoli. Continuum electronic states: The tiresia code. *Molecules*, 27(6):2026, 2022.</footnote>

<footnote n="15">¹⁵H. Wang and M. Thoss. Multilayer formulation of the multiconfiguration time-dependent Hartree theory. *J. Chem. Phys.*, 119:1289–1299, 2003.</footnote>

<footnote n="16">¹⁶H. Wang and M. Thoss. Theoretical study of ultrafast photoinduced electron transfer processes in mixed-valence systems. *J. Phys. Chem. A*, 107:2126–2136, 2003.</footnote>

<footnote n="17">¹⁷Uwe Manthe. A multilayer multiconfigurational time-dependent Hartree approach for quantum dynamics on general potential energy surfaces. *J. Chem. Phys.*, 128:164116, 2008.</footnote>

<footnote n="18">¹⁸Oriol Vendrell and Hans-Dieter Meyer. Multilayer multiconfiguration time-dependent hartree method: Implementation and applications to a henon–heiles hamiltonian and to pyrazine. *The Journal of Chemical Physics*, 134(4), 2011.</footnote>

<footnote n="19">¹⁹Alicia Palacios and Fernando Martín. The quantum chemistry of attosecond molecular science. *Wiley Interdisciplinary Reviews: Computational Molecular Science*, 10(1):e1430, 2020.</footnote>

<footnote n="20">²⁰Mauro Nisoli, Piero Decleva, Francesca Calegari, Alicia Palacios, and Fernando Martín. Attosecond electron dynamics in molecules. *Chemical reviews*, 117(16):10760–10825, 2017.</footnote>

<footnote n="21">²¹G. te Velde, F. M. Bickelhaupt, E. J. Baerends, C. Fonseca Guerra, S. J. A. van Gisbergen, J. G. Snijders, and T. Ziegler. Chemistry with adf. *J. Comput. Chem.*, 22(9):931–967, 2001.</footnote>

<footnote n="22">²²R Van Leeuwen and EJ Baerends. Exchange-correlation potential with correct asymptotic behavior. *Physical Review A*, 49(4):2421, 1994.</footnote>

<footnote n="23">²³Daniele Toffoli, Mauro Stener, Giovanna Fronzoni, and Pietro Decleva. Convergence of the multicenter b-spline dft approach for the continuum. *Chemical physics*, 276(1):25–43, 2002.</footnote>

<footnote n="24">²⁴Francesco Aquilante, Jochen Autschbach, Alberto Baiardi, Stefano Battaglia, Veniamin A Borin, Liviu F Chibotaru, Irene Conti, Luca De Vico, Mickaël Delcey, Nicolas Ferré, et al. Modern quantum chemistry with [open] molcas. *The Journal of chemical physics*, 152(21), 2020.</footnote>
<footnote n="25">²⁵H. Köppel. Diabatic representation: methods for the construction of diabatic states. In W. Domcke, D. R. Yarkony, and H. Köppel, editors, *Conical Intersections*, pages 175–204, Singapore, 2004. World Scientific Co.</footnote>

<footnote n="26">²⁶David Yarkony, Wolfgang Domcke, and Horst Köppel. *Conical Intersections: Electronic Structure, Dynamics & Spectroscopy*. World Scientific, 2004.</footnote>

<footnote n="27">²⁷H. Köppel, W. Domcke, and L. S. Cederbaum. Multimode molecular dynamics beyond the Born-Oppenheimer approximation. *Adv. Chem. Phys.*, 57:59, 1984.</footnote>

<footnote n="28">²⁸H. Köppel, W. Domcke, and L. S. Cederbaum. The multi-mode vibronic-coupling approach. In W. Domcke, D. R. Yarkony, and H. Köppel, editors, *Conical Intersections*, pages 323–367, Singapore, 2004. World Scientific Co.</footnote>

<footnote n="29">²⁹Diptesh Dey, Joanne L Woodhouse, Marcus P Taylor, Helen H Fielding, and Graham A Worth. On the multiphoton ionisation photoelectron spectra of phenol. *Physical Chemistry Chemical Physics*, 26(4):3451–3461, 2024.</footnote>

<footnote n="30">³⁰H.-J. Werner, P. J. Knowles, G. Knizia, F. R. Manby, and M. Schütz. Molpro: a general-purpose quantum chemistry program package. *WIREs Comput Mol Sci*, 2:242–253, 2012.</footnote>

<footnote n="31">³¹H-D Meyer, Uwe Manthe, and Lorenz S Cederbaum. The multi-configurational time-dependent hartree approach. *Chemical Physics Letters*, 165(1):73–78, 1990.</footnote>

<footnote n="32">³²U. Manthe, H.-D. Meyer, and L. S. Cederbaum. Multiconfigurational time-dependent Hartree study of complex dynamics: Photodissociation of NO₂. *J. Chem. Phys.*, 97:9062–9071, 1992.</footnote>

<footnote n="33">³³M. H. Beck, A. Jäckle, G. A. Worth, and H.-D. Meyer. The multi-configuration time-dependent Hartree (MCTDH) method: A highly efficient algorithm for propagating wave packets. *Phys. Rep*, 324:1–105, 2000.</footnote>

<footnote n="34">³⁴Hans-Dieter Meyer, Fabien Gatti, and Graham A Worth. *Multidimensional quantum dynamics: MCTDH theory and applications*. John Wiley & Sons, 2009.</footnote>

<footnote n="35">³⁵Fabien Gatti, Benjamin Lasorne, Hans-Dieter Meyer, and André Nauts. *Applications of quantum dynamics in chemistry*, volume 98. Springer, 2017.</footnote>

<footnote n="36">³⁶David Mendive-Tapia, Hans-Dieter Meyer, and Oriol Vendrell. Optimal mode combination in the multiconfiguration time-dependent hartree method through multivariate statistics: Factor analysis and hierarchical clustering. *Journal of Chemical Theory and Computation*, 19(4):1144–1156, 2023.</footnote>
<footnote n="37"><sup>37</sup>D. Mendive-Tapia, T. Firmino, H.-D. Meyer, and F. Gatti. Towards a systematic conver-
gence of Multi-Layer (ML) multi-configuration time-dependent Hartree nuclear wavefunc-
tions: the ML-spawning algorithm. *Chem. Phys.*, 482:113–123, 2017.</footnote>

<footnote n="38"><sup>38</sup>D. Mendive-Tapia, E. Mangaud, T. Firmino, A. de la Lande, M. Desouter-Lecomte, H-D Meyer, and F. Gatti. Multidimensional quantum mechanical modeling of electron transfer and electronic coherence in plant cryptochromes: The role of initial bath conditions. *J. Phys. Chem. B*, 122:126–136, 2018.</footnote>

<footnote n="39"><sup>39</sup>G. A. Worth, M. H. Beck, A. Jäckle, and H.-D. Meyer. The MCTDH Package, Version 8.2, (2000). H.-D. Meyer, Version 8.3 (2002), Version 8.4 (2007), ML-MCTDH implemented in forthcoming Version 8.5 (2011). See http://mctdh.uni-hd.de/.</footnote>

<footnote n="40"><sup>40</sup>Gregg Jaeger. *Quantum information*. Springer, 2007.</footnote>

<footnote n="41"><sup>41</sup>John C Tully. Molecular dynamics with electronic transitions. *The Journal of Chemical Physics*, 93(2):1061–1071, 1990.</footnote>

<footnote n="42"><sup>42</sup>Giovanni Granucci and Maurizio Persico. Critical appraisal of the fewest switches algo-
rithm for surface hopping. *The Journal of chemical physics*, 126(13):134114, 2007.</footnote>

<footnote n="43"><sup>43</sup>Giovanni Granucci, Maurizio Persico, and Alberto Zoccante. Including quantum decoher-
ence in surface hopping. *The Journal of chemical physics*, 133(13):134111, 2010.</footnote>

<footnote n="44"><sup>44</sup>Chaoyuan Zhu, Shikha Nangia, Ahren W Jasper, and Donald G Truhlar. Coherent switch-
ing with decay of mixing: an improved treatment of electronic coherence for non-born–
oppenheimer trajectories. *The Journal of chemical physics*, 121(16):7658–7670, 2004.</footnote>

<footnote n="45"><sup>45</sup>Chaoyuan Zhu, Ahren W Jasper, and Donald G Truhlar. Non-born-oppenheimer liouville-
von neumann dynamics. evolution of a subsystem controlled by linear and population-
driven decay of mixing with decoherent and coherent switching. *Journal of chemical theory
and computation*, 1(4):527–540, 2005.</footnote>

<footnote n="46"><sup>46</sup>Lipeng Sun and William L Hase. Comparisons of classical and wigner sampling of transition state energy levels for quasiclassical trajectory chemical dynamics simulations. *The Journal of chemical physics*, 133(4):044313, 2010.</footnote>

<footnote n="47"><sup>47</sup>A. A. Granovsky. Firefly version 8, http://classic.chem.msu.su/gran/firefly/index.html, 2018.</footnote>

<footnote n="48"><sup>48</sup>Michael W Schmidt, Kim K Baldridge, Jerry A Boatz, Steven T Elbert, Mark S Gordon, Jan H Jensen, Shiro Koseki, Nikita Matsunaga, Kiet A Nguyen, Shujun Su, et al. General atomic and molecular electronic structure system. *Journal of computational chemistry*,</footnote>
14(11):1347-1363, 1993.

<footnote n="49">⁴⁹A I Kuleff, Jörg Breidbach, and L S Cederbaum. Multielectron wave-packet propagation: General theory and application. *The Journal of Chemical Physics*, 123:044111, 2005.</footnote>

<footnote n="50">⁵⁰Manuel Lara-Astiaso, Alicia Palacios, Piero Decleva, Ivano Tavernelli, and Fernando Martín. Role of electron-nuclear coupled dynamics on charge migration induced by attosecond pulses in glycine. *Chemical Physics Letters*, 683:357-364, 2017.</footnote>

<footnote n="51">⁵¹M Vacher, MJ Bearpark, and MA Robb. Oscillating charge migration between lone pairs persists without significant interaction with nuclear motion in the glycine and gly-gly-nh-ch3 radical cations. *The Journal of chemical physics*, 140:201102, 2014.</footnote>

<footnote n="52">⁵²P Cannington and N Ham. He(i) and he(i1) photoelectron spectra of glycine and related molecules. *J. Electron Spectrosc. Relat. Phenom.*, 139:139-151, 1983.</footnote>

<footnote n="53">⁵³Hannes Jónsson, Greg Mills, and Karsten W Jacobsen. Nudged elastic band method for finding minimum energy paths of transitions. 1998.</footnote>

<footnote n="54">⁵⁴Frank Neese. The orca program system. *Wiley Interdisciplinary Reviews: Computational Molecular Science*, 2(1):73-78, 2012.</footnote>

<footnote n="55">⁵⁵Frank Neese. Software update: the orca program system, version 4.0. *Wiley Interdisciplinary Reviews: Computational Molecular Science*, 8(1):e1327, 2018.</footnote>

<footnote n="56">⁵⁶MC Castrovilli, P Bolognesi, A Domaracka, P Rousseau, M C Castrovilli, S Richter, Chatterjee, F Wang, and L Avaldi. Insights into the dissociative ionization of glycine by pepico experiments. *Physical Chemistry Chemical Physics*, (20):22841-22848, 2018.</footnote>

<footnote n="57">⁵⁷David A Kreplin, Peter J Knowles, and Hans-Joachim Werner. Mcscf optimization revisited. ii. combined first-and second-order orbital optimization for large molecules. *The Journal of Chemical Physics*, 152(7), 2020.</footnote>

<footnote n="58">⁵⁸Thomas A Manz and Nidia Gabaldon Limas. Introducing ddec6 atomic population analysis: part 1. charge partitioning theory and methodology. *RSC advances*, 6(53):47771-47801, 2016.</footnote>